version 1.7, 2011/11/21 20:42:50
|
version 1.21, 2023/08/07 08:38:48
|
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 51
|
Line 51
|
*> the right singular vectors of [A], respectively. The matrices [U] and [V] |
*> the right singular vectors of [A], respectively. The matrices [U] and [V] |
*> are computed and stored in the arrays U and V, respectively. The diagonal |
*> are computed and stored in the arrays U and V, respectively. The diagonal |
*> of [SIGMA] is computed and stored in the array SVA. |
*> of [SIGMA] is computed and stored in the array SVA. |
|
*> DGEJSV can sometimes compute tiny singular values and their singular vectors much |
|
*> more accurately than other SVD routines, see below under Further Details. |
*> \endverbatim |
*> \endverbatim |
* |
* |
* Arguments: |
* Arguments: |
Line 80
|
Line 82
|
*> desirable, then this option is advisable. The input matrix A |
*> desirable, then this option is advisable. The input matrix A |
*> is preprocessed with QR factorization with FULL (row and |
*> is preprocessed with QR factorization with FULL (row and |
*> column) pivoting. |
*> column) pivoting. |
*> = 'G' Computation as with 'F' with an additional estimate of the |
*> = 'G': Computation as with 'F' with an additional estimate of the |
*> condition number of B, where A=D*B. If A has heavily weighted |
*> condition number of B, where A=D*B. If A has heavily weighted |
*> 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 131
|
Line 133
|
*> specified range. If A .NE. 0 is scaled so that the largest singular |
*> specified range. If A .NE. 0 is scaled so that the largest singular |
*> value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues |
*> value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues |
*> the licence to kill columns of A whose norm in c*A is less than |
*> the licence to kill columns of A whose norm in c*A is less than |
*> DSQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN, |
*> DSQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN, |
*> where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E'). |
*> where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E'). |
*> = 'N': Do not kill small columns of c*A. This option assumes that |
*> = 'N': Do not kill small columns of c*A. This option assumes that |
*> BLAS and QR factorizations and triangular solvers are |
*> BLAS and QR factorizations and triangular solvers are |
Line 228
|
Line 230
|
*> If JOBU = 'F', then U contains on exit the M-by-M matrix of |
*> If JOBU = 'F', then U contains on exit the M-by-M matrix of |
*> the left singular vectors, including an ONB |
*> the left singular vectors, including an ONB |
*> of the orthogonal complement of the Range(A). |
*> of the orthogonal complement of the Range(A). |
*> If JOBU = 'W' .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N), |
*> If JOBU = 'W' .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N), |
*> then U is used as workspace if the procedure |
*> then U is used as workspace if the procedure |
*> replaces A with A^t. In that case, [V] is computed |
*> replaces A with A^t. In that case, [V] is computed |
*> in U as left singular vectors of A^t and then |
*> in U as left singular vectors of A^t and then |
*> copied back to the V array. This 'W' option is just |
*> copied back to the V array. This 'W' option is just |
*> a reminder to the caller that in this case U is |
*> a reminder to the caller that in this case U is |
*> reserved as workspace of length N*N. |
*> reserved as workspace of length N*N. |
*> If JOBU = 'N' U is not referenced. |
*> If JOBU = 'N' U is not referenced, unless JOBT='T'. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDU |
*> \param[in] LDU |
Line 250
|
Line 252
|
*> V is DOUBLE PRECISION array, dimension ( LDV, N ) |
*> V is DOUBLE PRECISION array, dimension ( LDV, N ) |
*> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of |
*> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of |
*> the right singular vectors; |
*> the right singular vectors; |
*> If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N), |
*> If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N), |
*> then V is used as workspace if the pprocedure |
*> then V is used as workspace if the pprocedure |
*> replaces A with A^t. In that case, [U] is computed |
*> replaces A with A^t. In that case, [U] is computed |
*> in V as right singular vectors of A^t and then |
*> in V as right singular vectors of A^t and then |
*> copied back to the U array. This 'W' option is just |
*> copied back to the U array. This 'W' option is just |
*> a reminder to the caller that in this case V is |
*> a reminder to the caller that in this case V is |
*> reserved as workspace of length N*N. |
*> reserved as workspace of length N*N. |
*> If JOBV = 'N' V is not referenced. |
*> If JOBV = 'N' V is not referenced, unless JOBT='T'. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDV |
*> \param[in] LDV |
Line 269
|
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 > 0 .AND. M > 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().) |
*> WORK(2) = See the description of WORK(1). |
*> WORK(2) = See the description of WORK(1). |
*> WORK(3) = SCONDA is an estimate for the condition number of |
*> WORK(3) = SCONDA is an estimate for the condition number of |
*> column equilibrated A. (If JOBA .EQ. 'E' or 'G') |
*> column equilibrated A. (If JOBA = 'E' or 'G') |
*> SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1). |
*> SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1). |
*> It is computed using DPOCON. It holds |
*> It is computed using DPOCON. It holds |
*> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA |
*> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA |
Line 287
|
Line 289
|
*> singular values might be lost. |
*> singular values might be lost. |
*> |
*> |
*> If full SVD is needed, the following two condition numbers are |
*> If full SVD is needed, the following two condition numbers are |
*> useful for the analysis of the algorithm. They are provied for |
*> useful for the analysis of the algorithm. They are provided for |
*> a developer/implementer who is familiar with the details of |
*> a developer/implementer who is familiar with the details of |
*> the method. |
*> the method. |
*> |
*> |
Line 295
|
Line 297
|
*> triangular factor in the first QR factorization. |
*> triangular factor in the first QR factorization. |
*> WORK(5) = an estimate of the scaled condition number of the |
*> WORK(5) = an estimate of the scaled condition number of the |
*> triangular factor in the second QR factorization. |
*> triangular factor in the second QR factorization. |
*> The following two parameters are computed if JOBT .EQ. 'T'. |
*> The following two parameters are computed if JOBT = 'T'. |
*> They are provided for a developer/implementer who is familiar |
*> They are provided for a developer/implementer who is familiar |
*> with the details of the method. |
*> with the details of the method. |
*> |
*> |
Line 311
|
Line 313
|
*> Length of WORK to confirm proper allocation of work space. |
*> Length of WORK to confirm proper allocation of work space. |
*> 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 = 'N', JOBV = 'N') and |
*> -> .. no scaled condition estimate required (JOBE.EQ.'N'): |
*> -> .. no scaled condition estimate required (JOBE = '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 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 = 'V'), |
*> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7). |
*> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7). |
*> -> 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, DGELQ, |
*> 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(DGELQ), 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 |
*> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7). |
*> -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7). |
*> -> For optimal performance: |
*> -> For optimal performance: |
*> if JOBU.EQ.'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7), |
*> if JOBU = '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 JOBU = 'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7), |
*> 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 = 'U') or |
*> M*NB (for JOBU.EQ.'F'). |
*> M*NB (for JOBU = 'F'). |
*> |
*> |
*> If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and |
*> If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and |
*> -> if JOBV.EQ.'V' |
*> -> if JOBV = '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 = '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 360
|
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 |
*> of JOBA and JOBR. |
*> of JOBA and JOBR. |
*> IWORK(2) = the number of the computed nonzero singular values |
*> IWORK(2) = the number of the computed nonzero singular values |
*> IWORK(3) = if nonzero, a warning message: |
*> IWORK(3) = if nonzero, a warning message: |
*> If IWORK(3).EQ.1 then some of the column norms of A |
*> If IWORK(3) = 1 then some of the column norms of A |
*> were denormalized floats. The requested high accuracy |
*> were denormalized floats. The requested high accuracy |
*> is not warranted by the data. |
*> is not warranted by the data. |
*> \endverbatim |
*> \endverbatim |
Line 375
|
Line 377
|
*> \param[out] INFO |
*> \param[out] INFO |
*> \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 |
* |
* |
* 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 November 2011 |
*> \ingroup doubleGEsing |
* |
|
*> \ingroup doubleGEcomputational |
|
* |
* |
*> \par Further Details: |
*> \par Further Details: |
* ===================== |
* ===================== |
Line 426
|
Line 426
|
*> 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 474
|
Line 474
|
$ 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.4.0) -- |
* -- LAPACK computational routine -- |
* -- 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..-- |
* November 2011 |
|
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
IMPLICIT NONE |
IMPLICIT NONE |
Line 506
|
Line 505
|
$ NOSCAL, ROWPIV, RSVEC, TRANSP |
$ NOSCAL, ROWPIV, RSVEC, TRANSP |
* .. |
* .. |
* .. Intrinsic Functions .. |
* .. Intrinsic Functions .. |
INTRINSIC DABS, DLOG, DMAX1, DMIN1, DBLE, |
INTRINSIC DABS, DLOG, MAX, MIN, DBLE, IDNINT, DSIGN, DSQRT |
$ MAX0, MIN0, IDNINT, DSIGN, DSQRT |
|
* .. |
* .. |
* .. External Functions .. |
* .. External Functions .. |
DOUBLE PRECISION DLAMCH, DNRM2 |
DOUBLE PRECISION DLAMCH, DNRM2 |
Line 561
|
Line 559
|
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. MAX0(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. |
& (LWORK .LT. MAX0(7,4*N+N*N,2*M+N))) .OR. |
& (LWORK .LT. MAX(7,4*N+N*N,2*M+N))) .OR. |
& (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX0(7,2*M+N,4*N+1))) |
& (LSVEC .AND. (.NOT.RSVEC) .AND. (LWORK .LT. MAX(7,2*M+N,4*N+1))) |
& .OR. |
& .OR. |
& (RSVEC .AND. (.NOT.LSVEC) .AND. (LWORK .LT. MAX0(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.MAX0(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.MAX0(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))) |
& THEN |
& THEN |
INFO = - 17 |
INFO = - 17 |
ELSE |
ELSE |
Line 589
|
Line 587
|
* |
* |
* Quick return for void matrix (Y3K safe) |
* Quick return for void matrix (Y3K safe) |
* #:) |
* #:) |
IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN |
IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN |
|
IWORK(1:3) = 0 |
|
WORK(1:7) = 0 |
|
RETURN |
|
ENDIF |
* |
* |
* Determine whether the matrix U should be M x N or M x M |
* Determine whether the matrix U should be M x N or M x M |
* |
* |
Line 644
|
Line 646
|
AAPP = ZERO |
AAPP = ZERO |
AAQQ = BIG |
AAQQ = BIG |
DO 4781 p = 1, N |
DO 4781 p = 1, N |
AAPP = DMAX1( AAPP, SVA(p) ) |
AAPP = MAX( AAPP, SVA(p) ) |
IF ( SVA(p) .NE. ZERO ) AAQQ = DMIN1( AAQQ, SVA(p) ) |
IF ( SVA(p) .NE. ZERO ) AAQQ = MIN( AAQQ, SVA(p) ) |
4781 CONTINUE |
4781 CONTINUE |
* |
* |
* Quick return for zero M x N matrix |
* Quick return for zero M x N matrix |
Line 715
|
Line 717
|
IWORK(1) = 0 |
IWORK(1) = 0 |
IWORK(2) = 0 |
IWORK(2) = 0 |
END IF |
END IF |
|
IWORK(3) = 0 |
IF ( ERREST ) WORK(3) = ONE |
IF ( ERREST ) WORK(3) = ONE |
IF ( LSVEC .AND. RSVEC ) THEN |
IF ( LSVEC .AND. RSVEC ) THEN |
WORK(4) = ONE |
WORK(4) = ONE |
Line 749
|
Line 752
|
* in one pass through the vector |
* in one pass through the vector |
WORK(M+N+p) = XSC * SCALEM |
WORK(M+N+p) = XSC * SCALEM |
WORK(N+p) = XSC * (SCALEM*DSQRT(TEMP1)) |
WORK(N+p) = XSC * (SCALEM*DSQRT(TEMP1)) |
AATMAX = DMAX1( AATMAX, WORK(N+p) ) |
AATMAX = MAX( AATMAX, WORK(N+p) ) |
IF (WORK(N+p) .NE. ZERO) AATMIN = DMIN1(AATMIN,WORK(N+p)) |
IF (WORK(N+p) .NE. ZERO) AATMIN = MIN(AATMIN,WORK(N+p)) |
1950 CONTINUE |
1950 CONTINUE |
ELSE |
ELSE |
DO 1904 p = 1, M |
DO 1904 p = 1, M |
WORK(M+N+p) = SCALEM*DABS( A(p,IDAMAX(N,A(p,1),LDA)) ) |
WORK(M+N+p) = SCALEM*DABS( A(p,IDAMAX(N,A(p,1),LDA)) ) |
AATMAX = DMAX1( AATMAX, WORK(M+N+p) ) |
AATMAX = MAX( AATMAX, WORK(M+N+p) ) |
AATMIN = DMIN1( AATMIN, WORK(M+N+p) ) |
AATMIN = MIN( AATMIN, WORK(M+N+p) ) |
1904 CONTINUE |
1904 CONTINUE |
END IF |
END IF |
* |
* |
Line 926
|
Line 929
|
* (eg speed by replacing global with restricted window pivoting, such |
* (eg speed by replacing global with restricted window pivoting, such |
* as in SGEQPX from TOMS # 782). Good results will be obtained using |
* as in SGEQPX from TOMS # 782). Good results will be obtained using |
* SGEQPX with properly (!) chosen numerical parameters. |
* SGEQPX with properly (!) chosen numerical parameters. |
* Any improvement of DGEQP3 improves overal performance of DGEJSV. |
* Any improvement of DGEQP3 improves overall performance of DGEJSV. |
* |
* |
* A * P1 = Q1 * [ R1^t 0]^t: |
* A * P1 = Q1 * [ R1^t 0]^t: |
DO 1963 p = 1, N |
DO 1963 p = 1, N |
Line 947
|
Line 950
|
IF ( L2ABER ) THEN |
IF ( L2ABER ) THEN |
* Standard absolute error bound suffices. All sigma_i with |
* Standard absolute error bound suffices. All sigma_i with |
* sigma_i < N*EPSLN*||A|| are flushed to zero. This is an |
* sigma_i < N*EPSLN*||A|| are flushed to zero. This is an |
* agressive enforcement of lower numerical rank by introducing a |
* aggressive enforcement of lower numerical rank by introducing a |
* backward error of the order of N*EPSLN*||A||. |
* backward error of the order of N*EPSLN*||A||. |
TEMP1 = DSQRT(DBLE(N))*EPSLN |
TEMP1 = DSQRT(DBLE(N))*EPSLN |
DO 3001 p = 2, N |
DO 3001 p = 2, N |
Line 959
|
Line 962
|
3001 CONTINUE |
3001 CONTINUE |
3002 CONTINUE |
3002 CONTINUE |
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 aggressive). |
* 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. |
Line 994
|
Line 997
|
MAXPRJ = ONE |
MAXPRJ = ONE |
DO 3051 p = 2, N |
DO 3051 p = 2, N |
TEMP1 = DABS(A(p,p)) / SVA(IWORK(p)) |
TEMP1 = DABS(A(p,p)) / SVA(IWORK(p)) |
MAXPRJ = DMIN1( MAXPRJ, TEMP1 ) |
MAXPRJ = MIN( MAXPRJ, TEMP1 ) |
3051 CONTINUE |
3051 CONTINUE |
IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE. |
IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE. |
END IF |
END IF |
Line 1052
|
Line 1055
|
* Singular Values only |
* Singular Values only |
* |
* |
* .. transpose A(1:NR,1:N) |
* .. transpose A(1:NR,1:N) |
DO 1946 p = 1, MIN0( N-1, NR ) |
DO 1946 p = 1, MIN( N-1, NR ) |
CALL DCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 ) |
CALL DCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 ) |
1946 CONTINUE |
1946 CONTINUE |
* |
* |
Line 1288
|
Line 1291
|
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 opinion on the condition number |
* .. then assume worst case scenario |
* .. then assume worst case scenario |
* R1 is OK for inverse <=> CONDR1 .LT. DBLE(N) |
* R1 is OK for inverse <=> CONDR1 .LT. DBLE(N) |
* more conservative <=> CONDR1 .LT. DSQRT(DBLE(N)) |
* more conservative <=> CONDR1 .LT. DSQRT(DBLE(N)) |
Line 1308
|
Line 1311
|
XSC = DSQRT(SMALL)/EPSLN |
XSC = DSQRT(SMALL)/EPSLN |
DO 3959 p = 2, NR |
DO 3959 p = 2, NR |
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 * MIN(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 |
Line 1329
|
Line 1332
|
ELSE |
ELSE |
* |
* |
* .. ill-conditioned case: second QRF with pivoting |
* .. ill-conditioned case: second QRF with pivoting |
* Note that windowed pivoting would be equaly good |
* Note that windowed pivoting would be equally good |
* numerically, and more run-time efficient. So, in |
* numerically, and more run-time efficient. So, in |
* an optimal implementation, the next call to DGEQP3 |
* an optimal implementation, the next call to DGEQP3 |
* should be replaced with eg. CALL SGEQPX (ACM TOMS #782) |
* should be replaced with eg. CALL SGEQPX (ACM TOMS #782) |
Line 1347
|
Line 1350
|
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 * MIN(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 |
Line 1360
|
Line 1363
|
XSC = DSQRT(SMALL) |
XSC = DSQRT(SMALL) |
DO 8970 p = 2, NR |
DO 8970 p = 2, NR |
DO 8971 q = 1, p - 1 |
DO 8971 q = 1, p - 1 |
TEMP1 = XSC * DMIN1(DABS(V(p,p)),DABS(V(q,q))) |
TEMP1 = XSC * MIN(DABS(V(p,p)),DABS(V(q,q))) |
V(p,q) = - DSIGN( TEMP1, V(q,p) ) |
V(p,q) = - DSIGN( TEMP1, V(q,p) ) |
8971 CONTINUE |
8971 CONTINUE |
8970 CONTINUE |
8970 CONTINUE |
Line 1382
|
Line 1385
|
* |
* |
IF ( CONDR2 .GE. COND_OK ) THEN |
IF ( CONDR2 .GE. COND_OK ) THEN |
* .. save the Householder vectors used for Q3 |
* .. save the Householder vectors used for Q3 |
* (this overwrittes the copy of R2, as it will not be |
* (this overwrites the copy of R2, as it will not be |
* needed in this branch, but it does not overwritte the |
* needed in this branch, but it does not overwritte the |
* Huseholder vectors of Q2.). |
* Huseholder vectors of Q2.). |
CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N ) |
CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N ) |
Line 1632
|
Line 1635
|
* |
* |
* This branch deploys a preconditioned Jacobi SVD with explicitly |
* This branch deploys a preconditioned Jacobi SVD with explicitly |
* accumulated rotations. It is included as optional, mainly for |
* accumulated rotations. It is included as optional, mainly for |
* experimental purposes. It does perfom well, and can also be used. |
* experimental purposes. It does perform well, and can also be used. |
* In this implementation, this branch will be automatically activated |
* In this implementation, this branch will be automatically activated |
* if the condition number sigma_max(A) / sigma_min(A) is predicted |
* if the condition number sigma_max(A) / sigma_min(A) is predicted |
* to be greater than the overflow threshold. This is because the |
* to be greater than the overflow threshold. This is because the |
Line 1671
|
Line 1674
|
XSC = DSQRT(SMALL/EPSLN) |
XSC = DSQRT(SMALL/EPSLN) |
DO 9970 q = 2, NR |
DO 9970 q = 2, NR |
DO 9971 p = 1, q - 1 |
DO 9971 p = 1, q - 1 |
TEMP1 = XSC * DMIN1(DABS(U(p,p)),DABS(U(q,q))) |
TEMP1 = XSC * MIN(DABS(U(p,p)),DABS(U(q,q))) |
U(p,q) = - DSIGN( TEMP1, U(q,p) ) |
U(p,q) = - DSIGN( TEMP1, U(q,p) ) |
9971 CONTINUE |
9971 CONTINUE |
9970 CONTINUE |
9970 CONTINUE |