473 SUBROUTINE dgejsv( 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 DOUBLE PRECISION a( lda, * ), sva( n ), u( ldu, * ), v( ldv, * ),
490 CHARACTER*1 joba, jobp, jobr, jobt, jobu, jobv
496 DOUBLE PRECISION zero, one
497 parameter( zero = 0.0d0, one = 1.0d0 )
500 DOUBLE PRECISION 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 dabs, dlog, dmax1, dmin1, dble,
510 $ max0, min0, idnint, dsign, dsqrt
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(
'DGEJSV', - info )
592 IF ( ( m .EQ. 0 ) .OR. ( n .EQ. 0 ) )
RETURN
598 IF (
lsame( jobu,
'F' ) ) n1 = m
606 sfmin =
dlamch(
'SafeMinimum')
607 small = sfmin / epsln
617 scalem = one / dsqrt(dble(m)*dble(n))
623 CALL
dlassq( m, a(1,p), 1, aapp, aaqq )
624 IF ( aapp .GT. big )
THEN
626 CALL
xerbla(
'DGEJSV', -info )
630 IF ( ( aapp .LT. (big / aaqq) ) .AND. noscal )
THEN
634 sva(p) = aapp * ( aaqq * scalem )
637 CALL
dscal( p-1, scalem, sva, 1 )
642 IF ( noscal ) scalem = one
647 aapp = dmax1( aapp, sva(p) )
648 IF ( sva(p) .NE. zero ) aaqq = dmin1( aaqq, sva(p) )
653 IF ( aapp .EQ. zero )
THEN
654 IF ( lsvec ) CALL
dlaset(
'G', m, n1, zero, one, u, ldu )
655 IF ( rsvec ) CALL
dlaset(
'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
dlascl(
'G',0,0,sva(1),scalem, m,1,a(1,1),lda,ierr )
690 CALL
dlacpy(
'A', m, 1, a, lda, u, ldu )
692 IF ( n1 .NE. n )
THEN
693 CALL
dgeqrf( m, n, u,ldu, work, work(n+1),lwork-n,ierr )
694 CALL
dorgqr( m,n1,1, u,ldu,work,work(n+1),lwork-n,ierr )
695 CALL
dcopy( 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
dlassq( n, a(p,1), lda, xsc, temp1 )
750 work(m+n+p) = xsc * scalem
751 work(n+p) = xsc * (scalem*dsqrt(temp1))
752 aatmax = dmax1( aatmax, work(n+p) )
753 IF (work(n+p) .NE. zero) aatmin = dmin1(aatmin,work(n+p))
757 work(m+n+p) = scalem*dabs( a(p,
idamax(n,a(p,1),lda)) )
758 aatmax = dmax1( aatmax, work(m+n+p) )
759 aatmin = dmin1( aatmin, work(m+n+p) )
778 CALL
dlassq( n, sva, 1, xsc, temp1 )
783 big1 = ( ( sva(p) / xsc )**2 ) * temp1
784 IF ( big1 .NE. zero ) entra = entra + big1 * dlog(big1)
786 entra = - entra / dlog(dble(n))
796 big1 = ( ( work(p) / xsc )**2 ) * temp1
797 IF ( big1 .NE. zero ) entrat = entrat + big1 * dlog(big1)
799 entrat = - entrat / dlog(dble(m))
804 transp = ( entrat .LT. entra )
850 temp1 = dsqrt( big / dble(n) )
852 CALL
dlascl(
'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
dlascl(
'G', 0, 0, aapp, temp1, m, n, a, lda, ierr )
883 IF ( ( aaqq.LT.dsqrt(sfmin) ) .AND. lsvec .AND. rsvec )
THEN
888 IF ( aaqq .LT. xsc )
THEN
890 IF ( sva(p) .LT. xsc )
THEN
891 CALL
dlaset(
'A', m, 1, zero, zero, a(1,p), lda )
906 q =
idamax( m-p+1, work(m+n+p), 1 ) + p - 1
910 work(m+n+p) = work(m+n+q)
914 CALL
dlaswp( n, a, lda, 1, m-1, iwork(2*n+1), 1 )
936 CALL
dgeqp3( m,n,a,lda, iwork,work, work(n+1),lwork-n, ierr )
952 temp1 = dsqrt(dble(n))*epsln
954 IF ( dabs(a(p,p)) .GE. (temp1*dabs(a(1,1))) )
THEN
961 ELSE IF ( l2rank )
THEN
967 IF ( ( dabs(a(p,p)) .LT. (epsln*dabs(a(p-1,p-1))) ) .OR.
968 $ ( dabs(a(p,p)) .LT. small ) .OR.
969 $ ( l2kill .AND. (dabs(a(p,p)) .LT. temp1) ) ) go to 3402
984 IF ( ( dabs(a(p,p)) .LT. small ) .OR.
985 $ ( l2kill .AND. (dabs(a(p,p)) .LT. temp1) ) ) go to 3302
993 IF ( nr .EQ. n )
THEN
996 temp1 = dabs(a(p,p)) / sva(iwork(p))
997 maxprj = dmin1( maxprj, temp1 )
999 IF ( maxprj**2 .GE. one - dble(n)*epsln ) almort = .true.
1008 IF ( n .EQ. nr )
THEN
1011 CALL
dlacpy(
'U', n, n, a, lda, v, ldv )
1013 temp1 = sva(iwork(p))
1014 CALL
dscal( p, one/temp1, v(1,p), 1 )
1016 CALL
dpocon(
'U', n, v, ldv, one, temp1,
1017 $ work(n+1), iwork(2*n+m+1), ierr )
1018 ELSE IF ( lsvec )
THEN
1020 CALL
dlacpy(
'U', n, n, a, lda, u, ldu )
1022 temp1 = sva(iwork(p))
1023 CALL
dscal( p, one/temp1, u(1,p), 1 )
1025 CALL
dpocon(
'U', n, u, ldu, one, temp1,
1026 $ work(n+1), iwork(2*n+m+1), ierr )
1028 CALL
dlacpy(
'U', n, n, a, lda, work(n+1), n )
1030 temp1 = sva(iwork(p))
1031 CALL
dscal( p, one/temp1, work(n+(p-1)*n+1), 1 )
1034 CALL
dpocon(
'U', n, work(n+1), n, one, temp1,
1035 $ work(n+n*n+1), iwork(2*n+m+1), ierr )
1037 sconda = one / dsqrt(temp1)
1045 l2pert = l2pert .AND. ( dabs( a(1,1)/a(nr,nr) ) .GT. dsqrt(big1) )
1050 IF ( .NOT. ( rsvec .OR. lsvec ) )
THEN
1055 DO 1946 p = 1, min0( n-1, nr )
1056 CALL
dcopy( n-p, a(p,p+1), lda, a(p+1,p), 1 )
1071 IF ( .NOT. almort )
THEN
1075 xsc = epsln / dble(n)
1077 temp1 = xsc*dabs(a(q,q))
1079 IF ( ( (p.GT.q) .AND. (dabs(a(p,q)).LE.temp1) )
1080 $ .OR. ( p .LT. q ) )
1081 $ a(p,q) = dsign( temp1, a(p,q) )
1085 CALL
dlaset(
'U', nr-1,nr-1, zero,zero, a(1,2),lda )
1090 CALL
dgeqrf( n,nr, a,lda, work, work(n+1),lwork-n, ierr )
1093 DO 1948 p = 1, nr - 1
1094 CALL
dcopy( nr-p, a(p,p+1), lda, a(p+1,p), 1 )
1105 xsc = epsln / dble(n)
1107 temp1 = xsc*dabs(a(q,q))
1109 IF ( ( (p.GT.q) .AND. (dabs(a(p,q)).LE.temp1) )
1110 $ .OR. ( p .LT. q ) )
1111 $ a(p,q) = dsign( temp1, a(p,q) )
1115 CALL
dlaset(
'U', nr-1, nr-1, zero, zero, a(1,2), lda )
1122 CALL
dgesvj(
'L',
'NoU',
'NoV', nr, nr, a, lda, sva,
1123 $ n, v, ldv, work, lwork, info )
1126 numrank = idnint(work(2))
1129 ELSE IF ( rsvec .AND. ( .NOT. lsvec ) )
THEN
1137 CALL
dcopy( n-p+1, a(p,p), lda, v(p,p), 1 )
1139 CALL
dlaset(
'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
1141 CALL
dgesvj(
'L',
'U',
'N', n, nr, v,ldv, sva, nr, a,lda,
1142 $ work, lwork, info )
1144 numrank = idnint(work(2))
1151 CALL
dlaset(
'Lower', nr-1, nr-1, zero, zero, a(2,1), lda )
1152 CALL
dgelqf( nr, n, a, lda, work, work(n+1), lwork-n, ierr)
1153 CALL
dlacpy(
'Lower', nr, nr, a, lda, v, ldv )
1154 CALL
dlaset(
'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
1155 CALL
dgeqrf( nr, nr, v, ldv, work(n+1), work(2*n+1),
1158 CALL
dcopy( nr-p+1, v(p,p), ldv, v(p,p), 1 )
1160 CALL
dlaset(
'Upper', nr-1, nr-1, zero, zero, v(1,2), ldv )
1162 CALL
dgesvj(
'Lower',
'U',
'N', nr, nr, v,ldv, sva, nr, u,
1163 $ ldu, work(n+1), lwork, info )
1165 numrank = idnint(work(n+2))
1166 IF ( nr .LT. n )
THEN
1167 CALL
dlaset(
'A',n-nr, nr, zero,zero, v(nr+1,1), ldv )
1168 CALL
dlaset(
'A',nr, n-nr, zero,zero, v(1,nr+1), ldv )
1169 CALL
dlaset(
'A',n-nr,n-nr,zero,one, v(nr+1,nr+1), ldv )
1172 CALL
dormlq(
'Left',
'Transpose', n, n, nr, a, lda, work,
1173 $ v, ldv, work(n+1), lwork-n, ierr )
1178 CALL
dcopy( n, v(p,1), ldv, a(iwork(p),1), lda )
1180 CALL
dlacpy(
'All', n, n, a, lda, v, ldv )
1183 CALL
dlacpy(
'All', n, n, v, ldv, u, ldu )
1186 ELSE IF ( lsvec .AND. ( .NOT. rsvec ) )
THEN
1193 CALL
dcopy( n-p+1, a(p,p), lda, u(p,p), 1 )
1195 CALL
dlaset(
'Upper', nr-1, nr-1, zero, zero, u(1,2), ldu )
1197 CALL
dgeqrf( n, nr, u, ldu, work(n+1), work(2*n+1),
1200 DO 1967 p = 1, nr - 1
1201 CALL
dcopy( nr-p, u(p,p+1), ldu, u(p+1,p), 1 )
1203 CALL
dlaset(
'Upper', nr-1, nr-1, zero, zero, u(1,2), ldu )
1205 CALL
dgesvj(
'Lower',
'U',
'N', nr,nr, u, ldu, sva, nr, a,
1206 $ lda, work(n+1), lwork-n, info )
1208 numrank = idnint(work(n+2))
1210 IF ( nr .LT. m )
THEN
1211 CALL
dlaset(
'A', m-nr, nr,zero, zero, u(nr+1,1), ldu )
1212 IF ( nr .LT. n1 )
THEN
1213 CALL
dlaset(
'A',nr, n1-nr, zero, zero, u(1,nr+1), ldu )
1214 CALL
dlaset(
'A',m-nr,n1-nr,zero,one,u(nr+1,nr+1), ldu )
1218 CALL
dormqr(
'Left',
'No Tr', m, n1, n, a, lda, work, u,
1219 $ ldu, work(n+1), lwork-n, ierr )
1222 $ CALL
dlaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1225 xsc = one /
dnrm2( m, u(1,p), 1 )
1226 CALL
dscal( m, xsc, u(1,p), 1 )
1230 CALL
dlacpy(
'All', n, n, u, ldu, v, ldv )
1237 IF ( .NOT. jracc )
THEN
1239 IF ( .NOT. almort )
THEN
1249 CALL
dcopy( n-p+1, a(p,p), lda, v(p,p), 1 )
1267 temp1 = xsc*dabs( v(q,q) )
1269 IF ( ( p .GT. q ) .AND. ( dabs(v(p,q)) .LE. temp1 )
1270 $ .OR. ( p .LT. q ) )
1271 $ v(p,q) = dsign( temp1, v(p,q) )
1272 IF ( p .LT. q ) v(p,q) = - v(p,q)
1276 CALL
dlaset(
'U', nr-1, nr-1, zero, zero, v(1,2), ldv )
1283 CALL
dlacpy(
'L', nr, nr, v, ldv, work(2*n+1), nr )
1285 temp1 =
dnrm2(nr-p+1,work(2*n+(p-1)*nr+p),1)
1286 CALL
dscal(nr-p+1,one/temp1,work(2*n+(p-1)*nr+p),1)
1288 CALL
dpocon(
'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 / dsqrt(temp1)
1296 cond_ok = dsqrt(dble(nr))
1299 IF ( condr1 .LT. cond_ok )
THEN
1304 CALL
dgeqrf( n, nr, v, ldv, work(n+1), work(2*n+1),
1308 xsc = dsqrt(small)/epsln
1310 DO 3958 q = 1, p - 1
1311 temp1 = xsc * dmin1(dabs(v(p,p)),dabs(v(q,q)))
1312 IF ( dabs(v(q,p)) .LE. temp1 )
1313 $ v(q,p) = dsign( temp1, v(q,p) )
1319 $ CALL
dlacpy(
'A', n, nr, v, ldv, work(2*n+1), n )
1323 DO 1969 p = 1, nr - 1
1324 CALL
dcopy( nr-p, v(p,p+1), ldv, v(p+1,p), 1 )
1342 CALL
dgeqp3( 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 * dmin1(dabs(v(p,p)),dabs(v(q,q)))
1351 IF ( dabs(v(q,p)) .LE. temp1 )
1352 $ v(q,p) = dsign( temp1, v(q,p) )
1357 CALL
dlacpy(
'A', n, nr, v, ldv, work(2*n+1), n )
1362 DO 8971 q = 1, p - 1
1363 temp1 = xsc * dmin1(dabs(v(p,p)),dabs(v(q,q)))
1364 v(p,q) = - dsign( temp1, v(q,p) )
1368 CALL
dlaset(
'L',nr-1,nr-1,zero,zero,v(2,1),ldv )
1371 CALL
dgelqf( 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
dlacpy(
'L',nr,nr,v,ldv,work(2*n+n*nr+nr+1),nr )
1376 temp1 =
dnrm2( p, work(2*n+n*nr+nr+p), nr )
1377 CALL
dscal( p, one/temp1, work(2*n+n*nr+nr+p), nr )
1379 CALL
dpocon(
'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 / dsqrt(temp1)
1383 IF ( condr2 .GE. cond_ok )
THEN
1388 CALL
dlacpy(
'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) = - dsign( temp1, v(p,q) )
1405 CALL
dlaset(
'U', nr-1,nr-1, zero,zero, v(1,2), ldv )
1414 IF ( condr1 .LT. cond_ok )
THEN
1416 CALL
dgesvj(
'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 = idnint(work(2*n+n*nr+nr+2))
1421 CALL
dcopy( nr, v(1,p), 1, u(1,p), 1 )
1422 CALL
dscal( nr, sva(p), v(1,p), 1 )
1427 IF ( nr .EQ. n )
THEN
1432 CALL
dtrsm(
'L',
'U',
'N',
'N', nr,nr,one, a,lda, v,ldv )
1438 CALL
dtrsm(
'L',
'U',
'T',
'N',nr,nr,one,work(2*n+1),
1440 IF ( nr .LT. n )
THEN
1441 CALL
dlaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv)
1442 CALL
dlaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv)
1443 CALL
dlaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv)
1445 CALL
dormqr(
'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
dgesvj(
'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 = idnint(work(2*n+n*nr+nr+2))
1462 CALL
dcopy( nr, v(1,p), 1, u(1,p), 1 )
1463 CALL
dscal( nr, sva(p), u(1,p), 1 )
1465 CALL
dtrsm(
'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
dlaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
1477 CALL
dlaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
1478 CALL
dlaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
1480 CALL
dormqr(
'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
dgesvj(
'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 = idnint(work(2*n+n*nr+nr+2))
1498 IF ( nr .LT. n )
THEN
1499 CALL
dlaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
1500 CALL
dlaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
1501 CALL
dlaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
1503 CALL
dormqr(
'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
dormlq(
'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 = dsqrt(dble(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 /
dnrm2( n, v(1,q), 1 )
1533 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1534 $ CALL
dscal( n, xsc, v(1,q), 1 )
1538 IF ( nr .LT. m )
THEN
1539 CALL
dlaset(
'A', m-nr, nr, zero, zero, u(nr+1,1), ldu )
1540 IF ( nr .LT. n1 )
THEN
1541 CALL
dlaset(
'A',nr,n1-nr,zero,zero,u(1,nr+1),ldu)
1542 CALL
dlaset(
'A',m-nr,n1-nr,zero,one,u(nr+1,nr+1),ldu)
1549 CALL
dormqr(
'Left',
'No_Tr', m, n1, n, a, lda, work, u,
1550 $ ldu, work(n+1), lwork-n, ierr )
1553 temp1 = dsqrt(dble(m)) * epsln
1555 xsc = one /
dnrm2( m, u(1,p), 1 )
1556 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1557 $ CALL
dscal( m, xsc, u(1,p), 1 )
1564 $ CALL
dlaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1571 CALL
dlacpy(
'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)=-dsign(temp1,work(n+(p-1)*n+q))
1581 CALL
dlaset(
'Lower',n-1,n-1,zero,zero,work(n+2),n )
1584 CALL
dgesvj(
'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 = idnint(work(n+n*n+2))
1590 CALL
dcopy( n, work(n+(p-1)*n+1), 1, u(1,p), 1 )
1591 CALL
dscal( n, sva(p), work(n+(p-1)*n+1), 1 )
1594 CALL
dtrsm(
'Left',
'Upper',
'NoTrans',
'No UD', n, n,
1595 $ one, a, lda, work(n+1), n )
1597 CALL
dcopy( n, work(n+p), n, v(iwork(p),1), ldv )
1599 temp1 = dsqrt(dble(n))*epsln
1601 xsc = one /
dnrm2( n, v(1,p), 1 )
1602 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1603 $ CALL
dscal( n, xsc, v(1,p), 1 )
1608 IF ( n .LT. m )
THEN
1609 CALL
dlaset(
'A', m-n, n, zero, zero, u(n+1,1), ldu )
1610 IF ( n .LT. n1 )
THEN
1611 CALL
dlaset(
'A',n, n1-n, zero, zero, u(1,n+1),ldu )
1612 CALL
dlaset(
'A',m-n,n1-n, zero, one,u(n+1,n+1),ldu )
1615 CALL
dormqr(
'Left',
'No Tr', m, n1, n, a, lda, work, u,
1616 $ ldu, work(n+1), lwork-n, ierr )
1617 temp1 = dsqrt(dble(m))*epsln
1619 xsc = one /
dnrm2( m, u(1,p), 1 )
1620 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1621 $ CALL
dscal( m, xsc, u(1,p), 1 )
1625 $ CALL
dlaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1644 CALL
dcopy( n-p+1, a(p,p), lda, v(p,p), 1 )
1648 xsc = dsqrt(small/epsln)
1650 temp1 = xsc*dabs( v(q,q) )
1652 IF ( ( p .GT. q ) .AND. ( dabs(v(p,q)) .LE. temp1 )
1653 $ .OR. ( p .LT. q ) )
1654 $ v(p,q) = dsign( temp1, v(p,q) )
1655 IF ( p .LT. q ) v(p,q) = - v(p,q)
1659 CALL
dlaset(
'U', nr-1, nr-1, zero, zero, v(1,2), ldv )
1662 CALL
dgeqrf( n, nr, v, ldv, work(n+1), work(2*n+1),
1664 CALL
dlacpy(
'L', n, nr, v, ldv, work(2*n+1), n )
1667 CALL
dcopy( nr-p+1, v(p,p), ldv, u(p,p), 1 )
1671 xsc = dsqrt(small/epsln)
1673 DO 9971 p = 1, q - 1
1674 temp1 = xsc * dmin1(dabs(u(p,p)),dabs(u(q,q)))
1675 u(p,q) = - dsign( temp1, u(q,p) )
1679 CALL
dlaset(
'U', nr-1, nr-1, zero, zero, u(1,2), ldu )
1682 CALL
dgesvj(
'G',
'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 = idnint(work(2*n+n*nr+2))
1687 IF ( nr .LT. n )
THEN
1688 CALL
dlaset(
'A',n-nr,nr,zero,zero,v(nr+1,1),ldv )
1689 CALL
dlaset(
'A',nr,n-nr,zero,zero,v(1,nr+1),ldv )
1690 CALL
dlaset(
'A',n-nr,n-nr,zero,one,v(nr+1,nr+1),ldv )
1693 CALL
dormqr(
'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 = dsqrt(dble(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 /
dnrm2( n, v(1,q), 1 )
1709 IF ( (xsc .LT. (one-temp1)) .OR. (xsc .GT. (one+temp1)) )
1710 $ CALL
dscal( n, xsc, v(1,q), 1 )
1716 IF ( nr .LT. m )
THEN
1717 CALL
dlaset(
'A', m-nr, nr, zero, zero, u(nr+1,1), ldu )
1718 IF ( nr .LT. n1 )
THEN
1719 CALL
dlaset(
'A',nr, n1-nr, zero, zero, u(1,nr+1),ldu )
1720 CALL
dlaset(
'A',m-nr,n1-nr, zero, one,u(nr+1,nr+1),ldu )
1724 CALL
dormqr(
'Left',
'No Tr', m, n1, n, a, lda, work, u,
1725 $ ldu, work(n+1), lwork-n, ierr )
1728 $ CALL
dlaswp( n1, u, ldu, 1, m-1, iwork(2*n+1), -1 )
1735 CALL
dswap( n, u(1,p), 1, v(1,p), 1 )
1744 IF ( uscal2 .LE. (big/sva(1))*uscal1 )
THEN
1745 CALL
dlascl(
'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 dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
LOGICAL function lsame(CA, CB)
LSAME
subroutine dgejsv(JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, M, N, A, LDA, SVA, U, LDU, V, LDV, WORK, LWORK, IWORK, INFO)
DGEJSV
subroutine dgesvj(JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, LDV, WORK, LWORK, INFO)
DGESVJ
subroutine dlaswp(N, A, LDA, K1, K2, IPIV, INCX)
DLASWP performs a series of row interchanges on a general rectangular matrix.
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
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.
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
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 dpocon(UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO)
DPOCON
DOUBLE PRECISION function dnrm2(N, X, INCX)
DNRM2
subroutine dgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
DGEQP3
INTEGER function idamax(N, DX, INCX)
IDAMAX