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