473 SUBROUTINE sgejsv( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
474 $ m, n, a, lda, sva, u, ldu, v, ldv,
475 $ work, lwork, iwork, info )
484 INTEGER INFO, LDA, LDU, LDV, LWORK, M, N
487 REAL A( lda, * ), SVA( n ), U( ldu, * ), V( ldv, * ),
490 CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
497 parameter( zero = 0.0e0, one = 1.0e0 )
500 REAL AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK,
501 $ condr1, condr2, entra, entrat, epsln, maxprj, scalem,
502 $ sconda, sfmin, small, temp1, uscal1, uscal2, xsc
503 INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING
504 LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LSVEC,
505 $ l2aber, l2kill, l2pert, l2rank, l2tran,
506 $ noscal, rowpiv, rsvec, transp
509 INTRINSIC abs, alog, amax1, amin1, float,
510 $ max0, min0, nint, sign, sqrt
516 EXTERNAL isamax, lsame, slamch, snrm2
528 lsvec = lsame( jobu,
'U' ) .OR. lsame( jobu,
'F' )
529 jracc = lsame( jobv,
'J' )
530 rsvec = lsame( jobv,
'V' ) .OR. jracc
531 rowpiv = lsame( joba,
'F' ) .OR. lsame( joba,
'G' )
532 l2rank = lsame( joba,
'R' )
533 l2aber = lsame( joba,
'A' )
534 errest = lsame( joba,
'E' ) .OR. lsame( joba,
'G' )
535 l2tran = lsame( jobt,
'T' )
536 l2kill = lsame( jobr,
'R' )
537 defr = lsame( jobr,
'N' )
538 l2pert = lsame( jobp,
'P' )
540 IF ( .NOT.(rowpiv .OR. l2rank .OR. l2aber .OR.
541 $ errest .OR. lsame( joba,
'C' ) ))
THEN
543 ELSE IF ( .NOT.( lsvec .OR. lsame( jobu,
'N' ) .OR.
544 $ lsame( jobu,
'W' )) )
THEN
546 ELSE IF ( .NOT.( rsvec .OR. lsame( jobv,
'N' ) .OR.
547 $ lsame( jobv,
'W' )) .OR. ( jracc .AND. (.NOT.lsvec) ) )
THEN
549 ELSE IF ( .NOT. ( l2kill .OR. defr ) )
THEN
551 ELSE IF ( .NOT. ( l2tran .OR. lsame( jobt,
'N' ) ) )
THEN
553 ELSE IF ( .NOT. ( l2pert .OR. lsame( jobp,
'N' ) ) )
THEN
555 ELSE IF ( m .LT. 0 )
THEN
557 ELSE IF ( ( n .LT. 0 ) .OR. ( n .GT. m ) )
THEN
559 ELSE IF ( lda .LT. m )
THEN
561 ELSE IF ( lsvec .AND. ( ldu .LT. m ) )
THEN
563 ELSE IF ( rsvec .AND. ( ldv .LT. n ) )
THEN
565 ELSE IF ( (.NOT.(lsvec .OR. rsvec .OR. errest).AND.
566 $ (lwork .LT. max0(7,4*n+1,2*m+n))) .OR.
567 $ (.NOT.(lsvec .OR. rsvec) .AND. errest .AND.
568 $ (lwork .LT. max0(7,4*n+n*n,2*m+n))) .OR.
569 $ (lsvec .AND. (.NOT.rsvec) .AND. (lwork .LT. max0(7,2*m+n,4*n+1)))
571 $ (rsvec .AND. (.NOT.lsvec) .AND. (lwork .LT. max0(7,2*m+n,4*n+1)))
573 $ (lsvec .AND. rsvec .AND. (.NOT.jracc) .AND.
574 $ (lwork.LT.max0(2*m+n,6*n+2*n*n)))
575 $ .OR. (lsvec .AND. rsvec .AND. jracc .AND.
576 $ lwork.LT.max0(2*m+n,4*n+n*n,2*n+n*n+6)))
584 IF ( info .NE. 0 )
THEN
586 CALL
xerbla(
'SGEJSV', - info )
592 IF ( ( m .EQ. 0 ) .OR. ( n .EQ. 0 ) )
RETURN
598 IF ( lsame( jobu,
'F' ) ) n1 = m
605 epsln = slamch(
'Epsilon')
606 sfmin = slamch(
'SafeMinimum')
607 small = sfmin / epsln
617 scalem = one / sqrt(float(m)*float(n))
623 CALL
slassq( m, a(1,p), 1, aapp, aaqq )
624 IF ( aapp .GT. big )
THEN
626 CALL
xerbla(
'SGEJSV', -info )
630 IF ( ( aapp .LT. (big / aaqq) ) .AND. noscal )
THEN
634 sva(p) = aapp * ( aaqq * scalem )
637 CALL
sscal( p-1, scalem, sva, 1 )
642 IF ( noscal ) scalem = one
647 aapp = amax1( aapp, sva(p) )
648 IF ( sva(p) .NE. zero ) aaqq = amin1( aaqq, sva(p) )
653 IF ( aapp .EQ. zero )
THEN
654 IF ( lsvec ) CALL
slaset(
'G', m, n1, zero, one, u, ldu )
655 IF ( rsvec ) CALL
slaset(
'G', n, n, zero, one, v, ldv )
658 IF ( errest ) work(3) = one
659 IF ( lsvec .AND. rsvec )
THEN
678 IF ( aaqq .LE. sfmin )
THEN
689 CALL
slascl(
'G',0,0,sva(1),scalem, m,1,a(1,1),lda,ierr )
690 CALL
slacpy(
'A', m, 1, a, lda, u, ldu )
692 IF ( n1 .NE. n )
THEN
693 CALL
sgeqrf( m, n, u,ldu, work, work(n+1),lwork-n,ierr )
694 CALL
sorgqr( m,n1,1, u,ldu,work,work(n+1),lwork-n,ierr )
695 CALL
scopy( m, a(1,1), 1, u(1,1), 1 )
701 IF ( sva(1) .LT. (big*scalem) )
THEN
702 sva(1) = sva(1) / scalem
705 work(1) = one / scalem
707 IF ( sva(1) .NE. zero )
THEN
709 IF ( ( sva(1) / scalem) .GE. sfmin )
THEN
718 IF ( errest ) work(3) = one
719 IF ( lsvec .AND. rsvec )
THEN
732 l2tran = l2tran .AND. ( m .EQ. n )
736 IF ( rowpiv .OR. l2tran )
THEN
747 CALL
slassq( n, a(p,1), lda, xsc, temp1 )
750 work(m+n+p) = xsc * scalem
751 work(n+p) = xsc * (scalem*sqrt(temp1))
752 aatmax = amax1( aatmax, work(n+p) )
753 IF (work(n+p) .NE. zero) aatmin = amin1(aatmin,work(n+p))
757 work(m+n+p) = scalem*abs( a(p,isamax(n,a(p,1),lda)) )
758 aatmax = amax1( aatmax, work(m+n+p) )
759 aatmin = amin1( aatmin, work(m+n+p) )
778 CALL
slassq( n, sva, 1, xsc, temp1 )
783 big1 = ( ( sva(p) / xsc )**2 ) * temp1
784 IF ( big1 .NE. zero ) entra = entra + big1 * alog(big1)
786 entra = - entra / alog(float(n))
796 big1 = ( ( work(p) / xsc )**2 ) * temp1
797 IF ( big1 .NE. zero ) entrat = entrat + big1 * alog(big1)
799 entrat = - entrat / alog(float(m))
804 transp = ( entrat .LT. entra )
850 temp1 = sqrt( big / float(n) )
852 CALL
slascl(
'G', 0, 0, aapp, temp1, n, 1, sva, n, ierr )
853 IF ( aaqq .GT. (aapp * sfmin) )
THEN
854 aaqq = ( aaqq / aapp ) * temp1
856 aaqq = ( aaqq * temp1 ) / aapp
858 temp1 = temp1 * scalem
859 CALL
slascl(
'G', 0, 0, aapp, temp1, m, n, a, lda, ierr )
883 IF ( ( aaqq.LT.sqrt(sfmin) ) .AND. lsvec .AND. rsvec )
THEN
888 IF ( aaqq .LT. xsc )
THEN
890 IF ( sva(p) .LT. xsc )
THEN
891 CALL
slaset(
'A', m, 1, zero, zero, a(1,p), lda )
906 q = isamax( m-p+1, work(m+n+p), 1 ) + p - 1
910 work(m+n+p) = work(m+n+q)
914 CALL
slaswp( n, a, lda, 1, m-1, iwork(2*n+1), 1 )
936 CALL
sgeqp3( m,n,a,lda, iwork,work, work(n+1),lwork-n, ierr )
952 temp1 = sqrt(float(n))*epsln
954 IF ( abs(a(p,p)) .GE. (temp1*abs(a(1,1))) )
THEN
961 ELSE IF ( l2rank )
THEN
967 IF ( ( abs(a(p,p)) .LT. (epsln*abs(a(p-1,p-1))) ) .OR.
968 $ ( abs(a(p,p)) .LT. small ) .OR.
969 $ ( l2kill .AND. (abs(a(p,p)) .LT. temp1) ) ) go to 3402
984 IF ( ( abs(a(p,p)) .LT. small ) .OR.
985 $ ( l2kill .AND. (abs(a(p,p)) .LT. temp1) ) ) go to 3302
993 IF ( nr .EQ. n )
THEN
996 temp1 = abs(a(p,p)) / sva(iwork(p))
997 maxprj = amin1( maxprj, temp1 )
999 IF ( maxprj**2 .GE. one - float(n)*epsln ) almort = .true.
1008 IF ( n .EQ. nr )
THEN
1011 CALL
slacpy(
'U', n, n, a, lda, v, ldv )
1013 temp1 = sva(iwork(p))
1014 CALL
sscal( p, one/temp1, v(1,p), 1 )
1016 CALL
spocon(
'U', n, v, ldv, one, temp1,
1017 $ work(n+1), iwork(2*n+m+1), ierr )
1018 ELSE IF ( lsvec )
THEN
1020 CALL
slacpy(
'U', n, n, a, lda, u, ldu )
1022 temp1 = sva(iwork(p))
1023 CALL
sscal( p, one/temp1, u(1,p), 1 )
1025 CALL
spocon(
'U', n, u, ldu, one, temp1,
1026 $ work(n+1), iwork(2*n+m+1), ierr )
1028 CALL
slacpy(
'U', n, n, a, lda, work(n+1), n )
1030 temp1 = sva(iwork(p))
1031 CALL
sscal( p, one/temp1, work(n+(p-1)*n+1), 1 )
1034 CALL
spocon(
'U', n, work(n+1), n, one, temp1,
1035 $ work(n+n*n+1), iwork(2*n+m+1), ierr )
1037 sconda = one / sqrt(temp1)
1045 l2pert = l2pert .AND. ( abs( a(1,1)/a(nr,nr) ) .GT. sqrt(big1) )
1050 IF ( .NOT. ( rsvec .OR. lsvec ) )
THEN
1055 DO 1946 p = 1, min0( n-1, nr )
1056 CALL
scopy( n-p, a(p,p+1), lda, a(p+1,p), 1 )
1071 IF ( .NOT. almort )
THEN
1075 xsc = epsln / float(n)
1077 temp1 = xsc*abs(a(q,q))
1079 IF ( ( (p.GT.q) .AND. (abs(a(p,q)).LE.temp1) )
1080 $ .OR. ( p .LT. q ) )
1081 $ a(p,q) = sign( temp1, a(p,q) )
1085 CALL
slaset(
'U', nr-1,nr-1, zero,zero, a(1,2),lda )
1090 CALL
sgeqrf( n,nr, a,lda, work, work(n+1),lwork-n, ierr )
1093 DO 1948 p = 1, nr - 1
1094 CALL
scopy( nr-p, a(p,p+1), lda, a(p+1,p), 1 )
1105 xsc = epsln / float(n)
1107 temp1 = xsc*abs(a(q,q))
1109 IF ( ( (p.GT.q) .AND. (abs(a(p,q)).LE.temp1) )
1110 $ .OR. ( p .LT. q ) )
1111 $ a(p,q) = sign( temp1, a(p,q) )
1115 CALL
slaset(
'U', nr-1, nr-1, zero, zero, a(1,2), lda )
1122 CALL
sgesvj(
'L',
'NoU',
'NoV', nr, nr, a, lda, sva,
1123 $ n, v, ldv, work, lwork, info )
1126 numrank = nint(work(2))
1129 ELSE IF ( rsvec .AND. ( .NOT. lsvec ) )
THEN
1137 CALL
scopy( n-p+1, a(p,p), lda, v(p,p), 1 )
1139 CALL
slaset(
'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
1141 CALL
sgesvj(
'L',
'U',
'N', n, nr, v,ldv, sva, nr, a,lda,
1142 $ work, lwork, info )
1144 numrank = nint(work(2))
1151 CALL
slaset(
'Lower', nr-1, nr-1, zero, zero, a(2,1), lda )
1152 CALL
sgelqf( nr, n, a, lda, work, work(n+1), lwork-n, ierr)
1153 CALL
slacpy(
'Lower', nr, nr, a, lda, v, ldv )
1154 CALL
slaset(
'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
1155 CALL
sgeqrf( nr, nr, v, ldv, work(n+1), work(2*n+1),
1158 CALL
scopy( nr-p+1, v(p,p), ldv, v(p,p), 1 )
1160 CALL
slaset(
'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
1162 CALL
sgesvj(
'Lower',
'U',
'N', nr, nr, v,ldv, sva, nr, u,
1163 $ ldu, work(n+1), lwork-n, info )
1165 numrank = nint(work(n+2))
1166 IF ( nr .LT. n )
THEN
1167 CALL
slaset(
'A',n-nr, nr, zero,zero, v(nr+1,1), ldv )
1168 CALL
slaset(
'A',nr, n-nr, zero,zero, v(1,nr+1), ldv )
1169 CALL
slaset(
'A',n-nr,n-nr,zero,one, v(nr+1,nr+1), ldv )
1172 CALL
sormlq(
'Left',
'Transpose', n, n, nr, a, lda, work,
1173 $ v, ldv, work(n+1), lwork-n, ierr )
1178 CALL
scopy( n, v(p,1), ldv, a(iwork(p),1), lda )
1180 CALL
slacpy(
'All', n, n, a, lda, v, ldv )
1183 CALL
slacpy(
'All', n, n, v, ldv, u, ldu )
1186 ELSE IF ( lsvec .AND. ( .NOT. rsvec ) )
THEN
1193 CALL
scopy( n-p+1, a(p,p), lda, u(p,p), 1 )
1195 CALL
slaset(
'Upper', nr-1, nr-1, zero, zero, u(1,2), ldu )
1197 CALL
sgeqrf( n, nr, u, ldu, work(n+1), work(2*n+1),
1200 DO 1967 p = 1, nr - 1
1201 CALL
scopy( nr-p, u(p,p+1), ldu, u(p+1,p), 1 )
1203 CALL
slaset(
'Upper', nr-1, nr-1, zero, zero, u(1,2), ldu )
1205 CALL
sgesvj(
'Lower',
'U',
'N', nr,nr, u, ldu, sva, nr, a,
1206 $ lda, work(n+1), lwork-n, info )
1208 numrank = nint(work(n+2))
1210 IF ( nr .LT. m )
THEN
1211 CALL
slaset(
'A', m-nr, nr,zero, zero, u(nr+1,1), ldu )
1212 IF ( nr .LT. n1 )
THEN
1213 CALL
slaset(
'A',nr, n1-nr, zero, zero, u(1,nr+1), ldu )
1214 CALL
slaset(
'A',m-nr,n1-nr,zero,one,u(nr+1,nr+1), ldu )
1218 CALL
sormqr(
'Left',
'No Tr', m, n1, n, a, lda, work, u,
1219 $ ldu, work(n+1), lwork-n, ierr )
1222 $ CALL
slaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1225 xsc = one / snrm2( m, u(1,p), 1 )
1226 CALL
sscal( m, xsc, u(1,p), 1 )
1230 CALL
slacpy(
'All', n, n, u, ldu, v, ldv )
1237 IF ( .NOT. jracc )
THEN
1239 IF ( .NOT. almort )
THEN
1249 CALL
scopy( n-p+1, a(p,p), lda, v(p,p), 1 )
1267 temp1 = xsc*abs( v(q,q) )
1269 IF ( ( p .GT. q ) .AND. ( abs(v(p,q)) .LE. temp1 )
1270 $ .OR. ( p .LT. q ) )
1271 $ v(p,q) = sign( temp1, v(p,q) )
1272 IF ( p .LT. q ) v(p,q) = - v(p,q)
1276 CALL
slaset(
'U', nr-1, nr-1, zero, zero, v(1,2), ldv )
1283 CALL
slacpy(
'L', nr, nr, v, ldv, work(2*n+1), nr )
1285 temp1 = snrm2(nr-p+1,work(2*n+(p-1)*nr+p),1)
1286 CALL
sscal(nr-p+1,one/temp1,work(2*n+(p-1)*nr+p),1)
1288 CALL
spocon(
'Lower',nr,work(2*n+1),nr,one,temp1,
1289 $ work(2*n+nr*nr+1),iwork(m+2*n+1),ierr)
1290 condr1 = one / sqrt(temp1)
1296 cond_ok = sqrt(float(nr))
1299 IF ( condr1 .LT. cond_ok )
THEN
1304 CALL
sgeqrf( n, nr, v, ldv, work(n+1), work(2*n+1),
1308 xsc = sqrt(small)/epsln
1310 DO 3958 q = 1, p - 1
1311 temp1 = xsc * amin1(abs(v(p,p)),abs(v(q,q)))
1312 IF ( abs(v(q,p)) .LE. temp1 )
1313 $ v(q,p) = sign( temp1, v(q,p) )
1319 $ CALL
slacpy(
'A', n, nr, v, ldv, work(2*n+1), n )
1323 DO 1969 p = 1, nr - 1
1324 CALL
scopy( nr-p, v(p,p+1), ldv, v(p+1,p), 1 )
1342 CALL
sgeqp3( n, nr, v, ldv, iwork(n+1), work(n+1),
1343 $ work(2*n+1), lwork-2*n, ierr )
1349 DO 3968 q = 1, p - 1
1350 temp1 = xsc * amin1(abs(v(p,p)),abs(v(q,q)))
1351 IF ( abs(v(q,p)) .LE. temp1 )
1352 $ v(q,p) = sign( temp1, v(q,p) )
1357 CALL
slacpy(
'A', n, nr, v, ldv, work(2*n+1), n )
1362 DO 8971 q = 1, p - 1
1363 temp1 = xsc * amin1(abs(v(p,p)),abs(v(q,q)))
1364 v(p,q) = - sign( temp1, v(q,p) )
1368 CALL
slaset(
'L',nr-1,nr-1,zero,zero,v(2,1),ldv )
1371 CALL
sgelqf( nr, nr, v, ldv, work(2*n+n*nr+1),
1372 $ work(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr, ierr )
1374 CALL
slacpy(
'L',nr,nr,v,ldv,work(2*n+n*nr+nr+1),nr )
1376 temp1 = snrm2( p, work(2*n+n*nr+nr+p), nr )
1377 CALL
sscal( p, one/temp1, work(2*n+n*nr+nr+p), nr )
1379 CALL
spocon(
'L',nr,work(2*n+n*nr+nr+1),nr,one,temp1,
1380 $ work(2*n+n*nr+nr+nr*nr+1),iwork(m+2*n+1),ierr )
1381 condr2 = one / sqrt(temp1)
1383 IF ( condr2 .GE. cond_ok )
THEN
1388 CALL
slacpy(
'U', nr, nr, v, ldv, work(2*n+1), n )
1398 temp1 = xsc * v(q,q)
1399 DO 4969 p = 1, q - 1
1401 v(p,q) = - sign( temp1, v(p,q) )
1405 CALL
slaset(
'U', nr-1,nr-1, zero,zero, v(1,2), ldv )
1414 IF ( condr1 .LT. cond_ok )
THEN
1416 CALL
sgesvj(
'L',
'U',
'N',nr,nr,v,ldv,sva,nr,u,
1417 $ ldu,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,info )
1418 scalem = work(2*n+n*nr+nr+1)
1419 numrank = nint(work(2*n+n*nr+nr+2))
1421 CALL
scopy( nr, v(1,p), 1, u(1,p), 1 )
1422 CALL
sscal( nr, sva(p), v(1,p), 1 )
1427 IF ( nr .EQ. n )
THEN
1432 CALL
strsm(
'L',
'U',
'N',
'N', nr,nr,one, a,lda, v,ldv )
1438 CALL
strsm(
'L',
'U',
'T',
'N',nr,nr,one,work(2*n+1),
1440 IF ( nr .LT. n )
THEN
1441 CALL
slaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv)
1442 CALL
slaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv)
1443 CALL
slaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv)
1445 CALL
sormqr(
'L',
'N',n,n,nr,work(2*n+1),n,work(n+1),
1446 $ v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr)
1449 ELSE IF ( condr2 .LT. cond_ok )
THEN
1457 CALL
sgesvj(
'L',
'U',
'N', nr, nr, v, ldv, sva, nr, u,
1458 $ ldu, work(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr, info )
1459 scalem = work(2*n+n*nr+nr+1)
1460 numrank = nint(work(2*n+n*nr+nr+2))
1462 CALL
scopy( nr, v(1,p), 1, u(1,p), 1 )
1463 CALL
sscal( nr, sva(p), u(1,p), 1 )
1465 CALL
strsm(
'L',
'U',
'N',
'N',nr,nr,one,work(2*n+1),n,u,ldu)
1469 work(2*n+n*nr+nr+iwork(n+p)) = u(p,q)
1472 u(p,q) = work(2*n+n*nr+nr+p)
1475 IF ( nr .LT. n )
THEN
1476 CALL
slaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
1477 CALL
slaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
1478 CALL
slaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
1480 CALL
sormqr(
'L',
'N',n,n,nr,work(2*n+1),n,work(n+1),
1481 $ v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
1494 CALL
sgesvj(
'L',
'U',
'V', nr, nr, v, ldv, sva, nr, u,
1495 $ ldu, work(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr, info )
1496 scalem = work(2*n+n*nr+nr+1)
1497 numrank = nint(work(2*n+n*nr+nr+2))
1498 IF ( nr .LT. n )
THEN
1499 CALL
slaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
1500 CALL
slaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
1501 CALL
slaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
1503 CALL
sormqr(
'L',
'N',n,n,nr,work(2*n+1),n,work(n+1),
1504 $ v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
1506 CALL
sormlq(
'L',
'T', nr, nr, nr, work(2*n+1), n,
1507 $ work(2*n+n*nr+1), u, ldu, work(2*n+n*nr+nr+1),
1508 $ lwork-2*n-n*nr-nr, ierr )
1511 work(2*n+n*nr+nr+iwork(n+p)) = u(p,q)
1514 u(p,q) = work(2*n+n*nr+nr+p)
1524 temp1 = sqrt(float(n)) * epsln
1527 work(2*n+n*nr+nr+iwork(p)) = v(p,q)
1530 v(p,q) = work(2*n+n*nr+nr+p)
1532 xsc = one / snrm2( n, v(1,q), 1 )
1533 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1534 $ CALL
sscal( n, xsc, v(1,q), 1 )
1538 IF ( nr .LT. m )
THEN
1539 CALL
slaset(
'A', m-nr, nr, zero, zero, u(nr+1,1), ldu )
1540 IF ( nr .LT. n1 )
THEN
1541 CALL
slaset(
'A',nr,n1-nr,zero,zero,u(1,nr+1),ldu)
1542 CALL
slaset(
'A',m-nr,n1-nr,zero,one,u(nr+1,nr+1),ldu)
1549 CALL
sormqr(
'Left',
'No_Tr', m, n1, n, a, lda, work, u,
1550 $ ldu, work(n+1), lwork-n, ierr )
1553 temp1 = sqrt(float(m)) * epsln
1555 xsc = one / snrm2( m, u(1,p), 1 )
1556 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1557 $ CALL
sscal( m, xsc, u(1,p), 1 )
1564 $ CALL
slaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1571 CALL
slacpy(
'Upper', n, n, a, lda, work(n+1), n )
1575 temp1 = xsc * work( n + (p-1)*n + p )
1576 DO 5971 q = 1, p - 1
1577 work(n+(q-1)*n+p)=-sign(temp1,work(n+(p-1)*n+q))
1581 CALL
slaset(
'Lower',n-1,n-1,zero,zero,work(n+2),n )
1584 CALL
sgesvj(
'Upper',
'U',
'N', n, n, work(n+1), n, sva,
1585 $ n, u, ldu, work(n+n*n+1), lwork-n-n*n, info )
1587 scalem = work(n+n*n+1)
1588 numrank = nint(work(n+n*n+2))
1590 CALL
scopy( n, work(n+(p-1)*n+1), 1, u(1,p), 1 )
1591 CALL
sscal( n, sva(p), work(n+(p-1)*n+1), 1 )
1594 CALL
strsm(
'Left',
'Upper',
'NoTrans',
'No UD', n, n,
1595 $ one, a, lda, work(n+1), n )
1597 CALL
scopy( n, work(n+p), n, v(iwork(p),1), ldv )
1599 temp1 = sqrt(float(n))*epsln
1601 xsc = one / snrm2( n, v(1,p), 1 )
1602 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1603 $ CALL
sscal( n, xsc, v(1,p), 1 )
1608 IF ( n .LT. m )
THEN
1609 CALL
slaset(
'A', m-n, n, zero, zero, u(n+1,1), ldu )
1610 IF ( n .LT. n1 )
THEN
1611 CALL
slaset(
'A',n, n1-n, zero, zero, u(1,n+1),ldu )
1612 CALL
slaset(
'A',m-n,n1-n, zero, one,u(n+1,n+1),ldu )
1615 CALL
sormqr(
'Left',
'No Tr', m, n1, n, a, lda, work, u,
1616 $ ldu, work(n+1), lwork-n, ierr )
1617 temp1 = sqrt(float(m))*epsln
1619 xsc = one / snrm2( m, u(1,p), 1 )
1620 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1621 $ CALL
sscal( m, xsc, u(1,p), 1 )
1625 $ CALL
slaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1644 CALL
scopy( n-p+1, a(p,p), lda, v(p,p), 1 )
1648 xsc = sqrt(small/epsln)
1650 temp1 = xsc*abs( v(q,q) )
1652 IF ( ( p .GT. q ) .AND. ( abs(v(p,q)) .LE. temp1 )
1653 $ .OR. ( p .LT. q ) )
1654 $ v(p,q) = sign( temp1, v(p,q) )
1655 IF ( p .LT. q ) v(p,q) = - v(p,q)
1659 CALL
slaset(
'U', nr-1, nr-1, zero, zero, v(1,2), ldv )
1662 CALL
sgeqrf( n, nr, v, ldv, work(n+1), work(2*n+1),
1664 CALL
slacpy(
'L', n, nr, v, ldv, work(2*n+1), n )
1667 CALL
scopy( nr-p+1, v(p,p), ldv, u(p,p), 1 )
1671 xsc = sqrt(small/epsln)
1673 DO 9971 p = 1, q - 1
1674 temp1 = xsc * amin1(abs(u(p,p)),abs(u(q,q)))
1675 u(p,q) = - sign( temp1, u(q,p) )
1679 CALL
slaset(
'U', nr-1, nr-1, zero, zero, u(1,2), ldu )
1682 CALL
sgesvj(
'L',
'U',
'V', nr, nr, u, ldu, sva,
1683 $ n, v, ldv, work(2*n+n*nr+1), lwork-2*n-n*nr, info )
1684 scalem = work(2*n+n*nr+1)
1685 numrank = nint(work(2*n+n*nr+2))
1687 IF ( nr .LT. n )
THEN
1688 CALL
slaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
1689 CALL
slaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
1690 CALL
slaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
1693 CALL
sormqr(
'L',
'N',n,n,nr,work(2*n+1),n,work(n+1),
1694 $ v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
1700 temp1 = sqrt(float(n)) * epsln
1703 work(2*n+n*nr+nr+iwork(p)) = v(p,q)
1706 v(p,q) = work(2*n+n*nr+nr+p)
1708 xsc = one / snrm2( n, v(1,q), 1 )
1709 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1710 $ CALL
sscal( n, xsc, v(1,q), 1 )
1716 IF ( nr .LT. m )
THEN
1717 CALL
slaset(
'A', m-nr, nr, zero, zero, u(nr+1,1), ldu )
1718 IF ( nr .LT. n1 )
THEN
1719 CALL
slaset(
'A',nr, n1-nr, zero, zero, u(1,nr+1),ldu )
1720 CALL
slaset(
'A',m-nr,n1-nr, zero, one,u(nr+1,nr+1),ldu )
1724 CALL
sormqr(
'Left',
'No Tr', m, n1, n, a, lda, work, u,
1725 $ ldu, work(n+1), lwork-n, ierr )
1728 $ CALL
slaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1735 CALL
sswap( n, u(1,p), 1, v(1,p), 1 )
1744 IF ( uscal2 .LE. (big/sva(1))*uscal1 )
THEN
1745 CALL
slascl(
'G', 0, 0, uscal1, uscal2, nr, 1, sva, n, ierr )
1750 IF ( nr .LT. n )
THEN
1756 work(1) = uscal2 * scalem
1758 IF ( errest ) work(3) = sconda
1759 IF ( lsvec .AND. rsvec )
THEN
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 sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
subroutine sgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
SGEQP3
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sgejsv(JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, M, N, A, LDA, SVA, U, LDU, V, LDV, WORK, LWORK, IWORK, INFO)
SGEJSV
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine spocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
SPOCON
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine sgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
SGESVJ
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
subroutine slaswp(N, A, LDA, K1, K2, IPIV, INCX)
SLASWP performs a series of row interchanges on a general rectangular matrix.
subroutine sscal(N, SA, SX, INCX)
SSCAL