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
513 DOUBLE PRECISION DLAMCH, DNRM2
516 EXTERNAL idamax, lsame, dlamch, dnrm2
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
605 epsln = dlamch(
'Epsilon')
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
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.
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
subroutine dgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
DGEQP3