236 SUBROUTINE dgsvj1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
237 $ eps, sfmin, tol, nsweep, work, lwork, info )
245 DOUBLE PRECISION eps, sfmin, tol
246 INTEGER info, lda, ldv, lwork, m, mv, n, n1, nsweep
250 DOUBLE PRECISION a( lda, * ), d( n ), sva( n ), v( ldv, * ),
257 DOUBLE PRECISION zero, half, one
258 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
261 DOUBLE PRECISION aapp, aapp0, aapq, aaqq, apoaq, aqoap, big,
262 $ bigtheta, cs, large, mxaapq, mxsinj, rootbig,
263 $ rooteps, rootsfmin, roottol, small, sn, t,
264 $ temp1, theta, thsign
265 INTEGER blskip, emptsw, i, ibr, igl, ierr, ijblsk,
266 $ iswrot, jbc, jgl, kbl, mvl, notrot, nblc, nblr,
267 $ p, pskipped, q, rowskip, swband
268 LOGICAL applv, rotok, rsvec
271 DOUBLE PRECISION fastr( 5 )
274 INTRINSIC dabs, dmax1, dble, min0, dsign, dsqrt
289 applv =
lsame( jobv,
'A' )
290 rsvec =
lsame( jobv,
'V' )
291 IF( .NOT.( rsvec .OR. applv .OR.
lsame( jobv,
'N' ) ) )
THEN
293 ELSE IF( m.LT.0 )
THEN
295 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) )
THEN
297 ELSE IF( n1.LT.0 )
THEN
299 ELSE IF( lda.LT.m )
THEN
301 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) )
THEN
303 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
304 $ ( applv.AND.( ldv.LT.mv ) ) )
THEN
306 ELSE IF( tol.LE.eps )
THEN
308 ELSE IF( nsweep.LT.0 )
THEN
310 ELSE IF( lwork.LT.m )
THEN
318 CALL
xerbla(
'DGSVJ1', -info )
324 ELSE IF( applv )
THEN
327 rsvec = rsvec .OR. applv
329 rooteps = dsqrt( eps )
330 rootsfmin = dsqrt( sfmin )
333 rootbig = one / rootsfmin
334 large = big / dsqrt( dble( m*n ) )
335 bigtheta = one / rooteps
336 roottol = dsqrt( tol )
350 IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
354 nblc = ( n-n1 ) / kbl
355 IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
356 blskip = ( kbl**2 ) + 1
359 rowskip = min0( 5, kbl )
375 DO 1993 i = 1, nsweep
385 DO 2000 ibr = 1, nblr
387 igl = ( ibr-1 )*kbl + 1
393 igl = ( ibr-1 )*kbl + 1
395 DO 2010 jbc = 1, nblc
397 jgl = n1 + ( jbc-1 )*kbl + 1
402 DO 2100 p = igl, min0( igl+kbl-1, n1 )
406 IF( aapp.GT.zero )
THEN
410 DO 2200 q = jgl, min0( jgl+kbl-1, n )
414 IF( aaqq.GT.zero )
THEN
421 IF( aaqq.GE.one )
THEN
422 IF( aapp.GE.aaqq )
THEN
423 rotok = ( small*aapp ).LE.aaqq
425 rotok = ( small*aaqq ).LE.aapp
427 IF( aapp.LT.( big / aaqq ) )
THEN
428 aapq = (
ddot( m, a( 1, p ), 1, a( 1,
429 $ q ), 1 )*d( p )*d( q ) / aaqq )
432 CALL
dcopy( m, a( 1, p ), 1, work, 1 )
433 CALL
dlascl(
'G', 0, 0, aapp, d( p ),
434 $ m, 1, work, lda, ierr )
435 aapq =
ddot( m, work, 1, a( 1, q ),
439 IF( aapp.GE.aaqq )
THEN
440 rotok = aapp.LE.( aaqq / small )
442 rotok = aaqq.LE.( aapp / small )
444 IF( aapp.GT.( small / aaqq ) )
THEN
445 aapq = (
ddot( m, a( 1, p ), 1, a( 1,
446 $ q ), 1 )*d( p )*d( q ) / aaqq )
449 CALL
dcopy( m, a( 1, q ), 1, work, 1 )
450 CALL
dlascl(
'G', 0, 0, aaqq, d( q ),
451 $ m, 1, work, lda, ierr )
452 aapq =
ddot( m, work, 1, a( 1, p ),
457 mxaapq = dmax1( mxaapq, dabs( aapq ) )
461 IF( dabs( aapq ).GT.tol )
THEN
471 theta = -half*dabs(aqoap-apoaq) / aapq
472 IF( aaqq.GT.aapp0 )theta = -theta
474 IF( dabs( theta ).GT.bigtheta )
THEN
476 fastr( 3 ) = t*d( p ) / d( q )
477 fastr( 4 ) = -t*d( q ) / d( p )
478 CALL
drotm( m, a( 1, p ), 1,
479 $ a( 1, q ), 1, fastr )
480 IF( rsvec )CALL
drotm( mvl,
484 sva( q ) = aaqq*dsqrt( dmax1( zero,
485 $ one+t*apoaq*aapq ) )
486 aapp = aapp*dsqrt( dmax1( zero,
487 $ one-t*aqoap*aapq ) )
488 mxsinj = dmax1( mxsinj, dabs( t ) )
493 thsign = -dsign( one, aapq )
494 IF( aaqq.GT.aapp0 )thsign = -thsign
495 t = one / ( theta+thsign*
496 $ dsqrt( one+theta*theta ) )
497 cs = dsqrt( one / ( one+t*t ) )
499 mxsinj = dmax1( mxsinj, dabs( sn ) )
500 sva( q ) = aaqq*dsqrt( dmax1( zero,
501 $ one+t*apoaq*aapq ) )
502 aapp = aapp*dsqrt( dmax1( zero,
503 $ one-t*aqoap*aapq ) )
505 apoaq = d( p ) / d( q )
506 aqoap = d( q ) / d( p )
507 IF( d( p ).GE.one )
THEN
509 IF( d( q ).GE.one )
THEN
511 fastr( 4 ) = -t*aqoap
514 CALL
drotm( m, a( 1, p ), 1,
517 IF( rsvec )CALL
drotm( mvl,
518 $ v( 1, p ), 1, v( 1, q ),
521 CALL
daxpy( m, -t*aqoap,
524 CALL
daxpy( m, cs*sn*apoaq,
528 CALL
daxpy( mvl, -t*aqoap,
540 IF( d( q ).GE.one )
THEN
541 CALL
daxpy( m, t*apoaq,
544 CALL
daxpy( m, -cs*sn*aqoap,
548 CALL
daxpy( mvl, t*apoaq,
559 IF( d( p ).GE.d( q ) )
THEN
560 CALL
daxpy( m, -t*aqoap,
563 CALL
daxpy( m, cs*sn*apoaq,
579 CALL
daxpy( m, t*apoaq,
590 $ t*apoaq, v( 1, p ),
603 IF( aapp.GT.aaqq )
THEN
604 CALL
dcopy( m, a( 1, p ), 1, work,
606 CALL
dlascl(
'G', 0, 0, aapp, one,
607 $ m, 1, work, lda, ierr )
608 CALL
dlascl(
'G', 0, 0, aaqq, one,
609 $ m, 1, a( 1, q ), lda,
611 temp1 = -aapq*d( p ) / d( q )
612 CALL
daxpy( m, temp1, work, 1,
614 CALL
dlascl(
'G', 0, 0, one, aaqq,
615 $ m, 1, a( 1, q ), lda,
617 sva( q ) = aaqq*dsqrt( dmax1( zero,
619 mxsinj = dmax1( mxsinj, sfmin )
621 CALL
dcopy( m, a( 1, q ), 1, work,
623 CALL
dlascl(
'G', 0, 0, aaqq, one,
624 $ m, 1, work, lda, ierr )
625 CALL
dlascl(
'G', 0, 0, aapp, one,
626 $ m, 1, a( 1, p ), lda,
628 temp1 = -aapq*d( q ) / d( p )
629 CALL
daxpy( m, temp1, work, 1,
631 CALL
dlascl(
'G', 0, 0, one, aapp,
632 $ m, 1, a( 1, p ), lda,
634 sva( p ) = aapp*dsqrt( dmax1( zero,
636 mxsinj = dmax1( mxsinj, sfmin )
643 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
645 IF( ( aaqq.LT.rootbig ) .AND.
646 $ ( aaqq.GT.rootsfmin ) )
THEN
647 sva( q ) =
dnrm2( m, a( 1, q ), 1 )*
652 CALL
dlassq( m, a( 1, q ), 1, t,
654 sva( q ) = t*dsqrt( aaqq )*d( q )
657 IF( ( aapp / aapp0 )**2.LE.rooteps )
THEN
658 IF( ( aapp.LT.rootbig ) .AND.
659 $ ( aapp.GT.rootsfmin ) )
THEN
660 aapp =
dnrm2( m, a( 1, p ), 1 )*
665 CALL
dlassq( m, a( 1, p ), 1, t,
667 aapp = t*dsqrt( aapp )*d( p )
675 pskipped = pskipped + 1
680 pskipped = pskipped + 1
685 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
691 IF( ( i.LE.swband ) .AND.
692 $ ( pskipped.GT.rowskip ) )
THEN
706 IF( aapp.EQ.zero )notrot = notrot +
707 $ min0( jgl+kbl-1, n ) - jgl + 1
708 IF( aapp.LT.zero )notrot = 0
718 DO 2012 p = igl, min0( igl+kbl-1, n )
719 sva( p ) = dabs( sva( p ) )
726 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
728 sva( n ) =
dnrm2( m, a( 1, n ), 1 )*d( n )
732 CALL
dlassq( m, a( 1, n ), 1, t, aapp )
733 sva( n ) = t*dsqrt( aapp )*d( n )
738 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
739 $ ( iswrot.LE.n ) ) )swband = i
741 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.dble( n )*tol ) .AND.
742 $ ( dble( n )*mxaapq*mxsinj.LT.tol ) )
THEN
747 IF( notrot.GE.emptsw )go to 1994
766 q =
idamax( n-p+1, sva( p ), 1 ) + p - 1
774 CALL
dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
775 IF( rsvec )CALL
dswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
subroutine drotm(N, DX, INCX, DY, INCY, DPARAM)
DROTM
LOGICAL function lsame(CA, CB)
LSAME
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
DOUBLE PRECISION function ddot(N, DX, INCX, DY, INCY)
DDOT
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dgsvj1(JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
DGSVJ1 pre-processor for the routine sgesvj, applies Jacobi rotations targeting only particular pivot...
DOUBLE PRECISION function dnrm2(N, X, INCX)
DNRM2
INTEGER function idamax(N, DX, INCX)
IDAMAX