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
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.
logical function lsame(CA, CB)
LSAME
subroutine dscal(N, DA, DX, INCX)
DSCAL
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
double precision function dlamch(CMACH)
DLAMCH
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
integer function idamax(N, DX, INCX)
IDAMAX
double precision function dnrm2(N, X, INCX)
DNRM2
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...