version 1.15, 2016/08/27 15:34:21
|
version 1.18, 2018/05/29 06:55:16
|
Line 2
|
Line 2
|
* |
* |
* =========== DOCUMENTATION =========== |
* =========== DOCUMENTATION =========== |
* |
* |
* Online html documentation available at |
* Online html documentation available at |
* http://www.netlib.org/lapack/explore-html/ |
* http://www.netlib.org/lapack/explore-html/ |
* |
* |
*> \htmlonly |
*> \htmlonly |
*> Download DGEJSV + dependencies |
*> Download DGEJSV + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgejsv.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgejsv.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgejsv.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgejsv.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgejsv.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgejsv.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
Line 21
|
Line 21
|
* 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 ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* IMPLICIT NONE |
* IMPLICIT NONE |
* INTEGER INFO, LDA, LDU, LDV, LWORK, M, N |
* INTEGER INFO, LDA, LDU, LDV, LWORK, M, N |
Line 32
|
Line 32
|
* INTEGER IWORK( * ) |
* INTEGER IWORK( * ) |
* CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV |
* CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
Line 87
|
Line 87
|
*> rows, then using this condition number gives too pessimistic |
*> rows, then using this condition number gives too pessimistic |
*> error bound. |
*> error bound. |
*> = 'A': Small singular values are the noise and the matrix is treated |
*> = 'A': Small singular values are the noise and the matrix is treated |
*> as numerically rank defficient. The error in the computed |
*> as numerically rank deficient. The error in the computed |
*> singular values is bounded by f(m,n)*epsilon*||A||. |
*> singular values is bounded by f(m,n)*epsilon*||A||. |
*> The computed SVD A = U * S * V^t restores A up to |
*> The computed SVD A = U * S * V^t restores A up to |
*> f(m,n)*epsilon*||A||. |
*> f(m,n)*epsilon*||A||. |
Line 271
|
Line 271
|
*> |
*> |
*> \param[out] WORK |
*> \param[out] WORK |
*> \verbatim |
*> \verbatim |
*> WORK is DOUBLE PRECISION array, dimension at least LWORK. |
*> WORK is DOUBLE PRECISION array, dimension (LWORK) |
*> On exit, if N.GT.0 .AND. M.GT.0 (else not referenced), |
*> 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 319
|
Line 319
|
*> ->> 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 DGEQP3 and DGEQRF. |
*> block size for DGEQP3 and DGEQRF. |
*> In general, optimal LWORK is computed as |
*> In general, optimal LWORK is computed as |
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), 7). |
*> 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+4*N,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 |
*> ->> For optimal performance (blocked code) the optimal value |
*> is LWORK >= max(2*M+N,3*N+(N+1)*NB, N*N+4*N, 7). |
*> 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 |
*> In general, the optimal length LWORK is computed as |
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), |
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF), |
*> N+N*N+LWORK(DPOCON),7). |
*> 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'), |
Line 335
|
Line 335
|
*> -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7), |
*> -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7), |
*> where NB is the optimal block size for DGEQP3, DGEQRF, DGELQF, |
*> where NB is the optimal block size for DGEQP3, DGEQRF, DGELQF, |
*> DORMLQ. In general, the optimal length LWORK is computed as |
*> DORMLQ. In general, the optimal length LWORK is computed as |
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON), |
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3), N+LWORK(DPOCON), |
*> N+LWORK(DGELQF), 2*N+LWORK(DGEQRF), N+LWORK(DORMLQ)). |
*> N+LWORK(DGELQF), 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 |
Line 346
|
Line 346
|
*> where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR. |
*> where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR. |
*> In general, the optimal length LWORK is computed as |
*> In general, the optimal length LWORK is computed as |
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON), |
*> LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON), |
*> 2*N+LWORK(DGEQRF), N+LWORK(DORMQR)). |
*> 2*N+LWORK(DGEQRF), N+LWORK(DORMQR)). |
*> Here LWORK(DORMQR) equals N*NB (for JOBU.EQ.'U') or |
*> Here LWORK(DORMQR) equals N*NB (for JOBU.EQ.'U') or |
*> M*NB (for JOBU.EQ.'F'). |
*> M*NB (for JOBU.EQ.'F'). |
*> |
*> |
*> If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and |
*> If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and |
*> -> if JOBV.EQ.'V' |
*> -> if JOBV.EQ.'V' |
*> the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N). |
*> the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N). |
*> -> if JOBV.EQ.'J' the minimal requirement is |
*> -> if JOBV.EQ.'J' the minimal requirement is |
*> LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6). |
*> LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6). |
*> -> For optimal performance, LWORK should be additionally |
*> -> For optimal performance, LWORK should be additionally |
*> larger than N+M*NB, where NB is the optimal block size |
*> larger than N+M*NB, where NB is the optimal block size |
Line 362
|
Line 362
|
*> |
*> |
*> \param[out] IWORK |
*> \param[out] IWORK |
*> \verbatim |
*> \verbatim |
*> IWORK is INTEGER array, dimension M+3*N. |
*> IWORK is INTEGER array, dimension (M+3*N). |
*> On exit, |
*> On exit, |
*> IWORK(1) = the numerical rank determined after the initial |
*> IWORK(1) = the numerical rank determined after the initial |
*> QR factorization with pivoting. See the descriptions |
*> QR factorization with pivoting. See the descriptions |
Line 378
|
Line 378
|
*> \verbatim |
*> \verbatim |
*> INFO is INTEGER |
*> INFO is INTEGER |
*> < 0 : if INFO = -i, then the i-th argument had an illegal value. |
*> < 0 : if INFO = -i, then the i-th argument had an illegal value. |
*> = 0 : successfull exit; |
*> = 0 : successful exit; |
*> > 0 : DGEJSV did not converge in the maximal allowed number |
*> > 0 : DGEJSV did not converge in the maximal allowed number |
*> of sweeps. The computed values may be inaccurate. |
*> of sweeps. The computed values may be inaccurate. |
*> \endverbatim |
*> \endverbatim |
Line 386
|
Line 386
|
* Authors: |
* Authors: |
* ======== |
* ======== |
* |
* |
*> \author Univ. of Tennessee |
*> \author Univ. of Tennessee |
*> \author Univ. of California Berkeley |
*> \author Univ. of California Berkeley |
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
* |
*> \date June 2016 |
*> \date June 2016 |
* |
* |
Line 428
|
Line 428
|
*> The rank revealing QR factorization (in this code: DGEQP3) should be |
*> The rank revealing QR factorization (in this code: DGEQP3) should be |
*> implemented as in [3]. We have a new version of DGEQP3 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 deficient 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 initial QRF with |
*> column pivoting can be preprocessed by the QRF without pivoting. That |
*> column pivoting can be preprocessed by the QRF without pivoting. That |
*> well known trick is not used in DGEJSV because in some cases heavy row |
*> well known trick is not used in DGEJSV because in some cases heavy row |
*> weighting can be treated with complete pivoting. The overhead in cases |
*> weighting can be treated with complete pivoting. The overhead in cases |
Line 476
|
Line 476
|
$ 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 computational routine (version 3.6.1) -- |
* -- LAPACK computational routine (version 3.7.1) -- |
* -- 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..-- |
* June 2016 |
* June 2016 |
Line 562
|
Line 562
|
ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN |
ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN |
INFO = - 13 |
INFO = - 13 |
ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN |
ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN |
INFO = - 14 |
INFO = - 15 |
ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND. |
ELSE IF ( (.NOT.(LSVEC .OR. RSVEC .OR. ERREST).AND. |
& (LWORK .LT. MAX(7,4*N+1,2*M+N))) .OR. |
& (LWORK .LT. MAX(7,4*N+1,2*M+N))) .OR. |
& (.NOT.(LSVEC .OR. RSVEC) .AND. ERREST .AND. |
& (.NOT.(LSVEC .OR. RSVEC) .AND. ERREST .AND. |
Line 571
|
Line 571
|
& .OR. |
& .OR. |
& (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX(7,2*M+N,4*N+1))) |
& (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX(7,2*M+N,4*N+1))) |
& .OR. |
& .OR. |
& (LSVEC .AND. RSVEC .AND. (.NOT.JRACC) .AND. |
& (LSVEC .AND. RSVEC .AND. (.NOT.JRACC) .AND. |
& (LWORK.LT.MAX(2*M+N,6*N+2*N*N))) |
& (LWORK.LT.MAX(2*M+N,6*N+2*N*N))) |
& .OR. (LSVEC .AND. RSVEC .AND. JRACC .AND. |
& .OR. (LSVEC .AND. RSVEC .AND. JRACC .AND. |
& LWORK.LT.MAX(2*M+N,4*N+N*N,2*N+N*N+6))) |
& LWORK.LT.MAX(2*M+N,4*N+N*N,2*N+N*N+6))) |
Line 967
|
Line 967
|
ELSE IF ( L2RANK ) THEN |
ELSE IF ( L2RANK ) THEN |
* .. similarly as above, only slightly more gentle (less agressive). |
* .. similarly as above, only slightly more gentle (less agressive). |
* Sudden drop on the diagonal of R1 is used as the criterion for |
* Sudden drop on the diagonal of R1 is used as the criterion for |
* close-to-rank-defficient. |
* close-to-rank-deficient. |
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. |