475 SUBROUTINE sgejsv( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP,
476 $ m, n, a, lda, sva, u, ldu, v, ldv,
477 $ work, lwork, iwork, info )
486 INTEGER INFO, LDA, LDU, LDV, LWORK, M, N
489 REAL A( lda, * ), SVA( n ), U( ldu, * ), V( ldv, * ),
492 CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
499 parameter( zero = 0.0e0, one = 1.0e0 )
502 REAL AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK,
503 $ condr1, condr2, entra, entrat, epsln, maxprj, scalem,
504 $ sconda, sfmin, small, temp1, uscal1, uscal2, xsc
505 INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING
506 LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LSVEC,
507 $ l2aber, l2kill, l2pert, l2rank, l2tran,
508 $ noscal, rowpiv, rsvec, transp
511 INTRINSIC abs, alog, max, min, float, nint, sign, sqrt
517 EXTERNAL isamax, lsame, slamch, snrm2
529 lsvec = lsame( jobu,
'U' ) .OR. lsame( jobu,
'F' )
530 jracc = lsame( jobv,
'J' )
531 rsvec = lsame( jobv,
'V' ) .OR. jracc
532 rowpiv = lsame( joba,
'F' ) .OR. lsame( joba,
'G' )
533 l2rank = lsame( joba,
'R' )
534 l2aber = lsame( joba,
'A' )
535 errest = lsame( joba,
'E' ) .OR. lsame( joba,
'G' )
536 l2tran = lsame( jobt,
'T' )
537 l2kill = lsame( jobr,
'R' )
538 defr = lsame( jobr,
'N' )
539 l2pert = lsame( jobp,
'P' )
541 IF ( .NOT.(rowpiv .OR. l2rank .OR. l2aber .OR.
542 $ errest .OR. lsame( joba,
'C' ) ))
THEN
544 ELSE IF ( .NOT.( lsvec .OR. lsame( jobu,
'N' ) .OR.
545 $ lsame( jobu,
'W' )) )
THEN
547 ELSE IF ( .NOT.( rsvec .OR. lsame( jobv,
'N' ) .OR.
548 $ lsame( jobv,
'W' )) .OR. ( jracc .AND. (.NOT.lsvec) ) )
THEN
550 ELSE IF ( .NOT. ( l2kill .OR. defr ) )
THEN
552 ELSE IF ( .NOT. ( l2tran .OR. lsame( jobt,
'N' ) ) )
THEN
554 ELSE IF ( .NOT. ( l2pert .OR. lsame( jobp,
'N' ) ) )
THEN
556 ELSE IF ( m .LT. 0 )
THEN
558 ELSE IF ( ( n .LT. 0 ) .OR. ( n .GT. m ) )
THEN
560 ELSE IF ( lda .LT. m )
THEN
562 ELSE IF ( lsvec .AND. ( ldu .LT. m ) )
THEN
564 ELSE IF ( rsvec .AND. ( ldv .LT. n ) )
THEN
566 ELSE IF ( (.NOT.(lsvec .OR. rsvec .OR. errest).AND.
567 $ (lwork .LT. max(7,4*n+1,2*m+n))) .OR.
568 $ (.NOT.(lsvec .OR. rsvec) .AND. errest .AND.
569 $ (lwork .LT. max(7,4*n+n*n,2*m+n))) .OR.
570 $ (lsvec .AND. (.NOT.rsvec) .AND. (lwork .LT. max(7,2*m+n,4*n+1)))
572 $ (rsvec .AND. (.NOT.lsvec) .AND. (lwork .LT. max(7,2*m+n,4*n+1)))
574 $ (lsvec .AND. rsvec .AND. (.NOT.jracc) .AND.
575 $ (lwork.LT.max(2*m+n,6*n+2*n*n)))
576 $ .OR. (lsvec .AND. rsvec .AND. jracc .AND.
577 $ lwork.LT.max(2*m+n,4*n+n*n,2*n+n*n+6)))
585 IF ( info .NE. 0 )
THEN
587 CALL xerbla(
'SGEJSV', - info )
593 IF ( ( m .EQ. 0 ) .OR. ( n .EQ. 0 ) )
THEN
603 IF ( lsame( jobu,
'F' ) ) n1 = m
610 epsln = slamch(
'Epsilon')
611 sfmin = slamch(
'SafeMinimum')
612 small = sfmin / epsln
622 scalem = one / sqrt(float(m)*float(n))
628 CALL slassq( m, a(1,p), 1, aapp, aaqq )
629 IF ( aapp .GT. big )
THEN
631 CALL xerbla(
'SGEJSV', -info )
635 IF ( ( aapp .LT. (big / aaqq) ) .AND. noscal )
THEN
639 sva(p) = aapp * ( aaqq * scalem )
642 CALL sscal( p-1, scalem, sva, 1 )
647 IF ( noscal ) scalem = one
652 aapp = max( aapp, sva(p) )
653 IF ( sva(p) .NE. zero ) aaqq = min( aaqq, sva(p) )
658 IF ( aapp .EQ. zero )
THEN
659 IF ( lsvec )
CALL slaset(
'G', m, n1, zero, one, u, ldu )
660 IF ( rsvec )
CALL slaset(
'G', n, n, zero, one, v, ldv )
663 IF ( errest ) work(3) = one
664 IF ( lsvec .AND. rsvec )
THEN
683 IF ( aaqq .LE. sfmin )
THEN
694 CALL slascl(
'G',0,0,sva(1),scalem, m,1,a(1,1),lda,ierr )
695 CALL slacpy(
'A', m, 1, a, lda, u, ldu )
697 IF ( n1 .NE. n )
THEN
698 CALL sgeqrf( m, n, u,ldu, work, work(n+1),lwork-n,ierr )
699 CALL sorgqr( m,n1,1, u,ldu,work,work(n+1),lwork-n,ierr )
700 CALL scopy( m, a(1,1), 1, u(1,1), 1 )
706 IF ( sva(1) .LT. (big*scalem) )
THEN
707 sva(1) = sva(1) / scalem
710 work(1) = one / scalem
712 IF ( sva(1) .NE. zero )
THEN
714 IF ( ( sva(1) / scalem) .GE. sfmin )
THEN
724 IF ( errest ) work(3) = one
725 IF ( lsvec .AND. rsvec )
THEN
738 l2tran = l2tran .AND. ( m .EQ. n )
742 IF ( rowpiv .OR. l2tran )
THEN
753 CALL slassq( n, a(p,1), lda, xsc, temp1 )
756 work(m+n+p) = xsc * scalem
757 work(n+p) = xsc * (scalem*sqrt(temp1))
758 aatmax = max( aatmax, work(n+p) )
759 IF (work(n+p) .NE. zero) aatmin = min(aatmin,work(n+p))
763 work(m+n+p) = scalem*abs( a(p,isamax(n,a(p,1),lda)) )
764 aatmax = max( aatmax, work(m+n+p) )
765 aatmin = min( aatmin, work(m+n+p) )
784 CALL slassq( n, sva, 1, xsc, temp1 )
789 big1 = ( ( sva(p) / xsc )**2 ) * temp1
790 IF ( big1 .NE. zero ) entra = entra + big1 * alog(big1)
792 entra = - entra / alog(float(n))
802 big1 = ( ( work(p) / xsc )**2 ) * temp1
803 IF ( big1 .NE. zero ) entrat = entrat + big1 * alog(big1)
805 entrat = - entrat / alog(float(m))
810 transp = ( entrat .LT. entra )
856 temp1 = sqrt( big / float(n) )
858 CALL slascl(
'G', 0, 0, aapp, temp1, n, 1, sva, n, ierr )
859 IF ( aaqq .GT. (aapp * sfmin) )
THEN
860 aaqq = ( aaqq / aapp ) * temp1
862 aaqq = ( aaqq * temp1 ) / aapp
864 temp1 = temp1 * scalem
865 CALL slascl(
'G', 0, 0, aapp, temp1, m, n, a, lda, ierr )
889 IF ( ( aaqq.LT.sqrt(sfmin) ) .AND. lsvec .AND. rsvec )
THEN
894 IF ( aaqq .LT. xsc )
THEN
896 IF ( sva(p) .LT. xsc )
THEN
897 CALL slaset(
'A', m, 1, zero, zero, a(1,p), lda )
912 q = isamax( m-p+1, work(m+n+p), 1 ) + p - 1
916 work(m+n+p) = work(m+n+q)
920 CALL slaswp( n, a, lda, 1, m-1, iwork(2*n+1), 1 )
942 CALL sgeqp3( m,n,a,lda, iwork,work, work(n+1),lwork-n, ierr )
958 temp1 = sqrt(float(n))*epsln
960 IF ( abs(a(p,p)) .GE. (temp1*abs(a(1,1))) )
THEN
967 ELSE IF ( l2rank )
THEN
973 IF ( ( abs(a(p,p)) .LT. (epsln*abs(a(p-1,p-1))) ) .OR.
974 $ ( abs(a(p,p)) .LT. small ) .OR.
975 $ ( l2kill .AND. (abs(a(p,p)) .LT. temp1) ) )
GO TO 3402
990 IF ( ( abs(a(p,p)) .LT. small ) .OR.
991 $ ( l2kill .AND. (abs(a(p,p)) .LT. temp1) ) )
GO TO 3302
999 IF ( nr .EQ. n )
THEN
1002 temp1 = abs(a(p,p)) / sva(iwork(p))
1003 maxprj = min( maxprj, temp1 )
1005 IF ( maxprj**2 .GE. one - float(n)*epsln ) almort = .true.
1014 IF ( n .EQ. nr )
THEN
1017 CALL slacpy(
'U', n, n, a, lda, v, ldv )
1019 temp1 = sva(iwork(p))
1020 CALL sscal( p, one/temp1, v(1,p), 1 )
1022 CALL spocon(
'U', n, v, ldv, one, temp1,
1023 $ work(n+1), iwork(2*n+m+1), ierr )
1024 ELSE IF ( lsvec )
THEN
1026 CALL slacpy(
'U', n, n, a, lda, u, ldu )
1028 temp1 = sva(iwork(p))
1029 CALL sscal( p, one/temp1, u(1,p), 1 )
1031 CALL spocon(
'U', n, u, ldu, one, temp1,
1032 $ work(n+1), iwork(2*n+m+1), ierr )
1034 CALL slacpy(
'U', n, n, a, lda, work(n+1), n )
1036 temp1 = sva(iwork(p))
1037 CALL sscal( p, one/temp1, work(n+(p-1)*n+1), 1 )
1040 CALL spocon(
'U', n, work(n+1), n, one, temp1,
1041 $ work(n+n*n+1), iwork(2*n+m+1), ierr )
1043 sconda = one / sqrt(temp1)
1051 l2pert = l2pert .AND. ( abs( a(1,1)/a(nr,nr) ) .GT. sqrt(big1) )
1056 IF ( .NOT. ( rsvec .OR. lsvec ) )
THEN
1061 DO 1946 p = 1, min( n-1, nr )
1062 CALL scopy( n-p, a(p,p+1), lda, a(p+1,p), 1 )
1077 IF ( .NOT. almort )
THEN
1081 xsc = epsln / float(n)
1083 temp1 = xsc*abs(a(q,q))
1085 IF ( ( (p.GT.q) .AND. (abs(a(p,q)).LE.temp1) )
1086 $ .OR. ( p .LT. q ) )
1087 $ a(p,q) = sign( temp1, a(p,q) )
1091 CALL slaset(
'U', nr-1,nr-1, zero,zero, a(1,2),lda )
1096 CALL sgeqrf( n,nr, a,lda, work, work(n+1),lwork-n, ierr )
1099 DO 1948 p = 1, nr - 1
1100 CALL scopy( nr-p, a(p,p+1), lda, a(p+1,p), 1 )
1111 xsc = epsln / float(n)
1113 temp1 = xsc*abs(a(q,q))
1115 IF ( ( (p.GT.q) .AND. (abs(a(p,q)).LE.temp1) )
1116 $ .OR. ( p .LT. q ) )
1117 $ a(p,q) = sign( temp1, a(p,q) )
1121 CALL slaset(
'U', nr-1, nr-1, zero, zero, a(1,2), lda )
1128 CALL sgesvj(
'L',
'NoU',
'NoV', nr, nr, a, lda, sva,
1129 $ n, v, ldv, work, lwork, info )
1132 numrank = nint(work(2))
1135 ELSE IF ( rsvec .AND. ( .NOT. lsvec ) )
THEN
1143 CALL scopy( n-p+1, a(p,p), lda, v(p,p), 1 )
1145 CALL slaset(
'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
1147 CALL sgesvj(
'L',
'U',
'N', n, nr, v,ldv, sva, nr, a,lda,
1148 $ work, lwork, info )
1150 numrank = nint(work(2))
1157 CALL slaset(
'Lower', nr-1, nr-1, zero, zero, a(2,1), lda )
1158 CALL sgelqf( nr, n, a, lda, work, work(n+1), lwork-n, ierr)
1159 CALL slacpy(
'Lower', nr, nr, a, lda, v, ldv )
1160 CALL slaset(
'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
1161 CALL sgeqrf( nr, nr, v, ldv, work(n+1), work(2*n+1),
1164 CALL scopy( nr-p+1, v(p,p), ldv, v(p,p), 1 )
1166 CALL slaset(
'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
1168 CALL sgesvj(
'Lower',
'U',
'N', nr, nr, v,ldv, sva, nr, u,
1169 $ ldu, work(n+1), lwork-n, info )
1171 numrank = nint(work(n+2))
1172 IF ( nr .LT. n )
THEN
1173 CALL slaset(
'A',n-nr, nr, zero,zero, v(nr+1,1), ldv )
1174 CALL slaset(
'A',nr, n-nr, zero,zero, v(1,nr+1), ldv )
1175 CALL slaset(
'A',n-nr,n-nr,zero,one, v(nr+1,nr+1), ldv )
1178 CALL sormlq(
'Left',
'Transpose', n, n, nr, a, lda, work,
1179 $ v, ldv, work(n+1), lwork-n, ierr )
1184 CALL scopy( n, v(p,1), ldv, a(iwork(p),1), lda )
1186 CALL slacpy(
'All', n, n, a, lda, v, ldv )
1189 CALL slacpy(
'All', n, n, v, ldv, u, ldu )
1192 ELSE IF ( lsvec .AND. ( .NOT. rsvec ) )
THEN
1199 CALL scopy( n-p+1, a(p,p), lda, u(p,p), 1 )
1201 CALL slaset(
'Upper', nr-1, nr-1, zero, zero, u(1,2), ldu )
1203 CALL sgeqrf( n, nr, u, ldu, work(n+1), work(2*n+1),
1206 DO 1967 p = 1, nr - 1
1207 CALL scopy( nr-p, u(p,p+1), ldu, u(p+1,p), 1 )
1209 CALL slaset(
'Upper', nr-1, nr-1, zero, zero, u(1,2), ldu )
1211 CALL sgesvj(
'Lower',
'U',
'N', nr,nr, u, ldu, sva, nr, a,
1212 $ lda, work(n+1), lwork-n, info )
1214 numrank = nint(work(n+2))
1216 IF ( nr .LT. m )
THEN
1217 CALL slaset(
'A', m-nr, nr,zero, zero, u(nr+1,1), ldu )
1218 IF ( nr .LT. n1 )
THEN
1219 CALL slaset(
'A',nr, n1-nr, zero, zero, u(1,nr+1), ldu )
1220 CALL slaset(
'A',m-nr,n1-nr,zero,one,u(nr+1,nr+1), ldu )
1224 CALL sormqr(
'Left',
'No Tr', m, n1, n, a, lda, work, u,
1225 $ ldu, work(n+1), lwork-n, ierr )
1228 $
CALL slaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1231 xsc = one / snrm2( m, u(1,p), 1 )
1232 CALL sscal( m, xsc, u(1,p), 1 )
1236 CALL slacpy(
'All', n, n, u, ldu, v, ldv )
1243 IF ( .NOT. jracc )
THEN
1245 IF ( .NOT. almort )
THEN
1255 CALL scopy( n-p+1, a(p,p), lda, v(p,p), 1 )
1273 temp1 = xsc*abs( v(q,q) )
1275 IF ( ( p .GT. q ) .AND. ( abs(v(p,q)) .LE. temp1 )
1276 $ .OR. ( p .LT. q ) )
1277 $ v(p,q) = sign( temp1, v(p,q) )
1278 IF ( p .LT. q ) v(p,q) = - v(p,q)
1282 CALL slaset(
'U', nr-1, nr-1, zero, zero, v(1,2), ldv )
1289 CALL slacpy(
'L', nr, nr, v, ldv, work(2*n+1), nr )
1291 temp1 = snrm2(nr-p+1,work(2*n+(p-1)*nr+p),1)
1292 CALL sscal(nr-p+1,one/temp1,work(2*n+(p-1)*nr+p),1)
1294 CALL spocon(
'Lower',nr,work(2*n+1),nr,one,temp1,
1295 $ work(2*n+nr*nr+1),iwork(m+2*n+1),ierr)
1296 condr1 = one / sqrt(temp1)
1302 cond_ok = sqrt(float(nr))
1305 IF ( condr1 .LT. cond_ok )
THEN
1310 CALL sgeqrf( n, nr, v, ldv, work(n+1), work(2*n+1),
1314 xsc = sqrt(small)/epsln
1316 DO 3958 q = 1, p - 1
1317 temp1 = xsc * min(abs(v(p,p)),abs(v(q,q)))
1318 IF ( abs(v(q,p)) .LE. temp1 )
1319 $ v(q,p) = sign( temp1, v(q,p) )
1325 $
CALL slacpy(
'A', n, nr, v, ldv, work(2*n+1), n )
1329 DO 1969 p = 1, nr - 1
1330 CALL scopy( nr-p, v(p,p+1), ldv, v(p+1,p), 1 )
1348 CALL sgeqp3( n, nr, v, ldv, iwork(n+1), work(n+1),
1349 $ work(2*n+1), lwork-2*n, ierr )
1355 DO 3968 q = 1, p - 1
1356 temp1 = xsc * min(abs(v(p,p)),abs(v(q,q)))
1357 IF ( abs(v(q,p)) .LE. temp1 )
1358 $ v(q,p) = sign( temp1, v(q,p) )
1363 CALL slacpy(
'A', n, nr, v, ldv, work(2*n+1), n )
1368 DO 8971 q = 1, p - 1
1369 temp1 = xsc * min(abs(v(p,p)),abs(v(q,q)))
1370 v(p,q) = - sign( temp1, v(q,p) )
1374 CALL slaset(
'L',nr-1,nr-1,zero,zero,v(2,1),ldv )
1377 CALL sgelqf( nr, nr, v, ldv, work(2*n+n*nr+1),
1378 $ work(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr, ierr )
1380 CALL slacpy(
'L',nr,nr,v,ldv,work(2*n+n*nr+nr+1),nr )
1382 temp1 = snrm2( p, work(2*n+n*nr+nr+p), nr )
1383 CALL sscal( p, one/temp1, work(2*n+n*nr+nr+p), nr )
1385 CALL spocon(
'L',nr,work(2*n+n*nr+nr+1),nr,one,temp1,
1386 $ work(2*n+n*nr+nr+nr*nr+1),iwork(m+2*n+1),ierr )
1387 condr2 = one / sqrt(temp1)
1389 IF ( condr2 .GE. cond_ok )
THEN
1394 CALL slacpy(
'U', nr, nr, v, ldv, work(2*n+1), n )
1404 temp1 = xsc * v(q,q)
1405 DO 4969 p = 1, q - 1
1407 v(p,q) = - sign( temp1, v(p,q) )
1411 CALL slaset(
'U', nr-1,nr-1, zero,zero, v(1,2), ldv )
1420 IF ( condr1 .LT. cond_ok )
THEN
1422 CALL sgesvj(
'L',
'U',
'N',nr,nr,v,ldv,sva,nr,u,
1423 $ ldu,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,info )
1424 scalem = work(2*n+n*nr+nr+1)
1425 numrank = nint(work(2*n+n*nr+nr+2))
1427 CALL scopy( nr, v(1,p), 1, u(1,p), 1 )
1428 CALL sscal( nr, sva(p), v(1,p), 1 )
1433 IF ( nr .EQ. n )
THEN
1438 CALL strsm(
'L',
'U',
'N',
'N', nr,nr,one, a,lda, v,ldv )
1444 CALL strsm(
'L',
'U',
'T',
'N',nr,nr,one,work(2*n+1),
1446 IF ( nr .LT. n )
THEN
1447 CALL slaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv)
1448 CALL slaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv)
1449 CALL slaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv)
1451 CALL sormqr(
'L',
'N',n,n,nr,work(2*n+1),n,work(n+1),
1452 $ v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr)
1455 ELSE IF ( condr2 .LT. cond_ok )
THEN
1463 CALL sgesvj(
'L',
'U',
'N', nr, nr, v, ldv, sva, nr, u,
1464 $ ldu, work(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr, info )
1465 scalem = work(2*n+n*nr+nr+1)
1466 numrank = nint(work(2*n+n*nr+nr+2))
1468 CALL scopy( nr, v(1,p), 1, u(1,p), 1 )
1469 CALL sscal( nr, sva(p), u(1,p), 1 )
1471 CALL strsm(
'L',
'U',
'N',
'N',nr,nr,one,work(2*n+1),n,u,ldu)
1475 work(2*n+n*nr+nr+iwork(n+p)) = u(p,q)
1478 u(p,q) = work(2*n+n*nr+nr+p)
1481 IF ( nr .LT. n )
THEN
1482 CALL slaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
1483 CALL slaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
1484 CALL slaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
1486 CALL sormqr(
'L',
'N',n,n,nr,work(2*n+1),n,work(n+1),
1487 $ v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
1500 CALL sgesvj(
'L',
'U',
'V', nr, nr, v, ldv, sva, nr, u,
1501 $ ldu, work(2*n+n*nr+nr+1), lwork-2*n-n*nr-nr, info )
1502 scalem = work(2*n+n*nr+nr+1)
1503 numrank = nint(work(2*n+n*nr+nr+2))
1504 IF ( nr .LT. n )
THEN
1505 CALL slaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
1506 CALL slaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
1507 CALL slaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
1509 CALL sormqr(
'L',
'N',n,n,nr,work(2*n+1),n,work(n+1),
1510 $ v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
1512 CALL sormlq(
'L',
'T', nr, nr, nr, work(2*n+1), n,
1513 $ work(2*n+n*nr+1), u, ldu, work(2*n+n*nr+nr+1),
1514 $ lwork-2*n-n*nr-nr, ierr )
1517 work(2*n+n*nr+nr+iwork(n+p)) = u(p,q)
1520 u(p,q) = work(2*n+n*nr+nr+p)
1530 temp1 = sqrt(float(n)) * epsln
1533 work(2*n+n*nr+nr+iwork(p)) = v(p,q)
1536 v(p,q) = work(2*n+n*nr+nr+p)
1538 xsc = one / snrm2( n, v(1,q), 1 )
1539 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1540 $
CALL sscal( n, xsc, v(1,q), 1 )
1544 IF ( nr .LT. m )
THEN
1545 CALL slaset(
'A', m-nr, nr, zero, zero, u(nr+1,1), ldu )
1546 IF ( nr .LT. n1 )
THEN
1547 CALL slaset(
'A',nr,n1-nr,zero,zero,u(1,nr+1),ldu)
1548 CALL slaset(
'A',m-nr,n1-nr,zero,one,u(nr+1,nr+1),ldu)
1555 CALL sormqr(
'Left',
'No_Tr', m, n1, n, a, lda, work, u,
1556 $ ldu, work(n+1), lwork-n, ierr )
1559 temp1 = sqrt(float(m)) * epsln
1561 xsc = one / snrm2( m, u(1,p), 1 )
1562 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1563 $
CALL sscal( m, xsc, u(1,p), 1 )
1570 $
CALL slaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1577 CALL slacpy(
'Upper', n, n, a, lda, work(n+1), n )
1581 temp1 = xsc * work( n + (p-1)*n + p )
1582 DO 5971 q = 1, p - 1
1583 work(n+(q-1)*n+p)=-sign(temp1,work(n+(p-1)*n+q))
1587 CALL slaset(
'Lower',n-1,n-1,zero,zero,work(n+2),n )
1590 CALL sgesvj(
'Upper',
'U',
'N', n, n, work(n+1), n, sva,
1591 $ n, u, ldu, work(n+n*n+1), lwork-n-n*n, info )
1593 scalem = work(n+n*n+1)
1594 numrank = nint(work(n+n*n+2))
1596 CALL scopy( n, work(n+(p-1)*n+1), 1, u(1,p), 1 )
1597 CALL sscal( n, sva(p), work(n+(p-1)*n+1), 1 )
1600 CALL strsm(
'Left',
'Upper',
'NoTrans',
'No UD', n, n,
1601 $ one, a, lda, work(n+1), n )
1603 CALL scopy( n, work(n+p), n, v(iwork(p),1), ldv )
1605 temp1 = sqrt(float(n))*epsln
1607 xsc = one / snrm2( n, v(1,p), 1 )
1608 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1609 $
CALL sscal( n, xsc, v(1,p), 1 )
1614 IF ( n .LT. m )
THEN
1615 CALL slaset(
'A', m-n, n, zero, zero, u(n+1,1), ldu )
1616 IF ( n .LT. n1 )
THEN
1617 CALL slaset(
'A',n, n1-n, zero, zero, u(1,n+1),ldu )
1618 CALL slaset(
'A',m-n,n1-n, zero, one,u(n+1,n+1),ldu )
1621 CALL sormqr(
'Left',
'No Tr', m, n1, n, a, lda, work, u,
1622 $ ldu, work(n+1), lwork-n, ierr )
1623 temp1 = sqrt(float(m))*epsln
1625 xsc = one / snrm2( m, u(1,p), 1 )
1626 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1627 $
CALL sscal( m, xsc, u(1,p), 1 )
1631 $
CALL slaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1650 CALL scopy( n-p+1, a(p,p), lda, v(p,p), 1 )
1654 xsc = sqrt(small/epsln)
1656 temp1 = xsc*abs( v(q,q) )
1658 IF ( ( p .GT. q ) .AND. ( abs(v(p,q)) .LE. temp1 )
1659 $ .OR. ( p .LT. q ) )
1660 $ v(p,q) = sign( temp1, v(p,q) )
1661 IF ( p .LT. q ) v(p,q) = - v(p,q)
1665 CALL slaset(
'U', nr-1, nr-1, zero, zero, v(1,2), ldv )
1668 CALL sgeqrf( n, nr, v, ldv, work(n+1), work(2*n+1),
1670 CALL slacpy(
'L', n, nr, v, ldv, work(2*n+1), n )
1673 CALL scopy( nr-p+1, v(p,p), ldv, u(p,p), 1 )
1677 xsc = sqrt(small/epsln)
1679 DO 9971 p = 1, q - 1
1680 temp1 = xsc * min(abs(u(p,p)),abs(u(q,q)))
1681 u(p,q) = - sign( temp1, u(q,p) )
1685 CALL slaset(
'U', nr-1, nr-1, zero, zero, u(1,2), ldu )
1688 CALL sgesvj(
'L',
'U',
'V', nr, nr, u, ldu, sva,
1689 $ n, v, ldv, work(2*n+n*nr+1), lwork-2*n-n*nr, info )
1690 scalem = work(2*n+n*nr+1)
1691 numrank = nint(work(2*n+n*nr+2))
1693 IF ( nr .LT. n )
THEN
1694 CALL slaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
1695 CALL slaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
1696 CALL slaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
1699 CALL sormqr(
'L',
'N',n,n,nr,work(2*n+1),n,work(n+1),
1700 $ v,ldv,work(2*n+n*nr+nr+1),lwork-2*n-n*nr-nr,ierr )
1706 temp1 = sqrt(float(n)) * epsln
1709 work(2*n+n*nr+nr+iwork(p)) = v(p,q)
1712 v(p,q) = work(2*n+n*nr+nr+p)
1714 xsc = one / snrm2( n, v(1,q), 1 )
1715 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1716 $
CALL sscal( n, xsc, v(1,q), 1 )
1722 IF ( nr .LT. m )
THEN
1723 CALL slaset(
'A', m-nr, nr, zero, zero, u(nr+1,1), ldu )
1724 IF ( nr .LT. n1 )
THEN
1725 CALL slaset(
'A',nr, n1-nr, zero, zero, u(1,nr+1),ldu )
1726 CALL slaset(
'A',m-nr,n1-nr, zero, one,u(nr+1,nr+1),ldu )
1730 CALL sormqr(
'Left',
'No Tr', m, n1, n, a, lda, work, u,
1731 $ ldu, work(n+1), lwork-n, ierr )
1734 $
CALL slaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1741 CALL sswap( n, u(1,p), 1, v(1,p), 1 )
1750 IF ( uscal2 .LE. (big/sva(1))*uscal1 )
THEN
1751 CALL slascl(
'G', 0, 0, uscal1, uscal2, nr, 1, sva, n, ierr )
1756 IF ( nr .LT. n )
THEN
1762 work(1) = uscal2 * scalem
1764 IF ( errest ) work(3) = sconda
1765 IF ( lsvec .AND. rsvec )
THEN
subroutine sgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
SGESVJ
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 sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
subroutine sgejsv(JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, M, N, A, LDA, SVA, U, LDU, V, LDV, WORK, LWORK, IWORK, INFO)
SGEJSV
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMLQ
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 spocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
SPOCON
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 sgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
SGEQP3
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine sgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGELQF
subroutine slassq(N, X, INCX, SCALE, SUMSQ)
SLASSQ updates a sum of squares represented in scaled form.