321 SUBROUTINE sgesvj( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
322 $ ldv, work, lwork, info )
330 INTEGER INFO, LDA, LDV, LWORK, M, MV, N
331 CHARACTER*1 JOBA, JOBU, JOBV
334 REAL A( lda, * ), SVA( n ), V( ldv, * ),
342 parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
344 parameter( nsweep = 30 )
347 REAL AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
348 $ bigtheta, cs, ctol, epsln, large, mxaapq,
349 $ mxsinj, rootbig, rooteps, rootsfmin, roottol,
350 $ skl, sfmin, small, sn, t, temp1, theta,
352 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
353 $ iswrot, jbc, jgl, kbl, lkahead, mvl, n2, n34,
354 $ n4, nbl, notrot, p, pskipped, q, rowskip,
356 LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
357 $ rsvec, uctol, upper
363 INTRINSIC abs, amax1, amin1, float, min0, sign, sqrt
391 lsvec = lsame( jobu,
'U' )
392 uctol = lsame( jobu,
'C' )
393 rsvec = lsame( jobv,
'V' )
394 applv = lsame( jobv,
'A' )
395 upper = lsame( joba,
'U' )
396 lower = lsame( joba,
'L' )
398 IF( .NOT.( upper .OR. lower .OR. lsame( joba,
'G' ) ) )
THEN
400 ELSE IF( .NOT.( lsvec .OR. uctol .OR. lsame( jobu,
'N' ) ) )
THEN
402 ELSE IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv,
'N' ) ) )
THEN
404 ELSE IF( m.LT.0 )
THEN
406 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) )
THEN
408 ELSE IF( lda.LT.m )
THEN
410 ELSE IF( mv.LT.0 )
THEN
412 ELSE IF( ( rsvec .AND. ( ldv.LT.n ) ) .OR.
413 $ ( applv .AND. ( ldv.LT.mv ) ) )
THEN
415 ELSE IF( uctol .AND. ( work( 1 ).LE.one ) )
THEN
417 ELSE IF( lwork.LT.max0( m+n, 6 ) )
THEN
425 CALL xerbla(
'SGESVJ', -info )
431 IF( ( m.EQ.0 ) .OR. ( n.EQ.0 ) )
RETURN
445 IF( lsvec .OR. rsvec .OR. applv )
THEN
446 ctol = sqrt( float( m ) )
454 epsln = slamch(
'Epsilon' )
455 rooteps = sqrt( epsln )
456 sfmin = slamch(
'SafeMinimum' )
457 rootsfmin = sqrt( sfmin )
458 small = sfmin / epsln
459 big = slamch(
'Overflow' )
461 rootbig = one / rootsfmin
462 large = big / sqrt( float( m*n ) )
463 bigtheta = one / rooteps
466 roottol = sqrt( tol )
468 IF( float( m )*epsln.GE.one )
THEN
470 CALL xerbla(
'SGESVJ', -info )
478 CALL slaset(
'A', mvl, n, zero, one, v, ldv )
479 ELSE IF( applv )
THEN
482 rsvec = rsvec .OR. applv
493 skl = one / sqrt( float( m )*float( n ) )
502 CALL slassq( m-p+1, a( p, p ), 1, aapp, aaqq )
503 IF( aapp.GT.big )
THEN
505 CALL xerbla(
'SGESVJ', -info )
509 IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale )
THEN
513 sva( p ) = aapp*( aaqq*skl )
517 sva( q ) = sva( q )*skl
522 ELSE IF( upper )
THEN
527 CALL slassq( p, a( 1, p ), 1, aapp, aaqq )
528 IF( aapp.GT.big )
THEN
530 CALL xerbla(
'SGESVJ', -info )
534 IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale )
THEN
538 sva( p ) = aapp*( aaqq*skl )
542 sva( q ) = sva( q )*skl
552 CALL slassq( m, a( 1, p ), 1, aapp, aaqq )
553 IF( aapp.GT.big )
THEN
555 CALL xerbla(
'SGESVJ', -info )
559 IF( ( aapp.LT.( big / aaqq ) ) .AND. noscale )
THEN
563 sva( p ) = aapp*( aaqq*skl )
567 sva( q ) = sva( q )*skl
574 IF( noscale )skl = one
583 IF( sva( p ).NE.zero )aaqq = amin1( aaqq, sva( p ) )
584 aapp = amax1( aapp, sva( p ) )
589 IF( aapp.EQ.zero )
THEN
590 IF( lsvec )
CALL slaset(
'G', m, n, zero, one, a, lda )
603 IF( lsvec )
CALL slascl(
'G', 0, 0, sva( 1 ), skl, m, 1,
604 $ a( 1, 1 ), lda, ierr )
605 work( 1 ) = one / skl
606 IF( sva( 1 ).GE.sfmin )
THEN
621 sn = sqrt( sfmin / epsln )
622 temp1 = sqrt( big / float( n ) )
623 IF( ( aapp.LE.sn ) .OR. ( aaqq.GE.temp1 ) .OR.
624 $ ( ( sn.LE.aaqq ) .AND. ( aapp.LE.temp1 ) ) )
THEN
625 temp1 = amin1( big, temp1 / aapp )
628 ELSE IF( ( aaqq.LE.sn ) .AND. ( aapp.LE.temp1 ) )
THEN
629 temp1 = amin1( sn / aaqq, big / ( aapp*sqrt( float( n ) ) ) )
632 ELSE IF( ( aaqq.GE.sn ) .AND. ( aapp.GE.temp1 ) )
THEN
633 temp1 = amax1( sn / aaqq, temp1 / aapp )
636 ELSE IF( ( aaqq.LE.sn ) .AND. ( aapp.GE.temp1 ) )
THEN
637 temp1 = amin1( sn / aaqq, big / ( sqrt( float( n ) )*aapp ) )
646 IF( temp1.NE.one )
THEN
647 CALL slascl(
'G', 0, 0, one, temp1, n, 1, sva, n, ierr )
650 IF( skl.NE.one )
THEN
651 CALL slascl( joba, 0, 0, one, skl, m, n, a, lda, ierr )
657 emptsw = ( n*( n-1 ) ) / 2
685 IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
690 rowskip = min0( 5, kbl )
701 IF( ( lower .OR. upper ) .AND. ( n.GT.max0( 64, 4*kbl ) ) )
THEN
723 CALL sgsvj0( jobv, m-n34, n-n34, a( n34+1, n34+1 ), lda,
724 $ work( n34+1 ), sva( n34+1 ), mvl,
725 $ v( n34*q+1, n34+1 ), ldv, epsln, sfmin, tol,
726 $ 2, work( n+1 ), lwork-n, ierr )
728 CALL sgsvj0( jobv, m-n2, n34-n2, a( n2+1, n2+1 ), lda,
729 $ work( n2+1 ), sva( n2+1 ), mvl,
730 $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 2,
731 $ work( n+1 ), lwork-n, ierr )
733 CALL sgsvj1( jobv, m-n2, n-n2, n4, a( n2+1, n2+1 ), lda,
734 $ work( n2+1 ), sva( n2+1 ), mvl,
735 $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 1,
736 $ work( n+1 ), lwork-n, ierr )
738 CALL sgsvj0( jobv, m-n4, n2-n4, a( n4+1, n4+1 ), lda,
739 $ work( n4+1 ), sva( n4+1 ), mvl,
740 $ v( n4*q+1, n4+1 ), ldv, epsln, sfmin, tol, 1,
741 $ work( n+1 ), lwork-n, ierr )
743 CALL sgsvj0( jobv, m, n4, a, lda, work, sva, mvl, v, ldv,
744 $ epsln, sfmin, tol, 1, work( n+1 ), lwork-n,
747 CALL sgsvj1( jobv, m, n2, n4, a, lda, work, sva, mvl, v,
748 $ ldv, epsln, sfmin, tol, 1, work( n+1 ),
752 ELSE IF( upper )
THEN
755 CALL sgsvj0( jobv, n4, n4, a, lda, work, sva, mvl, v, ldv,
756 $ epsln, sfmin, tol, 2, work( n+1 ), lwork-n,
759 CALL sgsvj0( jobv, n2, n4, a( 1, n4+1 ), lda, work( n4+1 ),
760 $ sva( n4+1 ), mvl, v( n4*q+1, n4+1 ), ldv,
761 $ epsln, sfmin, tol, 1, work( n+1 ), lwork-n,
764 CALL sgsvj1( jobv, n2, n2, n4, a, lda, work, sva, mvl, v,
765 $ ldv, epsln, sfmin, tol, 1, work( n+1 ),
768 CALL sgsvj0( jobv, n2+n4, n4, a( 1, n2+1 ), lda,
769 $ work( n2+1 ), sva( n2+1 ), mvl,
770 $ v( n2*q+1, n2+1 ), ldv, epsln, sfmin, tol, 1,
771 $ work( n+1 ), lwork-n, ierr )
779 DO 1993 i = 1, nsweep
797 igl = ( ibr-1 )*kbl + 1
799 DO 1002 ir1 = 0, min0( lkahead, nbl-ibr )
803 DO 2001 p = igl, min0( igl+kbl-1, n-1 )
807 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
809 CALL sswap( m, a( 1, p ), 1, a( 1, q ), 1 )
810 IF( rsvec )
CALL sswap( mvl, v( 1, p ), 1,
816 work( p ) = work( q )
834 IF( ( sva( p ).LT.rootbig ) .AND.
835 $ ( sva( p ).GT.rootsfmin ) )
THEN
836 sva( p ) = snrm2( m, a( 1, p ), 1 )*work( p )
840 CALL slassq( m, a( 1, p ), 1, temp1, aapp )
841 sva( p ) = temp1*sqrt( aapp )*work( p )
848 IF( aapp.GT.zero )
THEN
852 DO 2002 q = p + 1, min0( igl+kbl-1, n )
856 IF( aaqq.GT.zero )
THEN
859 IF( aaqq.GE.one )
THEN
860 rotok = ( small*aapp ).LE.aaqq
861 IF( aapp.LT.( big / aaqq ) )
THEN
862 aapq = ( sdot( m, a( 1, p ), 1, a( 1,
863 $ q ), 1 )*work( p )*work( q ) /
866 CALL scopy( m, a( 1, p ), 1,
868 CALL slascl(
'G', 0, 0, aapp,
870 $ work( n+1 ), lda, ierr )
871 aapq = sdot( m, work( n+1 ), 1,
872 $ a( 1, q ), 1 )*work( q ) / aaqq
875 rotok = aapp.LE.( aaqq / small )
876 IF( aapp.GT.( small / aaqq ) )
THEN
877 aapq = ( sdot( m, a( 1, p ), 1, a( 1,
878 $ q ), 1 )*work( p )*work( q ) /
881 CALL scopy( m, a( 1, q ), 1,
883 CALL slascl(
'G', 0, 0, aaqq,
885 $ work( n+1 ), lda, ierr )
886 aapq = sdot( m, work( n+1 ), 1,
887 $ a( 1, p ), 1 )*work( p ) / aapp
891 mxaapq = amax1( mxaapq, abs( aapq ) )
895 IF( abs( aapq ).GT.tol )
THEN
910 theta = -half*abs( aqoap-apoaq ) / aapq
912 IF( abs( theta ).GT.bigtheta )
THEN
915 fastr( 3 ) = t*work( p ) / work( q )
916 fastr( 4 ) = -t*work( q ) /
918 CALL srotm( m, a( 1, p ), 1,
919 $ a( 1, q ), 1, fastr )
920 IF( rsvec )
CALL srotm( mvl,
924 sva( q ) = aaqq*sqrt( amax1( zero,
925 $ one+t*apoaq*aapq ) )
926 aapp = aapp*sqrt( amax1( zero,
927 $ one-t*aqoap*aapq ) )
928 mxsinj = amax1( mxsinj, abs( t ) )
934 thsign = -sign( one, aapq )
935 t = one / ( theta+thsign*
936 $ sqrt( one+theta*theta ) )
937 cs = sqrt( one / ( one+t*t ) )
940 mxsinj = amax1( mxsinj, abs( sn ) )
941 sva( q ) = aaqq*sqrt( amax1( zero,
942 $ one+t*apoaq*aapq ) )
943 aapp = aapp*sqrt( amax1( zero,
944 $ one-t*aqoap*aapq ) )
946 apoaq = work( p ) / work( q )
947 aqoap = work( q ) / work( p )
948 IF( work( p ).GE.one )
THEN
949 IF( work( q ).GE.one )
THEN
951 fastr( 4 ) = -t*aqoap
952 work( p ) = work( p )*cs
953 work( q ) = work( q )*cs
954 CALL srotm( m, a( 1, p ), 1,
957 IF( rsvec )
CALL srotm( mvl,
958 $ v( 1, p ), 1, v( 1, q ),
961 CALL saxpy( m, -t*aqoap,
964 CALL saxpy( m, cs*sn*apoaq,
967 work( p ) = work( p )*cs
968 work( q ) = work( q ) / cs
970 CALL saxpy( mvl, -t*aqoap,
980 IF( work( q ).GE.one )
THEN
981 CALL saxpy( m, t*apoaq,
984 CALL saxpy( m, -cs*sn*aqoap,
987 work( p ) = work( p ) / cs
988 work( q ) = work( q )*cs
990 CALL saxpy( mvl, t*apoaq,
999 IF( work( p ).GE.work( q ) )
1001 CALL saxpy( m, -t*aqoap,
1004 CALL saxpy( m, cs*sn*apoaq,
1007 work( p ) = work( p )*cs
1008 work( q ) = work( q ) / cs
1020 CALL saxpy( m, t*apoaq,
1027 work( p ) = work( p ) / cs
1028 work( q ) = work( q )*cs
1031 $ t*apoaq, v( 1, p ),
1045 CALL scopy( m, a( 1, p ), 1,
1047 CALL slascl(
'G', 0, 0, aapp, one, m,
1048 $ 1, work( n+1 ), lda,
1050 CALL slascl(
'G', 0, 0, aaqq, one, m,
1051 $ 1, a( 1, q ), lda, ierr )
1052 temp1 = -aapq*work( p ) / work( q )
1053 CALL saxpy( m, temp1, work( n+1 ), 1,
1055 CALL slascl(
'G', 0, 0, one, aaqq, m,
1056 $ 1, a( 1, q ), lda, ierr )
1057 sva( q ) = aaqq*sqrt( amax1( zero,
1059 mxsinj = amax1( mxsinj, sfmin )
1066 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
1068 IF( ( aaqq.LT.rootbig ) .AND.
1069 $ ( aaqq.GT.rootsfmin ) )
THEN
1070 sva( q ) = snrm2( m, a( 1, q ), 1 )*
1075 CALL slassq( m, a( 1, q ), 1, t,
1077 sva( q ) = t*sqrt( aaqq )*work( q )
1080 IF( ( aapp / aapp0 ).LE.rooteps )
THEN
1081 IF( ( aapp.LT.rootbig ) .AND.
1082 $ ( aapp.GT.rootsfmin ) )
THEN
1083 aapp = snrm2( m, a( 1, p ), 1 )*
1088 CALL slassq( m, a( 1, p ), 1, t,
1090 aapp = t*sqrt( aapp )*work( p )
1097 IF( ir1.EQ.0 )notrot = notrot + 1
1099 pskipped = pskipped + 1
1103 IF( ir1.EQ.0 )notrot = notrot + 1
1104 pskipped = pskipped + 1
1107 IF( ( i.LE.swband ) .AND.
1108 $ ( pskipped.GT.rowskip ) )
THEN
1109 IF( ir1.EQ.0 )aapp = -aapp
1124 IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
1125 $ notrot = notrot + min0( igl+kbl-1, n ) - p
1136 igl = ( ibr-1 )*kbl + 1
1138 DO 2010 jbc = ibr + 1, nbl
1140 jgl = ( jbc-1 )*kbl + 1
1145 DO 2100 p = igl, min0( igl+kbl-1, n )
1148 IF( aapp.GT.zero )
THEN
1152 DO 2200 q = jgl, min0( jgl+kbl-1, n )
1155 IF( aaqq.GT.zero )
THEN
1162 IF( aaqq.GE.one )
THEN
1163 IF( aapp.GE.aaqq )
THEN
1164 rotok = ( small*aapp ).LE.aaqq
1166 rotok = ( small*aaqq ).LE.aapp
1168 IF( aapp.LT.( big / aaqq ) )
THEN
1169 aapq = ( sdot( m, a( 1, p ), 1, a( 1,
1170 $ q ), 1 )*work( p )*work( q ) /
1173 CALL scopy( m, a( 1, p ), 1,
1175 CALL slascl(
'G', 0, 0, aapp,
1177 $ work( n+1 ), lda, ierr )
1178 aapq = sdot( m, work( n+1 ), 1,
1179 $ a( 1, q ), 1 )*work( q ) / aaqq
1182 IF( aapp.GE.aaqq )
THEN
1183 rotok = aapp.LE.( aaqq / small )
1185 rotok = aaqq.LE.( aapp / small )
1187 IF( aapp.GT.( small / aaqq ) )
THEN
1188 aapq = ( sdot( m, a( 1, p ), 1, a( 1,
1189 $ q ), 1 )*work( p )*work( q ) /
1192 CALL scopy( m, a( 1, q ), 1,
1194 CALL slascl(
'G', 0, 0, aaqq,
1196 $ work( n+1 ), lda, ierr )
1197 aapq = sdot( m, work( n+1 ), 1,
1198 $ a( 1, p ), 1 )*work( p ) / aapp
1202 mxaapq = amax1( mxaapq, abs( aapq ) )
1206 IF( abs( aapq ).GT.tol )
THEN
1216 theta = -half*abs( aqoap-apoaq ) / aapq
1217 IF( aaqq.GT.aapp0 )theta = -theta
1219 IF( abs( theta ).GT.bigtheta )
THEN
1221 fastr( 3 ) = t*work( p ) / work( q )
1222 fastr( 4 ) = -t*work( q ) /
1224 CALL srotm( m, a( 1, p ), 1,
1225 $ a( 1, q ), 1, fastr )
1226 IF( rsvec )
CALL srotm( mvl,
1230 sva( q ) = aaqq*sqrt( amax1( zero,
1231 $ one+t*apoaq*aapq ) )
1232 aapp = aapp*sqrt( amax1( zero,
1233 $ one-t*aqoap*aapq ) )
1234 mxsinj = amax1( mxsinj, abs( t ) )
1239 thsign = -sign( one, aapq )
1240 IF( aaqq.GT.aapp0 )thsign = -thsign
1241 t = one / ( theta+thsign*
1242 $ sqrt( one+theta*theta ) )
1243 cs = sqrt( one / ( one+t*t ) )
1245 mxsinj = amax1( mxsinj, abs( sn ) )
1246 sva( q ) = aaqq*sqrt( amax1( zero,
1247 $ one+t*apoaq*aapq ) )
1248 aapp = aapp*sqrt( amax1( zero,
1249 $ one-t*aqoap*aapq ) )
1251 apoaq = work( p ) / work( q )
1252 aqoap = work( q ) / work( p )
1253 IF( work( p ).GE.one )
THEN
1255 IF( work( q ).GE.one )
THEN
1256 fastr( 3 ) = t*apoaq
1257 fastr( 4 ) = -t*aqoap
1258 work( p ) = work( p )*cs
1259 work( q ) = work( q )*cs
1260 CALL srotm( m, a( 1, p ), 1,
1263 IF( rsvec )
CALL srotm( mvl,
1264 $ v( 1, p ), 1, v( 1, q ),
1267 CALL saxpy( m, -t*aqoap,
1270 CALL saxpy( m, cs*sn*apoaq,
1274 CALL saxpy( mvl, -t*aqoap,
1282 work( p ) = work( p )*cs
1283 work( q ) = work( q ) / cs
1286 IF( work( q ).GE.one )
THEN
1287 CALL saxpy( m, t*apoaq,
1290 CALL saxpy( m, -cs*sn*aqoap,
1294 CALL saxpy( mvl, t*apoaq,
1302 work( p ) = work( p ) / cs
1303 work( q ) = work( q )*cs
1305 IF( work( p ).GE.work( q ) )
1307 CALL saxpy( m, -t*aqoap,
1310 CALL saxpy( m, cs*sn*apoaq,
1313 work( p ) = work( p )*cs
1314 work( q ) = work( q ) / cs
1326 CALL saxpy( m, t*apoaq,
1333 work( p ) = work( p ) / cs
1334 work( q ) = work( q )*cs
1337 $ t*apoaq, v( 1, p ),
1350 IF( aapp.GT.aaqq )
THEN
1351 CALL scopy( m, a( 1, p ), 1,
1353 CALL slascl(
'G', 0, 0, aapp, one,
1354 $ m, 1, work( n+1 ), lda,
1356 CALL slascl(
'G', 0, 0, aaqq, one,
1357 $ m, 1, a( 1, q ), lda,
1359 temp1 = -aapq*work( p ) / work( q )
1360 CALL saxpy( m, temp1, work( n+1 ),
1362 CALL slascl(
'G', 0, 0, one, aaqq,
1363 $ m, 1, a( 1, q ), lda,
1365 sva( q ) = aaqq*sqrt( amax1( zero,
1367 mxsinj = amax1( mxsinj, sfmin )
1369 CALL scopy( m, a( 1, q ), 1,
1371 CALL slascl(
'G', 0, 0, aaqq, one,
1372 $ m, 1, work( n+1 ), lda,
1374 CALL slascl(
'G', 0, 0, aapp, one,
1375 $ m, 1, a( 1, p ), lda,
1377 temp1 = -aapq*work( q ) / work( p )
1378 CALL saxpy( m, temp1, work( n+1 ),
1380 CALL slascl(
'G', 0, 0, one, aapp,
1381 $ m, 1, a( 1, p ), lda,
1383 sva( p ) = aapp*sqrt( amax1( zero,
1385 mxsinj = amax1( mxsinj, sfmin )
1392 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
1394 IF( ( aaqq.LT.rootbig ) .AND.
1395 $ ( aaqq.GT.rootsfmin ) )
THEN
1396 sva( q ) = snrm2( m, a( 1, q ), 1 )*
1401 CALL slassq( m, a( 1, q ), 1, t,
1403 sva( q ) = t*sqrt( aaqq )*work( q )
1406 IF( ( aapp / aapp0 )**2.LE.rooteps )
THEN
1407 IF( ( aapp.LT.rootbig ) .AND.
1408 $ ( aapp.GT.rootsfmin ) )
THEN
1409 aapp = snrm2( m, a( 1, p ), 1 )*
1414 CALL slassq( m, a( 1, p ), 1, t,
1416 aapp = t*sqrt( aapp )*work( p )
1424 pskipped = pskipped + 1
1429 pskipped = pskipped + 1
1433 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
1439 IF( ( i.LE.swband ) .AND.
1440 $ ( pskipped.GT.rowskip ) )
THEN
1454 IF( aapp.EQ.zero )notrot = notrot +
1455 $ min0( jgl+kbl-1, n ) - jgl + 1
1456 IF( aapp.LT.zero )notrot = 0
1466 DO 2012 p = igl, min0( igl+kbl-1, n )
1467 sva( p ) = abs( sva( p ) )
1474 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
1476 sva( n ) = snrm2( m, a( 1, n ), 1 )*work( n )
1480 CALL slassq( m, a( 1, n ), 1, t, aapp )
1481 sva( n ) = t*sqrt( aapp )*work( n )
1486 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
1487 $ ( iswrot.LE.n ) ) )swband = i
1489 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt( float( n ) )*
1490 $ tol ) .AND. ( float( n )*mxaapq*mxsinj.LT.tol ) )
THEN
1494 IF( notrot.GE.emptsw )
GO TO 1994
1516 DO 5991 p = 1, n - 1
1517 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
1523 work( p ) = work( q )
1525 CALL sswap( m, a( 1, p ), 1, a( 1, q ), 1 )
1526 IF( rsvec )
CALL sswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
1528 IF( sva( p ).NE.zero )
THEN
1530 IF( sva( p )*skl.GT.sfmin )n2 = n2 + 1
1533 IF( sva( n ).NE.zero )
THEN
1535 IF( sva( n )*skl.GT.sfmin )n2 = n2 + 1
1540 IF( lsvec .OR. uctol )
THEN
1542 CALL sscal( m, work( p ) / sva( p ), a( 1, p ), 1 )
1551 CALL sscal( mvl, work( p ), v( 1, p ), 1 )
1555 temp1 = one / snrm2( mvl, v( 1, p ), 1 )
1556 CALL sscal( mvl, temp1, v( 1, p ), 1 )
1562 IF( ( ( skl.GT.one ) .AND. ( sva( 1 ).LT.( big / skl ) ) )
1563 $ .OR. ( ( skl.LT.one ) .AND. ( sva( max( n2, 1 ) ) .GT.
1564 $ ( sfmin / skl ) ) ) )
THEN
1566 sva( p ) = skl*sva( p )
1576 work( 2 ) = float( n4 )
1579 work( 3 ) = float( n2 )
1584 work( 4 ) = float( i )
subroutine slassq(N, X, INCX, SCALE, SUMSQ)
SLASSQ updates a sum of squares represented in scaled form.
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
subroutine sgsvj0(JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
SGSVJ0 pre-processor for the routine sgesvj.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine srotm(N, SX, INCX, SY, INCY, SPARAM)
SROTM
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine sgsvj1(JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
SGSVJ1 pre-processor for the routine sgesvj, applies Jacobi rotations targeting only particular pivot...
subroutine sgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
SGESVJ
subroutine sscal(N, SA, SX, INCX)
SSCAL