version 1.5, 2010/12/21 13:53:25
|
version 1.6, 2011/07/22 07:38:04
|
Line 1
|
Line 1
|
SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, |
SUBROUTINE DGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, |
& M, N, A, LDA, SVA, U, LDU, V, LDV, |
$ M, N, A, LDA, SVA, U, LDU, V, LDV, |
& WORK, LWORK, IWORK, INFO ) |
$ WORK, LWORK, IWORK, INFO ) |
* |
* |
* -- LAPACK routine (version 3.3.0) -- |
* -- LAPACK routine (version 3.3.1) -- |
* |
* |
* -- Contributed by Zlatko Drmac of the University of Zagreb and -- |
* -- Contributed by Zlatko Drmac of the University of Zagreb and -- |
* -- Kresimir Veselic of the Fernuniversitaet Hagen -- |
* -- Kresimir Veselic of the Fernuniversitaet Hagen -- |
* November 2010 |
* -- April 2011 -- |
* |
* |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
Line 22
|
Line 22
|
* .. |
* .. |
* .. Array Arguments .. |
* .. Array Arguments .. |
DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ), |
DOUBLE PRECISION A( LDA, * ), SVA( N ), U( LDU, * ), V( LDV, * ), |
& WORK( LWORK ) |
$ WORK( LWORK ) |
INTEGER IWORK( * ) |
INTEGER IWORK( * ) |
CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV |
CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV |
* .. |
* .. |
Line 213
|
Line 213
|
* If JOBV = 'V' or 'J' or 'W', then LDV >= N. |
* If JOBV = 'V' or 'J' or 'W', then LDV >= N. |
* |
* |
* WORK (workspace/output) DOUBLE PRECISION array, dimension at least LWORK. |
* WORK (workspace/output) DOUBLE PRECISION array, dimension at least LWORK. |
* On exit, |
* On exit, if N.GT.0 .AND. M.GT.0 (else not referenced), |
* WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such |
* WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such |
* that SCALE*SVA(1:N) are the computed singular values |
* that SCALE*SVA(1:N) are the computed singular values |
* of A. (See the description of SVA().) |
* of A. (See the description of SVA().) |
Line 252
|
Line 252
|
* LWORK depends on the job: |
* LWORK depends on the job: |
* |
* |
* If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and |
* If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and |
* -> .. no scaled condition estimate required ( JOBE.EQ.'N'): |
* -> .. no scaled condition estimate required (JOBE.EQ.'N'): |
* LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement. |
* LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement. |
* For optimal performance (blocked code) the optimal value |
* ->> For optimal performance (blocked code) the optimal value |
* is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal |
* is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal |
* block size for xGEQP3/xGEQRF. |
* block size for DGEQP3 and DGEQRF. |
|
* In general, optimal LWORK is computed as |
|
* LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 7). |
* -> .. an estimate of the scaled condition number of A is |
* -> .. an estimate of the scaled condition number of A is |
* required (JOBA='E', 'G'). In this case, LWORK is the maximum |
* required (JOBA='E', 'G'). In this case, LWORK is the maximum |
* of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4N,7). |
* of the above and N*N+4*N, i.e. LWORK >= max(2*M+N,N*N+4*N,7). |
|
* ->> For optimal performance (blocked code) the optimal value |
|
* is LWORK >= max(2*M+N,3*N+(N+1)*NB, N*N+4*N, 7). |
|
* In general, the optimal length LWORK is computed as |
|
* LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), |
|
* N+N*N+LWORK(DPOCON),7). |
* |
* |
* If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'), |
* If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'), |
* -> the minimal requirement is LWORK >= max(2*N+M,7). |
* -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7). |
* -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7), |
* -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7), |
* where NB is the optimal block size. |
* where NB is the optimal block size for DGEQP3, DGEQRF, DGELQ, |
|
* DORMLQ. In general, the optimal length LWORK is computed as |
|
* LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON), |
|
* N+LWORK(DGELQ), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)). |
* |
* |
* If SIGMA and the left singular vectors are needed |
* If SIGMA and the left singular vectors are needed |
* -> the minimal requirement is LWORK >= max(2*N+M,7). |
* -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7). |
* -> For optimal performance, LWORK >= max(2*N+M,2*N+N*NB,7), |
* -> For optimal performance: |
* where NB is the optimal block size. |
* if JOBU.EQ.'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7), |
* |
* if JOBU.EQ.'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7), |
* If full SVD is needed ( JOBU.EQ.'U' or 'F', JOBV.EQ.'V' ) and |
* where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR. |
* -> .. the singular vectors are computed without explicit |
* In general, the optimal length LWORK is computed as |
* accumulation of the Jacobi rotations, LWORK >= 6*N+2*N*N |
* LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON), |
* -> .. in the iterative part, the Jacobi rotations are |
* 2*N+LWORK(DGEQRF), N+LWORK(DORMQR)). |
* explicitly accumulated (option, see the description of JOBV), |
* Here LWORK(DORMQR) equals N*NB (for JOBU.EQ.'U') or |
* then the minimal requirement is LWORK >= max(M+3*N+N*N,7). |
* M*NB (for JOBU.EQ.'F'). |
* For better performance, if NB is the optimal block size, |
* |
* LWORK >= max(3*N+N*N+M,3*N+N*N+N*NB,7). |
* If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and |
|
* -> if JOBV.EQ.'V' |
|
* the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N). |
|
* -> if JOBV.EQ.'J' the minimal requirement is |
|
* LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6). |
|
* -> For optimal performance, LWORK should be additionally |
|
* larger than N+M*NB, where NB is the optimal block size |
|
* for DORMQR. |
* |
* |
* IWORK (workspace/output) INTEGER array, dimension M+3*N. |
* IWORK (workspace/output) INTEGER array, dimension M+3*N. |
* On exit, |
* On exit, |
Line 300
|
Line 317
|
* Further Details |
* Further Details |
* =============== |
* =============== |
* |
* |
* DGEJSV implements a preconditioned Jacobi SVD algorithm. It uses SGEQP3, |
* DGEJSV implements a preconditioned Jacobi SVD algorithm. It uses DGEQP3, |
* SGEQRF, and SGELQF as preprocessors and preconditioners. Optionally, an |
* DGEQRF, and DGELQF as preprocessors and preconditioners. Optionally, an |
* additional row pivoting can be used as a preprocessor, which in some |
* additional row pivoting can be used as a preprocessor, which in some |
* cases results in much higher accuracy. An example is matrix A with the |
* cases results in much higher accuracy. An example is matrix A with the |
* structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned |
* structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned |
Line 319
|
Line 336
|
* the singular values in the range of normalized IEEE numbers is that the |
* the singular values in the range of normalized IEEE numbers is that the |
* spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not |
* spectral condition number kappa(A)=sigma_max(A)/sigma_min(A) does not |
* overflow. This code (DGEJSV) is best used in this restricted range, |
* overflow. This code (DGEJSV) is best used in this restricted range, |
* meaning that singular values of magnitude below ||A||_2 / SLAMCH('O') are |
* meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are |
* returned as zeros. See JOBR for details on this. |
* returned as zeros. See JOBR for details on this. |
* Further, this implementation is somewhat slower than the one described |
* Further, this implementation is somewhat slower than the one described |
* in [1,2] due to replacement of some non-LAPACK components, and because |
* in [1,2] due to replacement of some non-LAPACK components, and because |
* the choice of some tuning parameters in the iterative part (DGESVJ) is |
* the choice of some tuning parameters in the iterative part (DGESVJ) is |
* left to the implementer on a particular machine. |
* left to the implementer on a particular machine. |
* The rank revealing QR factorization (in this code: SGEQP3) should be |
* The rank revealing QR factorization (in this code: DGEQP3) should be |
* implemented as in [3]. We have a new version of SGEQP3 under development |
* implemented as in [3]. We have a new version of DGEQP3 under development |
* that is more robust than the current one in LAPACK, with a cleaner cut in |
* that is more robust than the current one in LAPACK, with a cleaner cut in |
* rank defficient cases. It will be available in the SIGMA library [4]. |
* rank defficient cases. It will be available in the SIGMA library [4]. |
* If M is much larger than N, it is obvious that the inital QRF with |
* If M is much larger than N, it is obvious that the inital QRF with |
Line 360
|
Line 377
|
* Department of Mathematics, University of Zagreb, 2008. |
* Department of Mathematics, University of Zagreb, 2008. |
* |
* |
* Bugs, examples and comments |
* Bugs, examples and comments |
* |
* |
* Please report all bugs and send interesting examples and/or comments to |
* Please report all bugs and send interesting examples and/or comments to |
* drmac@math.hr. Thank you. |
* drmac@math.hr. Thank you. |
* |
* |
* ========================================================================== |
* =========================================================================== |
* |
* |
* .. Local Parameters .. |
* .. Local Parameters .. |
DOUBLE PRECISION ZERO, ONE |
DOUBLE PRECISION ZERO, ONE |
Line 372
|
Line 389
|
* .. |
* .. |
* .. Local Scalars .. |
* .. Local Scalars .. |
DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK, |
DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, COND_OK, |
& CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN, MAXPRJ, SCALEM, |
$ CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN, MAXPRJ, SCALEM, |
& SCONDA, SFMIN, SMALL, TEMP1, USCAL1, USCAL2, XSC |
$ SCONDA, SFMIN, SMALL, TEMP1, USCAL1, USCAL2, XSC |
INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING |
INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING |
LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LSVEC, |
LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LSVEC, |
& L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN, |
$ L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN, |
& NOSCAL, ROWPIV, RSVEC, TRANSP |
$ NOSCAL, ROWPIV, RSVEC, TRANSP |
* .. |
* .. |
* .. Intrinsic Functions .. |
* .. Intrinsic Functions .. |
INTRINSIC DABS, DLOG, DMAX1, DMIN1, DBLE, |
INTRINSIC DABS, DLOG, DMAX1, DMIN1, DBLE, |
& MAX0, MIN0, IDNINT, DSIGN, DSQRT |
$ MAX0, MIN0, IDNINT, DSIGN, DSQRT |
* .. |
* .. |
* .. External Functions .. |
* .. External Functions .. |
DOUBLE PRECISION DLAMCH, DNRM2 |
DOUBLE PRECISION DLAMCH, DNRM2 |
Line 391
|
Line 408
|
* .. |
* .. |
* .. External Subroutines .. |
* .. External Subroutines .. |
EXTERNAL DCOPY, DGELQF, DGEQP3, DGEQRF, DLACPY, DLASCL, |
EXTERNAL DCOPY, DGELQF, DGEQP3, DGEQRF, DLACPY, DLASCL, |
& DLASET, DLASSQ, DLASWP, DORGQR, DORMLQ, |
$ DLASET, DLASSQ, DLASWP, DORGQR, DORMLQ, |
& DORMQR, DPOCON, DSCAL, DSWAP, DTRSM, XERBLA |
$ DORMQR, DPOCON, DSCAL, DSWAP, DTRSM, XERBLA |
* |
* |
EXTERNAL DGESVJ |
EXTERNAL DGESVJ |
* .. |
* .. |
Line 412
|
Line 429
|
L2PERT = LSAME( JOBP, 'P' ) |
L2PERT = LSAME( JOBP, 'P' ) |
* |
* |
IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR. |
IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR. |
& ERREST .OR. LSAME( JOBA, 'C' ) )) THEN |
$ ERREST .OR. LSAME( JOBA, 'C' ) )) THEN |
INFO = - 1 |
INFO = - 1 |
ELSE IF ( .NOT.( LSVEC .OR. LSAME( JOBU, 'N' ) .OR. |
ELSE IF ( .NOT.( LSVEC .OR. LSAME( JOBU, 'N' ) .OR. |
& LSAME( JOBU, 'W' )) ) THEN |
$ LSAME( JOBU, 'W' )) ) THEN |
INFO = - 2 |
INFO = - 2 |
ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR. |
ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR. |
& LSAME( JOBV, 'W' )) .OR. ( JRACC .AND. (.NOT.LSVEC) ) ) THEN |
$ LSAME( JOBV, 'W' )) .OR. ( JRACC .AND. (.NOT.LSVEC) ) ) THEN |
INFO = - 3 |
INFO = - 3 |
ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) ) THEN |
ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) ) THEN |
INFO = - 4 |
INFO = - 4 |
Line 438
|
Line 455
|
INFO = - 14 |
INFO = - 14 |
ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND. |
ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND. |
& (LWORK .LT. MAX0(7,4*N+1,2*M+N))) .OR. |
& (LWORK .LT. MAX0(7,4*N+1,2*M+N))) .OR. |
& (.NOT.(LSVEC .OR. LSVEC) .AND. ERREST .AND. |
& (.NOT.(LSVEC .OR. RSVEC) .AND. ERREST .AND. |
& (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR. |
& (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR. |
& (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR. |
& (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*M+N,4*N+1))) |
& (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*N+M))) .OR. |
& .OR. |
& (LSVEC .AND. RSVEC .AND. .NOT.JRACC .AND. (LWORK.LT.6*N+2*N*N)) |
& (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(7,2*M+N,4*N+1))) |
& .OR. (LSVEC.AND.RSVEC.AND.JRACC.AND.LWORK.LT.MAX0(7,M+3*N+N*N))) |
& .OR. |
|
& (LSVEC .AND. RSVEC .AND. (.NOT.JRACC) .AND. |
|
& (LWORK.LT.MAX0(2*M+N,6*N+2*N*N))) |
|
& .OR. (LSVEC .AND. RSVEC .AND. JRACC .AND. |
|
& LWORK.LT.MAX0(2*M+N,4*N+N*N,2*N+N*N+6))) |
& THEN |
& THEN |
INFO = - 17 |
INFO = - 17 |
ELSE |
ELSE |
Line 454
|
Line 475
|
IF ( INFO .NE. 0 ) THEN |
IF ( INFO .NE. 0 ) THEN |
* #:( |
* #:( |
CALL XERBLA( 'DGEJSV', - INFO ) |
CALL XERBLA( 'DGEJSV', - INFO ) |
|
RETURN |
END IF |
END IF |
* |
* |
* Quick return for void matrix (Y3K safe) |
* Quick return for void matrix (Y3K safe) |
Line 471
|
Line 493
|
* |
* |
*! NOTE: Make sure DLAMCH() does not fail on the target architecture. |
*! NOTE: Make sure DLAMCH() does not fail on the target architecture. |
* |
* |
|
|
EPSLN = DLAMCH('Epsilon') |
EPSLN = DLAMCH('Epsilon') |
SFMIN = DLAMCH('SafeMinimum') |
SFMIN = DLAMCH('SafeMinimum') |
SMALL = SFMIN / EPSLN |
SMALL = SFMIN / EPSLN |
Line 536
|
Line 557
|
END IF |
END IF |
IWORK(1) = 0 |
IWORK(1) = 0 |
IWORK(2) = 0 |
IWORK(2) = 0 |
|
IWORK(3) = 0 |
RETURN |
RETURN |
END IF |
END IF |
* |
* |
Line 834
|
Line 856
|
TEMP1 = DSQRT(SFMIN) |
TEMP1 = DSQRT(SFMIN) |
DO 3401 p = 2, N |
DO 3401 p = 2, N |
IF ( ( DABS(A(p,p)) .LT. (EPSLN*DABS(A(p-1,p-1))) ) .OR. |
IF ( ( DABS(A(p,p)) .LT. (EPSLN*DABS(A(p-1,p-1))) ) .OR. |
& ( DABS(A(p,p)) .LT. SMALL ) .OR. |
$ ( DABS(A(p,p)) .LT. SMALL ) .OR. |
& ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402 |
$ ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402 |
NR = NR + 1 |
NR = NR + 1 |
3401 CONTINUE |
3401 CONTINUE |
3402 CONTINUE |
3402 CONTINUE |
Line 851
|
Line 873
|
TEMP1 = DSQRT(SFMIN) |
TEMP1 = DSQRT(SFMIN) |
DO 3301 p = 2, N |
DO 3301 p = 2, N |
IF ( ( DABS(A(p,p)) .LT. SMALL ) .OR. |
IF ( ( DABS(A(p,p)) .LT. SMALL ) .OR. |
& ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302 |
$ ( L2KILL .AND. (DABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302 |
NR = NR + 1 |
NR = NR + 1 |
3301 CONTINUE |
3301 CONTINUE |
3302 CONTINUE |
3302 CONTINUE |
Line 883
|
Line 905
|
CALL DSCAL( p, ONE/TEMP1, V(1,p), 1 ) |
CALL DSCAL( p, ONE/TEMP1, V(1,p), 1 ) |
3053 CONTINUE |
3053 CONTINUE |
CALL DPOCON( 'U', N, V, LDV, ONE, TEMP1, |
CALL DPOCON( 'U', N, V, LDV, ONE, TEMP1, |
& WORK(N+1), IWORK(2*N+M+1), IERR ) |
$ WORK(N+1), IWORK(2*N+M+1), IERR ) |
ELSE IF ( LSVEC ) THEN |
ELSE IF ( LSVEC ) THEN |
* .. U is available as workspace |
* .. U is available as workspace |
CALL DLACPY( 'U', N, N, A, LDA, U, LDU ) |
CALL DLACPY( 'U', N, N, A, LDA, U, LDU ) |
Line 892
|
Line 914
|
CALL DSCAL( p, ONE/TEMP1, U(1,p), 1 ) |
CALL DSCAL( p, ONE/TEMP1, U(1,p), 1 ) |
3054 CONTINUE |
3054 CONTINUE |
CALL DPOCON( 'U', N, U, LDU, ONE, TEMP1, |
CALL DPOCON( 'U', N, U, LDU, ONE, TEMP1, |
& WORK(N+1), IWORK(2*N+M+1), IERR ) |
$ WORK(N+1), IWORK(2*N+M+1), IERR ) |
ELSE |
ELSE |
CALL DLACPY( 'U', N, N, A, LDA, WORK(N+1), N ) |
CALL DLACPY( 'U', N, N, A, LDA, WORK(N+1), N ) |
DO 3052 p = 1, N |
DO 3052 p = 1, N |
Line 901
|
Line 923
|
3052 CONTINUE |
3052 CONTINUE |
* .. the columns of R are scaled to have unit Euclidean lengths. |
* .. the columns of R are scaled to have unit Euclidean lengths. |
CALL DPOCON( 'U', N, WORK(N+1), N, ONE, TEMP1, |
CALL DPOCON( 'U', N, WORK(N+1), N, ONE, TEMP1, |
& WORK(N+N*N+1), IWORK(2*N+M+1), IERR ) |
$ WORK(N+N*N+1), IWORK(2*N+M+1), IERR ) |
END IF |
END IF |
SCONDA = ONE / DSQRT(TEMP1) |
SCONDA = ONE / DSQRT(TEMP1) |
* SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1). |
* SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1). |
Line 916
|
Line 938
|
* |
* |
* Phase 3: |
* Phase 3: |
* |
* |
|
|
IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN |
IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN |
* |
* |
* Singular Values only |
* Singular Values only |
Line 947
|
Line 968
|
TEMP1 = XSC*DABS(A(q,q)) |
TEMP1 = XSC*DABS(A(q,q)) |
DO 4949 p = 1, N |
DO 4949 p = 1, N |
IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) ) |
IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) ) |
& .OR. ( p .LT. q ) ) |
$ .OR. ( p .LT. q ) ) |
& A(p,q) = DSIGN( TEMP1, A(p,q) ) |
$ A(p,q) = DSIGN( TEMP1, A(p,q) ) |
4949 CONTINUE |
4949 CONTINUE |
4947 CONTINUE |
4947 CONTINUE |
ELSE |
ELSE |
Line 977
|
Line 998
|
TEMP1 = XSC*DABS(A(q,q)) |
TEMP1 = XSC*DABS(A(q,q)) |
DO 1949 p = 1, NR |
DO 1949 p = 1, NR |
IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) ) |
IF ( ( (p.GT.q) .AND. (DABS(A(p,q)).LE.TEMP1) ) |
& .OR. ( p .LT. q ) ) |
$ .OR. ( p .LT. q ) ) |
& A(p,q) = DSIGN( TEMP1, A(p,q) ) |
$ A(p,q) = DSIGN( TEMP1, A(p,q) ) |
1949 CONTINUE |
1949 CONTINUE |
1947 CONTINUE |
1947 CONTINUE |
ELSE |
ELSE |
Line 990
|
Line 1011
|
* the part which destroys triangular form (confusing?!)) |
* the part which destroys triangular form (confusing?!)) |
* |
* |
CALL DGESVJ( 'L', 'NoU', 'NoV', NR, NR, A, LDA, SVA, |
CALL DGESVJ( 'L', 'NoU', 'NoV', NR, NR, A, LDA, SVA, |
& N, V, LDV, WORK, LWORK, INFO ) |
$ N, V, LDV, WORK, LWORK, INFO ) |
* |
* |
SCALEM = WORK(1) |
SCALEM = WORK(1) |
NUMRANK = IDNINT(WORK(2)) |
NUMRANK = IDNINT(WORK(2)) |
Line 1009
|
Line 1030
|
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) |
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) |
* |
* |
CALL DGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA, |
CALL DGESVJ( 'L','U','N', N, NR, V,LDV, SVA, NR, A,LDA, |
& WORK, LWORK, INFO ) |
$ WORK, LWORK, INFO ) |
SCALEM = WORK(1) |
SCALEM = WORK(1) |
NUMRANK = IDNINT(WORK(2)) |
NUMRANK = IDNINT(WORK(2)) |
|
|
Line 1023
|
Line 1044
|
CALL DLACPY( 'Lower', NR, NR, A, LDA, V, LDV ) |
CALL DLACPY( 'Lower', NR, NR, A, LDA, V, LDV ) |
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) |
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) |
CALL DGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1), |
CALL DGEQRF( NR, NR, V, LDV, WORK(N+1), WORK(2*N+1), |
& LWORK-2*N, IERR ) |
$ LWORK-2*N, IERR ) |
DO 8998 p = 1, NR |
DO 8998 p = 1, NR |
CALL DCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 ) |
CALL DCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 ) |
8998 CONTINUE |
8998 CONTINUE |
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) |
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV ) |
* |
* |
CALL DGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U, |
CALL DGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U, |
& LDU, WORK(N+1), LWORK, INFO ) |
$ LDU, WORK(N+1), LWORK, INFO ) |
SCALEM = WORK(N+1) |
SCALEM = WORK(N+1) |
NUMRANK = IDNINT(WORK(N+2)) |
NUMRANK = IDNINT(WORK(N+2)) |
IF ( NR .LT. N ) THEN |
IF ( NR .LT. N ) THEN |
Line 1040
|
Line 1061
|
END IF |
END IF |
* |
* |
CALL DORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK, |
CALL DORMLQ( 'Left', 'Transpose', N, N, NR, A, LDA, WORK, |
& V, LDV, WORK(N+1), LWORK-N, IERR ) |
$ V, LDV, WORK(N+1), LWORK-N, IERR ) |
* |
* |
END IF |
END IF |
* |
* |
Line 1065
|
Line 1086
|
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU ) |
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU ) |
* |
* |
CALL DGEQRF( N, NR, U, LDU, WORK(N+1), WORK(2*N+1), |
CALL DGEQRF( N, NR, U, LDU, WORK(N+1), WORK(2*N+1), |
& LWORK-2*N, IERR ) |
$ LWORK-2*N, IERR ) |
* |
* |
DO 1967 p = 1, NR - 1 |
DO 1967 p = 1, NR - 1 |
CALL DCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 ) |
CALL DCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 ) |
Line 1073
|
Line 1094
|
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU ) |
CALL DLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, U(1,2), LDU ) |
* |
* |
CALL DGESVJ( 'Lower', 'U', 'N', NR,NR, U, LDU, SVA, NR, A, |
CALL DGESVJ( 'Lower', 'U', 'N', NR,NR, U, LDU, SVA, NR, A, |
& LDA, WORK(N+1), LWORK-N, INFO ) |
$ LDA, WORK(N+1), LWORK-N, INFO ) |
SCALEM = WORK(N+1) |
SCALEM = WORK(N+1) |
NUMRANK = IDNINT(WORK(N+2)) |
NUMRANK = IDNINT(WORK(N+2)) |
* |
* |
Line 1086
|
Line 1107
|
END IF |
END IF |
* |
* |
CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, |
CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, |
& LDU, WORK(N+1), LWORK-N, IERR ) |
$ LDU, WORK(N+1), LWORK-N, IERR ) |
* |
* |
IF ( ROWPIV ) |
IF ( ROWPIV ) |
& CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
$ CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
* |
* |
DO 1974 p = 1, N1 |
DO 1974 p = 1, N1 |
XSC = ONE / DNRM2( M, U(1,p), 1 ) |
XSC = ONE / DNRM2( M, U(1,p), 1 ) |
Line 1137
|
Line 1158
|
TEMP1 = XSC*DABS( V(q,q) ) |
TEMP1 = XSC*DABS( V(q,q) ) |
DO 2968 p = 1, N |
DO 2968 p = 1, N |
IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 ) |
IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 ) |
& .OR. ( p .LT. q ) ) |
$ .OR. ( p .LT. q ) ) |
& V(p,q) = DSIGN( TEMP1, V(p,q) ) |
$ V(p,q) = DSIGN( TEMP1, V(p,q) ) |
IF ( p. LT. q ) V(p,q) = - V(p,q) |
IF ( p .LT. q ) V(p,q) = - V(p,q) |
2968 CONTINUE |
2968 CONTINUE |
2969 CONTINUE |
2969 CONTINUE |
ELSE |
ELSE |
Line 1156
|
Line 1177
|
CALL DSCAL(NR-p+1,ONE/TEMP1,WORK(2*N+(p-1)*NR+p),1) |
CALL DSCAL(NR-p+1,ONE/TEMP1,WORK(2*N+(p-1)*NR+p),1) |
3950 CONTINUE |
3950 CONTINUE |
CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1, |
CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1, |
& WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR) |
$ WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR) |
CONDR1 = ONE / DSQRT(TEMP1) |
CONDR1 = ONE / DSQRT(TEMP1) |
* .. here need a second oppinion on the condition number |
* .. here need a second oppinion on the condition number |
* .. then assume worst case scenario |
* .. then assume worst case scenario |
Line 1172
|
Line 1193
|
* of a lower triangular matrix. |
* of a lower triangular matrix. |
* R1^t = Q2 * R2 |
* R1^t = Q2 * R2 |
CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), |
CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), |
& LWORK-2*N, IERR ) |
$ LWORK-2*N, IERR ) |
* |
* |
IF ( L2PERT ) THEN |
IF ( L2PERT ) THEN |
XSC = DSQRT(SMALL)/EPSLN |
XSC = DSQRT(SMALL)/EPSLN |
Line 1180
|
Line 1201
|
DO 3958 q = 1, p - 1 |
DO 3958 q = 1, p - 1 |
TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q))) |
TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q))) |
IF ( DABS(V(q,p)) .LE. TEMP1 ) |
IF ( DABS(V(q,p)) .LE. TEMP1 ) |
& V(q,p) = DSIGN( TEMP1, V(q,p) ) |
$ V(q,p) = DSIGN( TEMP1, V(q,p) ) |
3958 CONTINUE |
3958 CONTINUE |
3959 CONTINUE |
3959 CONTINUE |
END IF |
END IF |
* |
* |
IF ( NR .NE. N ) |
IF ( NR .NE. N ) |
|
$ CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N ) |
* .. save ... |
* .. save ... |
& CALL DLACPY( 'A', N, NR, V, LDV, WORK(2*N+1), N ) |
|
* |
* |
* .. this transposed copy should be better than naive |
* .. this transposed copy should be better than naive |
DO 1969 p = 1, NR - 1 |
DO 1969 p = 1, NR - 1 |
Line 1210
|
Line 1231
|
IWORK(N+p) = 0 |
IWORK(N+p) = 0 |
3003 CONTINUE |
3003 CONTINUE |
CALL DGEQP3( N, NR, V, LDV, IWORK(N+1), WORK(N+1), |
CALL DGEQP3( N, NR, V, LDV, IWORK(N+1), WORK(N+1), |
& WORK(2*N+1), LWORK-2*N, IERR ) |
$ WORK(2*N+1), LWORK-2*N, IERR ) |
** CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), |
** CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), |
** & LWORK-2*N, IERR ) |
** $ LWORK-2*N, IERR ) |
IF ( L2PERT ) THEN |
IF ( L2PERT ) THEN |
XSC = DSQRT(SMALL) |
XSC = DSQRT(SMALL) |
DO 3969 p = 2, NR |
DO 3969 p = 2, NR |
DO 3968 q = 1, p - 1 |
DO 3968 q = 1, p - 1 |
TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q))) |
TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q))) |
IF ( DABS(V(q,p)) .LE. TEMP1 ) |
IF ( DABS(V(q,p)) .LE. TEMP1 ) |
& V(q,p) = DSIGN( TEMP1, V(q,p) ) |
$ V(q,p) = DSIGN( TEMP1, V(q,p) ) |
3968 CONTINUE |
3968 CONTINUE |
3969 CONTINUE |
3969 CONTINUE |
END IF |
END IF |
Line 1239
|
Line 1260
|
END IF |
END IF |
* Now, compute R2 = L3 * Q3, the LQ factorization. |
* Now, compute R2 = L3 * Q3, the LQ factorization. |
CALL DGELQF( NR, NR, V, LDV, WORK(2*N+N*NR+1), |
CALL DGELQF( NR, NR, V, LDV, WORK(2*N+N*NR+1), |
& WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR ) |
$ WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR ) |
* .. and estimate the condition number |
* .. and estimate the condition number |
CALL DLACPY( 'L',NR,NR,V,LDV,WORK(2*N+N*NR+NR+1),NR ) |
CALL DLACPY( 'L',NR,NR,V,LDV,WORK(2*N+N*NR+NR+1),NR ) |
DO 4950 p = 1, NR |
DO 4950 p = 1, NR |
Line 1247
|
Line 1268
|
CALL DSCAL( p, ONE/TEMP1, WORK(2*N+N*NR+NR+p), NR ) |
CALL DSCAL( p, ONE/TEMP1, WORK(2*N+N*NR+NR+p), NR ) |
4950 CONTINUE |
4950 CONTINUE |
CALL DPOCON( 'L',NR,WORK(2*N+N*NR+NR+1),NR,ONE,TEMP1, |
CALL DPOCON( 'L',NR,WORK(2*N+N*NR+NR+1),NR,ONE,TEMP1, |
& WORK(2*N+N*NR+NR+NR*NR+1),IWORK(M+2*N+1),IERR ) |
$ WORK(2*N+N*NR+NR+NR*NR+1),IWORK(M+2*N+1),IERR ) |
CONDR2 = ONE / DSQRT(TEMP1) |
CONDR2 = ONE / DSQRT(TEMP1) |
* |
* |
IF ( CONDR2 .GE. COND_OK ) THEN |
IF ( CONDR2 .GE. COND_OK ) THEN |
Line 1284
|
Line 1305
|
IF ( CONDR1 .LT. COND_OK ) THEN |
IF ( CONDR1 .LT. COND_OK ) THEN |
* |
* |
CALL DGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U, |
CALL DGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U, |
& LDU,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,INFO ) |
$ LDU,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,INFO ) |
SCALEM = WORK(2*N+N*NR+NR+1) |
SCALEM = WORK(2*N+N*NR+NR+1) |
NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) |
NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) |
DO 3970 p = 1, NR |
DO 3970 p = 1, NR |
Line 1294
|
Line 1315
|
|
|
* .. pick the right matrix equation and solve it |
* .. pick the right matrix equation and solve it |
* |
* |
IF ( NR. EQ. N ) THEN |
IF ( NR .EQ. N ) THEN |
* :)) .. best case, R1 is inverted. The solution of this matrix |
* :)) .. best case, R1 is inverted. The solution of this matrix |
* equation is Q2*V2 = the product of the Jacobi rotations |
* equation is Q2*V2 = the product of the Jacobi rotations |
* used in DGESVJ, premultiplied with the orthogonal matrix |
* used in DGESVJ, premultiplied with the orthogonal matrix |
Line 1306
|
Line 1327
|
* used in DGESVJ. The Q-factor from the second QR |
* used in DGESVJ. The Q-factor from the second QR |
* factorization is then built in explicitly. |
* factorization is then built in explicitly. |
CALL DTRSM('L','U','T','N',NR,NR,ONE,WORK(2*N+1), |
CALL DTRSM('L','U','T','N',NR,NR,ONE,WORK(2*N+1), |
& N,V,LDV) |
$ N,V,LDV) |
IF ( NR .LT. N ) THEN |
IF ( NR .LT. N ) THEN |
CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV) |
CALL DLASET('A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV) |
CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV) |
CALL DLASET('A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV) |
CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV) |
CALL DLASET('A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV) |
END IF |
END IF |
CALL DORMQR('L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), |
CALL DORMQR('L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), |
& V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR) |
$ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR) |
END IF |
END IF |
* |
* |
ELSE IF ( CONDR2 .LT. COND_OK ) THEN |
ELSE IF ( CONDR2 .LT. COND_OK ) THEN |
Line 1325
|
Line 1346
|
* the lower triangular L3 from the LQ factorization of |
* the lower triangular L3 from the LQ factorization of |
* R2=L3*Q3), pre-multiplied with the transposed Q3. |
* R2=L3*Q3), pre-multiplied with the transposed Q3. |
CALL DGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U, |
CALL DGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U, |
& LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO ) |
$ LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO ) |
SCALEM = WORK(2*N+N*NR+NR+1) |
SCALEM = WORK(2*N+N*NR+NR+1) |
NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) |
NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) |
DO 3870 p = 1, NR |
DO 3870 p = 1, NR |
Line 1348
|
Line 1369
|
CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV ) |
CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV ) |
END IF |
END IF |
CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), |
CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), |
& V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) |
$ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) |
ELSE |
ELSE |
* Last line of defense. |
* Last line of defense. |
* #:( This is a rather pathological case: no scaled condition |
* #:( This is a rather pathological case: no scaled condition |
Line 1362
|
Line 1383
|
* Compute the full SVD of L3 using DGESVJ with explicit |
* Compute the full SVD of L3 using DGESVJ with explicit |
* accumulation of Jacobi rotations. |
* accumulation of Jacobi rotations. |
CALL DGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U, |
CALL DGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U, |
& LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO ) |
$ LDU, WORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, INFO ) |
SCALEM = WORK(2*N+N*NR+NR+1) |
SCALEM = WORK(2*N+N*NR+NR+1) |
NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) |
NUMRANK = IDNINT(WORK(2*N+N*NR+NR+2)) |
IF ( NR .LT. N ) THEN |
IF ( NR .LT. N ) THEN |
Line 1371
|
Line 1392
|
CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV ) |
CALL DLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV ) |
END IF |
END IF |
CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), |
CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), |
& V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) |
$ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) |
* |
* |
CALL DORMLQ( 'L', 'T', NR, NR, NR, WORK(2*N+1), N, |
CALL DORMLQ( 'L', 'T', NR, NR, NR, WORK(2*N+1), N, |
& WORK(2*N+N*NR+1), U, LDU, WORK(2*N+N*NR+NR+1), |
$ WORK(2*N+N*NR+1), U, LDU, WORK(2*N+N*NR+NR+1), |
& LWORK-2*N-N*NR-NR, IERR ) |
$ LWORK-2*N-N*NR-NR, IERR ) |
DO 773 q = 1, NR |
DO 773 q = 1, NR |
DO 772 p = 1, NR |
DO 772 p = 1, NR |
WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q) |
WORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q) |
Line 1401
|
Line 1422
|
973 CONTINUE |
973 CONTINUE |
XSC = ONE / DNRM2( N, V(1,q), 1 ) |
XSC = ONE / DNRM2( N, V(1,q), 1 ) |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
& CALL DSCAL( N, XSC, V(1,q), 1 ) |
$ CALL DSCAL( N, XSC, V(1,q), 1 ) |
1972 CONTINUE |
1972 CONTINUE |
* At this moment, V contains the right singular vectors of A. |
* At this moment, V contains the right singular vectors of A. |
* Next, assemble the left singular vector matrix U (M x N). |
* Next, assemble the left singular vector matrix U (M x N). |
Line 1417
|
Line 1438
|
* matrix U. This applies to all cases. |
* matrix U. This applies to all cases. |
* |
* |
CALL DORMQR( 'Left', 'No_Tr', M, N1, N, A, LDA, WORK, U, |
CALL DORMQR( 'Left', 'No_Tr', M, N1, N, A, LDA, WORK, U, |
& LDU, WORK(N+1), LWORK-N, IERR ) |
$ LDU, WORK(N+1), LWORK-N, IERR ) |
|
|
* The columns of U are normalized. The cost is O(M*N) flops. |
* The columns of U are normalized. The cost is O(M*N) flops. |
TEMP1 = DSQRT(DBLE(M)) * EPSLN |
TEMP1 = DSQRT(DBLE(M)) * EPSLN |
DO 1973 p = 1, NR |
DO 1973 p = 1, NR |
XSC = ONE / DNRM2( M, U(1,p), 1 ) |
XSC = ONE / DNRM2( M, U(1,p), 1 ) |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
& CALL DSCAL( M, XSC, U(1,p), 1 ) |
$ CALL DSCAL( M, XSC, U(1,p), 1 ) |
1973 CONTINUE |
1973 CONTINUE |
* |
* |
* If the initial QRF is computed with row pivoting, the left |
* If the initial QRF is computed with row pivoting, the left |
* singular vectors must be adjusted. |
* singular vectors must be adjusted. |
* |
* |
IF ( ROWPIV ) |
IF ( ROWPIV ) |
& CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
$ CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
* |
* |
ELSE |
ELSE |
* |
* |
Line 1452
|
Line 1473
|
END IF |
END IF |
* |
* |
CALL DGESVJ( 'Upper', 'U', 'N', N, N, WORK(N+1), N, SVA, |
CALL DGESVJ( 'Upper', 'U', 'N', N, N, WORK(N+1), N, SVA, |
& N, U, LDU, WORK(N+N*N+1), LWORK-N-N*N, INFO ) |
$ N, U, LDU, WORK(N+N*N+1), LWORK-N-N*N, INFO ) |
* |
* |
SCALEM = WORK(N+N*N+1) |
SCALEM = WORK(N+N*N+1) |
NUMRANK = IDNINT(WORK(N+N*N+2)) |
NUMRANK = IDNINT(WORK(N+N*N+2)) |
Line 1462
|
Line 1483
|
6970 CONTINUE |
6970 CONTINUE |
* |
* |
CALL DTRSM( 'Left', 'Upper', 'NoTrans', 'No UD', N, N, |
CALL DTRSM( 'Left', 'Upper', 'NoTrans', 'No UD', N, N, |
& ONE, A, LDA, WORK(N+1), N ) |
$ ONE, A, LDA, WORK(N+1), N ) |
DO 6972 p = 1, N |
DO 6972 p = 1, N |
CALL DCOPY( N, WORK(N+p), N, V(IWORK(p),1), LDV ) |
CALL DCOPY( N, WORK(N+p), N, V(IWORK(p),1), LDV ) |
6972 CONTINUE |
6972 CONTINUE |
Line 1470
|
Line 1491
|
DO 6971 p = 1, N |
DO 6971 p = 1, N |
XSC = ONE / DNRM2( N, V(1,p), 1 ) |
XSC = ONE / DNRM2( N, V(1,p), 1 ) |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
& CALL DSCAL( N, XSC, V(1,p), 1 ) |
$ CALL DSCAL( N, XSC, V(1,p), 1 ) |
6971 CONTINUE |
6971 CONTINUE |
* |
* |
* Assemble the left singular vector matrix U (M x N). |
* Assemble the left singular vector matrix U (M x N). |
Line 1483
|
Line 1504
|
END IF |
END IF |
END IF |
END IF |
CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, |
CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, |
& LDU, WORK(N+1), LWORK-N, IERR ) |
$ LDU, WORK(N+1), LWORK-N, IERR ) |
TEMP1 = DSQRT(DBLE(M))*EPSLN |
TEMP1 = DSQRT(DBLE(M))*EPSLN |
DO 6973 p = 1, N1 |
DO 6973 p = 1, N1 |
XSC = ONE / DNRM2( M, U(1,p), 1 ) |
XSC = ONE / DNRM2( M, U(1,p), 1 ) |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
& CALL DSCAL( M, XSC, U(1,p), 1 ) |
$ CALL DSCAL( M, XSC, U(1,p), 1 ) |
6973 CONTINUE |
6973 CONTINUE |
* |
* |
IF ( ROWPIV ) |
IF ( ROWPIV ) |
& CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
$ CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
* |
* |
END IF |
END IF |
* |
* |
Line 1520
|
Line 1541
|
TEMP1 = XSC*DABS( V(q,q) ) |
TEMP1 = XSC*DABS( V(q,q) ) |
DO 5968 p = 1, N |
DO 5968 p = 1, N |
IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 ) |
IF ( ( p .GT. q ) .AND. ( DABS(V(p,q)) .LE. TEMP1 ) |
& .OR. ( p .LT. q ) ) |
$ .OR. ( p .LT. q ) ) |
& V(p,q) = DSIGN( TEMP1, V(p,q) ) |
$ V(p,q) = DSIGN( TEMP1, V(p,q) ) |
IF ( p. LT. q ) V(p,q) = - V(p,q) |
IF ( p .LT. q ) V(p,q) = - V(p,q) |
5968 CONTINUE |
5968 CONTINUE |
5969 CONTINUE |
5969 CONTINUE |
ELSE |
ELSE |
Line 1530
|
Line 1551
|
END IF |
END IF |
|
|
CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), |
CALL DGEQRF( N, NR, V, LDV, WORK(N+1), WORK(2*N+1), |
& LWORK-2*N, IERR ) |
$ LWORK-2*N, IERR ) |
CALL DLACPY( 'L', N, NR, V, LDV, WORK(2*N+1), N ) |
CALL DLACPY( 'L', N, NR, V, LDV, WORK(2*N+1), N ) |
* |
* |
DO 7969 p = 1, NR |
DO 7969 p = 1, NR |
Line 1550
|
Line 1571
|
END IF |
END IF |
|
|
CALL DGESVJ( 'G', 'U', 'V', NR, NR, U, LDU, SVA, |
CALL DGESVJ( 'G', 'U', 'V', NR, NR, U, LDU, SVA, |
& N, V, LDV, WORK(2*N+N*NR+1), LWORK-2*N-N*NR, INFO ) |
$ N, V, LDV, WORK(2*N+N*NR+1), LWORK-2*N-N*NR, INFO ) |
SCALEM = WORK(2*N+N*NR+1) |
SCALEM = WORK(2*N+N*NR+1) |
NUMRANK = IDNINT(WORK(2*N+N*NR+2)) |
NUMRANK = IDNINT(WORK(2*N+N*NR+2)) |
|
|
Line 1561
|
Line 1582
|
END IF |
END IF |
|
|
CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), |
CALL DORMQR( 'L','N',N,N,NR,WORK(2*N+1),N,WORK(N+1), |
& V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) |
$ V,LDV,WORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) |
* |
* |
* Permute the rows of V using the (column) permutation from the |
* Permute the rows of V using the (column) permutation from the |
* first QRF. Also, scale the columns to make them unit in |
* first QRF. Also, scale the columns to make them unit in |
Line 1577
|
Line 1598
|
8973 CONTINUE |
8973 CONTINUE |
XSC = ONE / DNRM2( N, V(1,q), 1 ) |
XSC = ONE / DNRM2( N, V(1,q), 1 ) |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) |
& CALL DSCAL( N, XSC, V(1,q), 1 ) |
$ CALL DSCAL( N, XSC, V(1,q), 1 ) |
7972 CONTINUE |
7972 CONTINUE |
* |
* |
* At this moment, V contains the right singular vectors of A. |
* At this moment, V contains the right singular vectors of A. |
Line 1592
|
Line 1613
|
END IF |
END IF |
* |
* |
CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, |
CALL DORMQR( 'Left', 'No Tr', M, N1, N, A, LDA, WORK, U, |
& LDU, WORK(N+1), LWORK-N, IERR ) |
$ LDU, WORK(N+1), LWORK-N, IERR ) |
* |
* |
IF ( ROWPIV ) |
IF ( ROWPIV ) |
& CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
$ CALL DLASWP( N1, U, LDU, 1, M-1, IWORK(2*N+1), -1 ) |
* |
* |
* |
* |
END IF |
END IF |