335 SUBROUTINE dgesvj( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
336 $ ldv, work, lwork, info )
344 INTEGER INFO, LDA, LDV, LWORK, M, MV, N
345 CHARACTER*1 JOBA, JOBU, JOBV
348 DOUBLE PRECISION A( lda, * ), SVA( n ), V( ldv, * ),
355 DOUBLE PRECISION ZERO, HALF, ONE
356 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0)
358 parameter( nsweep = 30 )
361 DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
362 $ bigtheta, cs, ctol, epsln, large, mxaapq,
363 $ mxsinj, rootbig, rooteps, rootsfmin, roottol,
364 $ skl, sfmin, small, sn, t, temp1, theta,
366 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
367 $ iswrot, jbc, jgl, kbl, lkahead, mvl, n2, n34,
368 $ n4, nbl, notrot, p, pskipped, q, rowskip,
370 LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
371 $ rsvec, uctol, upper
374 DOUBLE PRECISION FASTR( 5 )
377 INTRINSIC dabs, dmax1, dmin1, dble, min0, dsign, dsqrt
382 DOUBLE PRECISION DDOT, DNRM2
387 DOUBLE PRECISION DLAMCH
405 lsvec = lsame( jobu,
'U' )
406 uctol = lsame( jobu,
'C' )
407 rsvec = lsame( jobv,
'V' )
408 applv = lsame( jobv,
'A' )
409 upper = lsame( joba,
'U' )
410 lower = lsame( joba,
'L' )
412 IF( .NOT.( upper .OR. lower .OR. lsame( joba,
'G' ) ) )
THEN
414 ELSE IF( .NOT.( lsvec .OR. uctol .OR. lsame( jobu,
'N' ) ) )
THEN
416 ELSE IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv,
'N' ) ) )
THEN
418 ELSE IF( m.LT.0 )
THEN
420 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) )
THEN
422 ELSE IF( lda.LT.m )
THEN
424 ELSE IF( mv.LT.0 )
THEN
426 ELSE IF( ( rsvec .AND. ( ldv.LT.n ) ) .OR.
427 $ ( applv .AND. ( ldv.LT.mv ) ) )
THEN
429 ELSE IF( uctol .AND. ( work( 1 ).LE.one ) )
THEN
431 ELSE IF( lwork.LT.max0( m+n, 6 ) )
THEN
439 CALL xerbla(
'DGESVJ', -info )
445 IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) )
RETURN
459 IF( lsvec .OR. rsvec .OR. applv )
THEN
460 ctol = dsqrt( dble( m ) )
468 epsln = dlamch(
'Epsilon' )
469 rooteps = dsqrt( epsln )
470 sfmin = dlamch(
'SafeMinimum' )
471 rootsfmin = dsqrt( sfmin )
472 small = sfmin / epsln
473 big = dlamch(
'Overflow' )
475 rootbig = one / rootsfmin
476 large = big / dsqrt( dble( m*n ) )
477 bigtheta = one / rooteps
480 roottol = dsqrt( tol )
482 IF( dble( m )*epsln.GE.one )
THEN
484 CALL xerbla(
'DGESVJ', -info )
492 CALL dlaset(
'A', mvl, n, zero, one, v, ldv )
493 ELSE IF( applv )
THEN
496 rsvec = rsvec .OR. applv
507 skl= one / dsqrt( dble( m )*dble( n ) )
516 CALL dlassq( m-p+1, a( p, p ), 1, aapp, aaqq )
517 IF( aapp.GT.big )
THEN
519 CALL xerbla(
'DGESVJ', -info )
523 IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale )
THEN
527 sva( p ) = aapp*( aaqq*skl)
531 sva( q ) = sva( q )*skl
536 ELSE IF( upper )
THEN
541 CALL dlassq( p, a( 1, p ), 1, aapp, aaqq )
542 IF( aapp.GT.big )
THEN
544 CALL xerbla(
'DGESVJ', -info )
548 IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale )
THEN
552 sva( p ) = aapp*( aaqq*skl)
556 sva( q ) = sva( q )*skl
566 CALL dlassq( m, a( 1, p ), 1, aapp, aaqq )
567 IF( aapp.GT.big )
THEN
569 CALL xerbla(
'DGESVJ', -info )
573 IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale )
THEN
577 sva( p ) = aapp*( aaqq*skl)
581 sva( q ) = sva( q )*skl
588 IF( noscale )skl= one
597 IF( sva( p ).NE.zero )aaqq = dmin1( aaqq, sva( p ) )
598 aapp = dmax1( aapp, sva( p ) )
603 IF( aapp.EQ.zero )
THEN
604 IF( lsvec )
CALL dlaset(
'G', m, n, zero, one, a, lda )
617 IF( lsvec )
CALL dlascl(
'G', 0, 0, sva( 1 ), skl, m, 1,
618 $ a( 1, 1 ), lda, ierr )
619 work( 1 ) = one / skl
620 IF( sva( 1 ).GE.sfmin )
THEN
635 sn = dsqrt( sfmin / epsln )
636 temp1 = dsqrt( big / dble( n ) )
637 IF( ( aapp.LE.sn ) .OR. ( aaqq.GE.temp1 ) .OR.
638 $ ( ( sn.LE.aaqq ) .AND. ( aapp.LE.temp1 ) ) )
THEN
639 temp1 = dmin1( big, temp1 / aapp )
642 ELSE IF( ( aaqq.LE.sn ) .AND. ( aapp.LE.temp1 ) )
THEN
643 temp1 = dmin1( sn / aaqq, big / ( aapp*dsqrt( dble( n ) ) ) )
646 ELSE IF( ( aaqq.GE.sn ) .AND. ( aapp.GE.temp1 ) )
THEN
647 temp1 = dmax1( sn / aaqq, temp1 / aapp )
650 ELSE IF( ( aaqq.LE.sn ) .AND. ( aapp.GE.temp1 ) )
THEN
651 temp1 = dmin1( sn / aaqq, big / ( dsqrt( dble( n ) )*aapp ) )
660 IF( temp1.NE.one )
THEN
661 CALL dlascl(
'G', 0, 0, one, temp1, n, 1, sva, n, ierr )
664 IF( skl.NE.one )
THEN
665 CALL dlascl( joba, 0, 0, one, skl, m, n, a, lda, ierr )
671 emptsw = ( n*( n-1 ) ) / 2
699 IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
704 rowskip = min0( 5, kbl )
715 IF( ( lower .OR. upper ) .AND. ( n.GT.max0( 64, 4*kbl ) ) )
THEN
737 CALL dgsvj0( jobv, m-n34, n-n34, a( n34+1, n34+1 ), lda,
738 $ work( n34+1 ), sva( n34+1 ), mvl,
739 $ v( n34*q+1, n34+1 ), ldv, epsln, sfmin, tol,
740 $ 2, work( n+1 ), lwork-n, ierr )
742 CALL dgsvj0( jobv, m-n2, n34-n2, a( n2+1, n2+1 ), lda,
743 $ work( n2+1 ), sva( n2+1 ), mvl,
744 $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 2,
745 $ work( n+1 ), lwork-n, ierr )
747 CALL dgsvj1( jobv, m-n2, n-n2, n4, a( n2+1, n2+1 ), lda,
748 $ work( n2+1 ), sva( n2+1 ), mvl,
749 $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 1,
750 $ work( n+1 ), lwork-n, ierr )
752 CALL dgsvj0( jobv, m-n4, n2-n4, a( n4+1, n4+1 ), lda,
753 $ work( n4+1 ), sva( n4+1 ), mvl,
754 $ v( n4*q+1, n4+1 ), ldv, epsln, sfmin, tol, 1,
755 $ work( n+1 ), lwork-n, ierr )
757 CALL dgsvj0( jobv, m, n4, a, lda, work, sva, mvl, v, ldv,
758 $ epsln, sfmin, tol, 1, work( n+1 ), lwork-n,
761 CALL dgsvj1( jobv, m, n2, n4, a, lda, work, sva, mvl, v,
762 $ ldv, epsln, sfmin, tol, 1, work( n+1 ),
766 ELSE IF( upper )
THEN
769 CALL dgsvj0( jobv, n4, n4, a, lda, work, sva, mvl, v, ldv,
770 $ epsln, sfmin, tol, 2, work( n+1 ), lwork-n,
773 CALL dgsvj0( jobv, n2, n4, a( 1, n4+1 ), lda, work( n4+1 ),
774 $ sva( n4+1 ), mvl, v( n4*q+1, n4+1 ), ldv,
775 $ epsln, sfmin, tol, 1, work( n+1 ), lwork-n,
778 CALL dgsvj1( jobv, n2, n2, n4, a, lda, work, sva, mvl, v,
779 $ ldv, epsln, sfmin, tol, 1, work( n+1 ),
782 CALL dgsvj0( jobv, n2+n4, n4, a( 1, n2+1 ), lda,
783 $ work( n2+1 ), sva( n2+1 ), mvl,
784 $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 1,
785 $ work( n+1 ), lwork-n, ierr )
793 DO 1993 i = 1, nsweep
811 igl = ( ibr-1 )*kbl + 1
813 DO 1002 ir1 = 0, min0( lkahead, nbl-ibr )
817 DO 2001 p = igl, min0( igl+kbl-1, n-1 )
821 q = idamax( n-p+1, sva( p ), 1 ) + p - 1
823 CALL dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
824 IF( rsvec )
CALL dswap( mvl, v( 1, p ), 1,
830 work( p ) = work( q )
848 IF( ( sva( p ).LT.rootbig ) .AND.
849 $ ( sva( p ).GT.rootsfmin ) )
THEN
850 sva( p ) = dnrm2( m, a( 1, p ), 1 )*work( p )
854 CALL dlassq( m, a( 1, p ), 1, temp1, aapp )
855 sva( p ) = temp1*dsqrt( aapp )*work( p )
862 IF( aapp.GT.zero )
THEN
866 DO 2002 q = p + 1, min0( igl+kbl-1, n )
870 IF( aaqq.GT.zero )
THEN
873 IF( aaqq.GE.one )
THEN
874 rotok = ( small*aapp ).LE.aaqq
875 IF( aapp.LT.( big / aaqq ) )
THEN
876 aapq = ( ddot( m, a( 1, p ), 1, a( 1,
877 $ q ), 1 )*work( p )*work( q ) /
880 CALL dcopy( m, a( 1, p ), 1,
882 CALL dlascl(
'G', 0, 0, aapp,
884 $ work( n+1 ), lda, ierr )
885 aapq = ddot( m, work( n+1 ), 1,
886 $ a( 1, q ), 1 )*work( q ) / aaqq
889 rotok = aapp.LE.( aaqq / small )
890 IF( aapp.GT.( small / aaqq ) )
THEN
891 aapq = ( ddot( m, a( 1, p ), 1, a( 1,
892 $ q ), 1 )*work( p )*work( q ) /
895 CALL dcopy( m, a( 1, q ), 1,
897 CALL dlascl(
'G', 0, 0, aaqq,
899 $ work( n+1 ), lda, ierr )
900 aapq = ddot( m, work( n+1 ), 1,
901 $ a( 1, p ), 1 )*work( p ) / aapp
905 mxaapq = dmax1( mxaapq, dabs( aapq ) )
909 IF( dabs( aapq ).GT.tol )
THEN
924 theta = -half*dabs(aqoap-apoaq)/aapq
926 IF( dabs( theta ).GT.bigtheta )
THEN
929 fastr( 3 ) = t*work( p ) / work( q )
930 fastr( 4 ) = -t*work( q ) /
932 CALL drotm( m, a( 1, p ), 1,
933 $ a( 1, q ), 1, fastr )
934 IF( rsvec )
CALL drotm( mvl,
938 sva( q ) = aaqq*dsqrt( dmax1( zero,
939 $ one+t*apoaq*aapq ) )
940 aapp = aapp*dsqrt( dmax1( zero,
941 $ one-t*aqoap*aapq ) )
942 mxsinj = dmax1( mxsinj, dabs( t ) )
948 thsign = -dsign( one, aapq )
949 t = one / ( theta+thsign*
950 $ dsqrt( one+theta*theta ) )
951 cs = dsqrt( one / ( one+t*t ) )
954 mxsinj = dmax1( mxsinj, dabs( sn ) )
955 sva( q ) = aaqq*dsqrt( dmax1( zero,
956 $ one+t*apoaq*aapq ) )
957 aapp = aapp*dsqrt( dmax1( zero,
958 $ one-t*aqoap*aapq ) )
960 apoaq = work( p ) / work( q )
961 aqoap = work( q ) / work( p )
962 IF( work( p ).GE.one )
THEN
963 IF( work( q ).GE.one )
THEN
965 fastr( 4 ) = -t*aqoap
966 work( p ) = work( p )*cs
967 work( q ) = work( q )*cs
968 CALL drotm( m, a( 1, p ), 1,
971 IF( rsvec )
CALL drotm( mvl,
972 $ v( 1, p ), 1, v( 1, q ),
975 CALL daxpy( m, -t*aqoap,
978 CALL daxpy( m, cs*sn*apoaq,
981 work( p ) = work( p )*cs
982 work( q ) = work( q ) / cs
984 CALL daxpy( mvl, -t*aqoap,
994 IF( work( q ).GE.one )
THEN
995 CALL daxpy( m, t*apoaq,
998 CALL daxpy( m, -cs*sn*aqoap,
1001 work( p ) = work( p ) / cs
1002 work( q ) = work( q )*cs
1004 CALL daxpy( mvl, t*apoaq,
1013 IF( work( p ).GE.work( q ) )
1015 CALL daxpy( m, -t*aqoap,
1018 CALL daxpy( m, cs*sn*apoaq,
1021 work( p ) = work( p )*cs
1022 work( q ) = work( q ) / cs
1034 CALL daxpy( m, t*apoaq,
1041 work( p ) = work( p ) / cs
1042 work( q ) = work( q )*cs
1045 $ t*apoaq, v( 1, p ),
1059 CALL dcopy( m, a( 1, p ), 1,
1061 CALL dlascl(
'G', 0, 0, aapp, one, m,
1062 $ 1, work( n+1 ), lda,
1064 CALL dlascl(
'G', 0, 0, aaqq, one, m,
1065 $ 1, a( 1, q ), lda, ierr )
1066 temp1 = -aapq*work( p ) / work( q )
1067 CALL daxpy( m, temp1, work( n+1 ), 1,
1069 CALL dlascl(
'G', 0, 0, one, aaqq, m,
1070 $ 1, a( 1, q ), lda, ierr )
1071 sva( q ) = aaqq*dsqrt( dmax1( zero,
1073 mxsinj = dmax1( mxsinj, sfmin )
1080 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
1082 IF( ( aaqq.LT.rootbig ) .AND.
1083 $ ( aaqq.GT.rootsfmin ) )
THEN
1084 sva( q ) = dnrm2( m, a( 1, q ), 1 )*
1089 CALL dlassq( m, a( 1, q ), 1, t,
1091 sva( q ) = t*dsqrt( aaqq )*work( q )
1094 IF( ( aapp / aapp0 ).LE.rooteps )
THEN
1095 IF( ( aapp.LT.rootbig ) .AND.
1096 $ ( aapp.GT.rootsfmin ) )
THEN
1097 aapp = dnrm2( m, a( 1, p ), 1 )*
1102 CALL dlassq( m, a( 1, p ), 1, t,
1104 aapp = t*dsqrt( aapp )*work( p )
1111 IF( ir1.EQ.0 )notrot = notrot + 1
1113 pskipped = pskipped + 1
1117 IF( ir1.EQ.0 )notrot = notrot + 1
1118 pskipped = pskipped + 1
1121 IF( ( i.LE.swband ) .AND.
1122 $ ( pskipped.GT.rowskip ) )
THEN
1123 IF( ir1.EQ.0 )aapp = -aapp
1138 IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
1139 $ notrot = notrot + min0( igl+kbl-1, n ) - p
1150 igl = ( ibr-1 )*kbl + 1
1152 DO 2010 jbc = ibr + 1, nbl
1154 jgl = ( jbc-1 )*kbl + 1
1159 DO 2100 p = igl, min0( igl+kbl-1, n )
1162 IF( aapp.GT.zero )
THEN
1166 DO 2200 q = jgl, min0( jgl+kbl-1, n )
1169 IF( aaqq.GT.zero )
THEN
1176 IF( aaqq.GE.one )
THEN
1177 IF( aapp.GE.aaqq )
THEN
1178 rotok = ( small*aapp ).LE.aaqq
1180 rotok = ( small*aaqq ).LE.aapp
1182 IF( aapp.LT.( big / aaqq ) )
THEN
1183 aapq = ( ddot( m, a( 1, p ), 1, a( 1,
1184 $ q ), 1 )*work( p )*work( q ) /
1187 CALL dcopy( m, a( 1, p ), 1,
1189 CALL dlascl(
'G', 0, 0, aapp,
1191 $ work( n+1 ), lda, ierr )
1192 aapq = ddot( m, work( n+1 ), 1,
1193 $ a( 1, q ), 1 )*work( q ) / aaqq
1196 IF( aapp.GE.aaqq )
THEN
1197 rotok = aapp.LE.( aaqq / small )
1199 rotok = aaqq.LE.( aapp / small )
1201 IF( aapp.GT.( small / aaqq ) )
THEN
1202 aapq = ( ddot( m, a( 1, p ), 1, a( 1,
1203 $ q ), 1 )*work( p )*work( q ) /
1206 CALL dcopy( m, a( 1, q ), 1,
1208 CALL dlascl(
'G', 0, 0, aaqq,
1210 $ work( n+1 ), lda, ierr )
1211 aapq = ddot( m, work( n+1 ), 1,
1212 $ a( 1, p ), 1 )*work( p ) / aapp
1216 mxaapq = dmax1( mxaapq, dabs( aapq ) )
1220 IF( dabs( aapq ).GT.tol )
THEN
1230 theta = -half*dabs(aqoap-apoaq)/aapq
1231 IF( aaqq.GT.aapp0 )theta = -theta
1233 IF( dabs( theta ).GT.bigtheta )
THEN
1235 fastr( 3 ) = t*work( p ) / work( q )
1236 fastr( 4 ) = -t*work( q ) /
1238 CALL drotm( m, a( 1, p ), 1,
1239 $ a( 1, q ), 1, fastr )
1240 IF( rsvec )
CALL drotm( mvl,
1244 sva( q ) = aaqq*dsqrt( dmax1( zero,
1245 $ one+t*apoaq*aapq ) )
1246 aapp = aapp*dsqrt( dmax1( zero,
1247 $ one-t*aqoap*aapq ) )
1248 mxsinj = dmax1( mxsinj, dabs( t ) )
1253 thsign = -dsign( one, aapq )
1254 IF( aaqq.GT.aapp0 )thsign = -thsign
1255 t = one / ( theta+thsign*
1256 $ dsqrt( one+theta*theta ) )
1257 cs = dsqrt( one / ( one+t*t ) )
1259 mxsinj = dmax1( mxsinj, dabs( sn ) )
1260 sva( q ) = aaqq*dsqrt( dmax1( zero,
1261 $ one+t*apoaq*aapq ) )
1262 aapp = aapp*dsqrt( dmax1( zero,
1263 $ one-t*aqoap*aapq ) )
1265 apoaq = work( p ) / work( q )
1266 aqoap = work( q ) / work( p )
1267 IF( work( p ).GE.one )
THEN
1269 IF( work( q ).GE.one )
THEN
1270 fastr( 3 ) = t*apoaq
1271 fastr( 4 ) = -t*aqoap
1272 work( p ) = work( p )*cs
1273 work( q ) = work( q )*cs
1274 CALL drotm( m, a( 1, p ), 1,
1277 IF( rsvec )
CALL drotm( mvl,
1278 $ v( 1, p ), 1, v( 1, q ),
1281 CALL daxpy( m, -t*aqoap,
1284 CALL daxpy( m, cs*sn*apoaq,
1288 CALL daxpy( mvl, -t*aqoap,
1296 work( p ) = work( p )*cs
1297 work( q ) = work( q ) / cs
1300 IF( work( q ).GE.one )
THEN
1301 CALL daxpy( m, t*apoaq,
1304 CALL daxpy( m, -cs*sn*aqoap,
1308 CALL daxpy( mvl, t*apoaq,
1316 work( p ) = work( p ) / cs
1317 work( q ) = work( q )*cs
1319 IF( work( p ).GE.work( q ) )
1321 CALL daxpy( m, -t*aqoap,
1324 CALL daxpy( m, cs*sn*apoaq,
1327 work( p ) = work( p )*cs
1328 work( q ) = work( q ) / cs
1340 CALL daxpy( m, t*apoaq,
1347 work( p ) = work( p ) / cs
1348 work( q ) = work( q )*cs
1351 $ t*apoaq, v( 1, p ),
1364 IF( aapp.GT.aaqq )
THEN
1365 CALL dcopy( m, a( 1, p ), 1,
1367 CALL dlascl(
'G', 0, 0, aapp, one,
1368 $ m, 1, work( n+1 ), lda,
1370 CALL dlascl(
'G', 0, 0, aaqq, one,
1371 $ m, 1, a( 1, q ), lda,
1373 temp1 = -aapq*work( p ) / work( q )
1374 CALL daxpy( m, temp1, work( n+1 ),
1376 CALL dlascl(
'G', 0, 0, one, aaqq,
1377 $ m, 1, a( 1, q ), lda,
1379 sva( q ) = aaqq*dsqrt( dmax1( zero,
1381 mxsinj = dmax1( mxsinj, sfmin )
1383 CALL dcopy( m, a( 1, q ), 1,
1385 CALL dlascl(
'G', 0, 0, aaqq, one,
1386 $ m, 1, work( n+1 ), lda,
1388 CALL dlascl(
'G', 0, 0, aapp, one,
1389 $ m, 1, a( 1, p ), lda,
1391 temp1 = -aapq*work( q ) / work( p )
1392 CALL daxpy( m, temp1, work( n+1 ),
1394 CALL dlascl(
'G', 0, 0, one, aapp,
1395 $ m, 1, a( 1, p ), lda,
1397 sva( p ) = aapp*dsqrt( dmax1( zero,
1399 mxsinj = dmax1( mxsinj, sfmin )
1406 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
1408 IF( ( aaqq.LT.rootbig ) .AND.
1409 $ ( aaqq.GT.rootsfmin ) )
THEN
1410 sva( q ) = dnrm2( m, a( 1, q ), 1 )*
1415 CALL dlassq( m, a( 1, q ), 1, t,
1417 sva( q ) = t*dsqrt( aaqq )*work( q )
1420 IF( ( aapp / aapp0 )**2.LE.rooteps )
THEN
1421 IF( ( aapp.LT.rootbig ) .AND.
1422 $ ( aapp.GT.rootsfmin ) )
THEN
1423 aapp = dnrm2( m, a( 1, p ), 1 )*
1428 CALL dlassq( m, a( 1, p ), 1, t,
1430 aapp = t*dsqrt( aapp )*work( p )
1438 pskipped = pskipped + 1
1443 pskipped = pskipped + 1
1447 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
1453 IF( ( i.LE.swband ) .AND.
1454 $ ( pskipped.GT.rowskip ) )
THEN
1468 IF( aapp.EQ.zero )notrot = notrot +
1469 $ min0( jgl+kbl-1, n ) - jgl + 1
1470 IF( aapp.LT.zero )notrot = 0
1480 DO 2012 p = igl, min0( igl+kbl-1, n )
1481 sva( p ) = dabs( sva( p ) )
1488 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
1490 sva( n ) = dnrm2( m, a( 1, n ), 1 )*work( n )
1494 CALL dlassq( m, a( 1, n ), 1, t, aapp )
1495 sva( n ) = t*dsqrt( aapp )*work( n )
1500 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
1501 $ ( iswrot.LE.n ) ) )swband = i
1503 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.dsqrt( dble( n ) )*
1504 $ tol ) .AND. ( dble( n )*mxaapq*mxsinj.LT.tol ) )
THEN
1508 IF( notrot.GE.emptsw )
GO TO 1994
1530 DO 5991 p = 1, n - 1
1531 q = idamax( n-p+1, sva( p ), 1 ) + p - 1
1537 work( p ) = work( q )
1539 CALL dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
1540 IF( rsvec )
CALL dswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
1542 IF( sva( p ).NE.zero )
THEN
1544 IF( sva( p )*skl.GT.sfmin )n2 = n2 + 1
1547 IF( sva( n ).NE.zero )
THEN
1549 IF( sva( n )*skl.GT.sfmin )n2 = n2 + 1
1554 IF( lsvec .OR. uctol )
THEN
1556 CALL dscal( m, work( p ) / sva( p ), a( 1, p ), 1 )
1565 CALL dscal( mvl, work( p ), v( 1, p ), 1 )
1569 temp1 = one / dnrm2( mvl, v( 1, p ), 1 )
1570 CALL dscal( mvl, temp1, v( 1, p ), 1 )
1576 IF( ( ( skl.GT.one ) .AND. ( sva( 1 ).LT.( big / skl) ) )
1577 $ .OR. ( ( skl.LT.one ) .AND. ( sva( max( n2, 1 ) ) .GT.
1578 $ ( sfmin / skl) ) ) )
THEN
1580 sva( p ) = skl*sva( p )
1590 work( 2 ) = dble( n4 )
1593 work( 3 ) = dble( n2 )
1598 work( 4 ) = dble( i )
subroutine drotm(N, DX, INCX, DY, INCY, DPARAM)
DROTM
subroutine dgsvj0(JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
DGSVJ0 pre-processor for the routine sgesvj.
subroutine dgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
DGESVJ
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 dscal(N, DA, DX, INCX)
DSCAL
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
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
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...