version 1.6, 2011/07/22 07:38:05
|
version 1.7, 2011/11/21 20:42:52
|
Line 1
|
Line 1
|
SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, |
*> \brief \b DGESVJ |
$ LDV, WORK, LWORK, INFO ) |
* |
|
* =========== DOCUMENTATION =========== |
* |
* |
* -- LAPACK routine (version 3.3.1) -- |
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
* |
* |
* -- Contributed by Zlatko Drmac of the University of Zagreb and -- |
*> \htmlonly |
* -- Kresimir Veselic of the Fernuniversitaet Hagen -- |
*> Download DGESVJ + dependencies |
* -- April 2011 -- |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvj.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvj.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvj.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, |
|
* LDV, WORK, LWORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* INTEGER INFO, LDA, LDV, LWORK, M, MV, N |
|
* CHARACTER*1 JOBA, JOBU, JOBV |
|
* .. |
|
* .. Array Arguments .. |
|
* DOUBLE PRECISION A( LDA, * ), SVA( N ), V( LDV, * ), |
|
* $ WORK( LWORK ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DGESVJ computes the singular value decomposition (SVD) of a real |
|
*> M-by-N matrix A, where M >= N. The SVD of A is written as |
|
*> [++] [xx] [x0] [xx] |
|
*> A = U * SIGMA * V^t, [++] = [xx] * [ox] * [xx] |
|
*> [++] [xx] |
|
*> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal |
|
*> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements |
|
*> of SIGMA are the singular values of A. The columns of U and V are the |
|
*> left and the right singular vectors of A, respectively. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] JOBA |
|
*> \verbatim |
|
*> JOBA is CHARACTER* 1 |
|
*> Specifies the structure of A. |
|
*> = 'L': The input matrix A is lower triangular; |
|
*> = 'U': The input matrix A is upper triangular; |
|
*> = 'G': The input matrix A is general M-by-N matrix, M >= N. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] JOBU |
|
*> \verbatim |
|
*> JOBU is CHARACTER*1 |
|
*> Specifies whether to compute the left singular vectors |
|
*> (columns of U): |
|
*> = 'U': The left singular vectors corresponding to the nonzero |
|
*> singular values are computed and returned in the leading |
|
*> columns of A. See more details in the description of A. |
|
*> The default numerical orthogonality threshold is set to |
|
*> approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E'). |
|
*> = 'C': Analogous to JOBU='U', except that user can control the |
|
*> level of numerical orthogonality of the computed left |
|
*> singular vectors. TOL can be set to TOL = CTOL*EPS, where |
|
*> CTOL is given on input in the array WORK. |
|
*> No CTOL smaller than ONE is allowed. CTOL greater |
|
*> than 1 / EPS is meaningless. The option 'C' |
|
*> can be used if M*EPS is satisfactory orthogonality |
|
*> of the computed left singular vectors, so CTOL=M could |
|
*> save few sweeps of Jacobi rotations. |
|
*> See the descriptions of A and WORK(1). |
|
*> = 'N': The matrix U is not computed. However, see the |
|
*> description of A. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] JOBV |
|
*> \verbatim |
|
*> JOBV is CHARACTER*1 |
|
*> Specifies whether to compute the right singular vectors, that |
|
*> is, the matrix V: |
|
*> = 'V' : the matrix V is computed and returned in the array V |
|
*> = 'A' : the Jacobi rotations are applied to the MV-by-N |
|
*> array V. In other words, the right singular vector |
|
*> matrix V is not computed explicitly, instead it is |
|
*> applied to an MV-by-N matrix initially stored in the |
|
*> first MV rows of V. |
|
*> = 'N' : the matrix V is not computed and the array V is not |
|
*> referenced |
|
*> \endverbatim |
|
*> |
|
*> \param[in] M |
|
*> \verbatim |
|
*> M is INTEGER |
|
*> The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] N |
|
*> \verbatim |
|
*> N is INTEGER |
|
*> The number of columns of the input matrix A. |
|
*> M >= N >= 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] A |
|
*> \verbatim |
|
*> A is DOUBLE PRECISION array, dimension (LDA,N) |
|
*> On entry, the M-by-N matrix A. |
|
*> On exit : |
|
*> If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' : |
|
*> If INFO .EQ. 0 : |
|
*> RANKA orthonormal columns of U are returned in the |
|
*> leading RANKA columns of the array A. Here RANKA <= N |
|
*> is the number of computed singular values of A that are |
|
*> above the underflow threshold DLAMCH('S'). The singular |
|
*> vectors corresponding to underflowed or zero singular |
|
*> values are not computed. The value of RANKA is returned |
|
*> in the array WORK as RANKA=NINT(WORK(2)). Also see the |
|
*> descriptions of SVA and WORK. The computed columns of U |
|
*> are mutually numerically orthogonal up to approximately |
|
*> TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'), |
|
*> see the description of JOBU. |
|
*> If INFO .GT. 0 : |
|
*> the procedure DGESVJ did not converge in the given number |
|
*> of iterations (sweeps). In that case, the computed |
|
*> columns of U may not be orthogonal up to TOL. The output |
|
*> U (stored in A), SIGMA (given by the computed singular |
|
*> values in SVA(1:N)) and V is still a decomposition of the |
|
*> input matrix A in the sense that the residual |
|
*> ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small. |
|
*> |
|
*> If JOBU .EQ. 'N' : |
|
*> If INFO .EQ. 0 : |
|
*> Note that the left singular vectors are 'for free' in the |
|
*> one-sided Jacobi SVD algorithm. However, if only the |
|
*> singular values are needed, the level of numerical |
|
*> orthogonality of U is not an issue and iterations are |
|
*> stopped when the columns of the iterated matrix are |
|
*> numerically orthogonal up to approximately M*EPS. Thus, |
|
*> on exit, A contains the columns of U scaled with the |
|
*> corresponding singular values. |
|
*> If INFO .GT. 0 : |
|
*> the procedure DGESVJ did not converge in the given number |
|
*> of iterations (sweeps). |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDA |
|
*> \verbatim |
|
*> LDA is INTEGER |
|
*> The leading dimension of the array A. LDA >= max(1,M). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] SVA |
|
*> \verbatim |
|
*> SVA is DOUBLE PRECISION array, dimension (N) |
|
*> On exit : |
|
*> If INFO .EQ. 0 : |
|
*> depending on the value SCALE = WORK(1), we have: |
|
*> If SCALE .EQ. ONE : |
|
*> SVA(1:N) contains the computed singular values of A. |
|
*> During the computation SVA contains the Euclidean column |
|
*> norms of the iterated matrices in the array A. |
|
*> If SCALE .NE. ONE : |
|
*> The singular values of A are SCALE*SVA(1:N), and this |
|
*> factored representation is due to the fact that some of the |
|
*> singular values of A might underflow or overflow. |
|
*> If INFO .GT. 0 : |
|
*> the procedure DGESVJ did not converge in the given number of |
|
*> iterations (sweeps) and SCALE*SVA(1:N) may not be accurate. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] MV |
|
*> \verbatim |
|
*> MV is INTEGER |
|
*> If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ |
|
*> is applied to the first MV rows of V. See the description of JOBV. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] V |
|
*> \verbatim |
|
*> V is DOUBLE PRECISION array, dimension (LDV,N) |
|
*> If JOBV = 'V', then V contains on exit the N-by-N matrix of |
|
*> the right singular vectors; |
|
*> If JOBV = 'A', then V contains the product of the computed right |
|
*> singular vector matrix and the initial matrix in |
|
*> the array V. |
|
*> If JOBV = 'N', then V is not referenced. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDV |
|
*> \verbatim |
|
*> LDV is INTEGER |
|
*> The leading dimension of the array V, LDV .GE. 1. |
|
*> If JOBV .EQ. 'V', then LDV .GE. max(1,N). |
|
*> If JOBV .EQ. 'A', then LDV .GE. max(1,MV) . |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] WORK |
|
*> \verbatim |
|
*> WORK is DOUBLE PRECISION array, dimension max(4,M+N). |
|
*> On entry : |
|
*> If JOBU .EQ. 'C' : |
|
*> WORK(1) = CTOL, where CTOL defines the threshold for convergence. |
|
*> The process stops if all columns of A are mutually |
|
*> orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). |
|
*> It is required that CTOL >= ONE, i.e. it is not |
|
*> allowed to force the routine to obtain orthogonality |
|
*> below EPS. |
|
*> On exit : |
|
*> WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) |
|
*> are the computed singular values of A. |
|
*> (See description of SVA().) |
|
*> WORK(2) = NINT(WORK(2)) is the number of the computed nonzero |
|
*> singular values. |
|
*> WORK(3) = NINT(WORK(3)) is the number of the computed singular |
|
*> values that are larger than the underflow threshold. |
|
*> WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi |
|
*> rotations needed for numerical convergence. |
|
*> WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. |
|
*> This is useful information in cases when DGESVJ did |
|
*> not converge, as it can be used to estimate whether |
|
*> the output is stil useful and for post festum analysis. |
|
*> WORK(6) = the largest absolute value over all sines of the |
|
*> Jacobi rotation angles in the last sweep. It can be |
|
*> useful for a post festum analysis. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LWORK |
|
*> \verbatim |
|
*> LWORK is INTEGER |
|
*> length of WORK, WORK >= MAX(6,M+N) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0 : successful exit. |
|
*> < 0 : if INFO = -i, then the i-th argument had an illegal value |
|
*> > 0 : DGESVJ did not converge in the maximal allowed number (30) |
|
*> of sweeps. The output may still be useful. See the |
|
*> description of WORK. |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date November 2011 |
|
* |
|
*> \ingroup doubleGEcomputational |
|
* |
|
*> \par Further Details: |
|
* ===================== |
|
*> |
|
*> \verbatim |
|
*> |
|
*> The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane |
|
*> rotations. The rotations are implemented as fast scaled rotations of |
|
*> Anda and Park [1]. In the case of underflow of the Jacobi angle, a |
|
*> modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses |
|
*> column interchanges of de Rijk [2]. The relative accuracy of the computed |
|
*> singular values and the accuracy of the computed singular vectors (in |
|
*> angle metric) is as guaranteed by the theory of Demmel and Veselic [3]. |
|
*> The condition number that determines the accuracy in the full rank case |
|
*> is essentially min_{D=diag} kappa(A*D), where kappa(.) is the |
|
*> spectral condition number. The best performance of this Jacobi SVD |
|
*> procedure is achieved if used in an accelerated version of Drmac and |
|
*> Veselic [5,6], and it is the kernel routine in the SIGMA library [7]. |
|
*> Some tunning parameters (marked with [TP]) are available for the |
|
*> implementer. |
|
*> The computational range for the nonzero singular values is the machine |
|
*> number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even |
|
*> denormalized singular values can be computed with the corresponding |
|
*> gradual loss of accurate digits. |
|
*> \endverbatim |
|
* |
|
*> \par Contributors: |
|
* ================== |
|
*> |
|
*> \verbatim |
|
*> |
|
*> ============ |
|
*> |
|
*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) |
|
*> \endverbatim |
|
* |
|
*> \par References: |
|
* ================ |
|
*> |
|
*> \verbatim |
|
*> |
|
*> [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling. |
|
*> SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174. |
|
*> [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the |
|
*> singular value decomposition on a vector computer. |
|
*> SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371. |
|
*> [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR. |
|
*> [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular |
|
*> value computation in floating point arithmetic. |
|
*> SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222. |
|
*> [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. |
|
*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. |
|
*> LAPACK Working note 169. |
|
*> [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. |
|
*> SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. |
|
*> LAPACK Working note 170. |
|
*> [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, |
|
*> QSVD, (H,K)-SVD computations. |
|
*> Department of Mathematics, University of Zagreb, 2008. |
|
*> \endverbatim |
|
* |
|
*> \par Bugs, examples and comments: |
|
* ================================= |
|
*> |
|
*> \verbatim |
|
*> =========================== |
|
*> Please report all bugs and send interesting test examples and comments to |
|
*> drmac@math.hr. Thank you. |
|
*> \endverbatim |
|
*> |
|
* ===================================================================== |
|
SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V, |
|
$ LDV, WORK, LWORK, INFO ) |
* |
* |
|
* -- LAPACK computational routine (version 3.4.0) -- |
* -- 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 |
* |
* |
* This routine is also part of SIGMA (version 1.23, October 23. 2008.) |
|
* SIGMA is a library of algorithms for highly accurate algorithms for |
|
* computation of SVD, PSVD, QSVD, (H,K)-SVD, and for solution of the |
|
* eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0. |
|
* |
|
IMPLICIT NONE |
|
* .. |
|
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
INTEGER INFO, LDA, LDV, LWORK, M, MV, N |
INTEGER INFO, LDA, LDV, LWORK, M, MV, N |
CHARACTER*1 JOBA, JOBU, JOBV |
CHARACTER*1 JOBA, JOBU, JOBV |
Line 26
|
Line 349
|
$ WORK( LWORK ) |
$ WORK( LWORK ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* DGESVJ computes the singular value decomposition (SVD) of a real |
|
* M-by-N matrix A, where M >= N. The SVD of A is written as |
|
* [++] [xx] [x0] [xx] |
|
* A = U * SIGMA * V^t, [++] = [xx] * [ox] * [xx] |
|
* [++] [xx] |
|
* where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal |
|
* matrix, and V is an N-by-N orthogonal matrix. The diagonal elements |
|
* of SIGMA are the singular values of A. The columns of U and V are the |
|
* left and the right singular vectors of A, respectively. |
|
* |
|
* Further Details |
|
* ~~~~~~~~~~~~~~~ |
|
* The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane |
|
* rotations. The rotations are implemented as fast scaled rotations of |
|
* Anda and Park [1]. In the case of underflow of the Jacobi angle, a |
|
* modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses |
|
* column interchanges of de Rijk [2]. The relative accuracy of the computed |
|
* singular values and the accuracy of the computed singular vectors (in |
|
* angle metric) is as guaranteed by the theory of Demmel and Veselic [3]. |
|
* The condition number that determines the accuracy in the full rank case |
|
* is essentially min_{D=diag} kappa(A*D), where kappa(.) is the |
|
* spectral condition number. The best performance of this Jacobi SVD |
|
* procedure is achieved if used in an accelerated version of Drmac and |
|
* Veselic [5,6], and it is the kernel routine in the SIGMA library [7]. |
|
* Some tunning parameters (marked with [TP]) are available for the |
|
* implementer. |
|
* The computational range for the nonzero singular values is the machine |
|
* number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even |
|
* denormalized singular values can be computed with the corresponding |
|
* gradual loss of accurate digits. |
|
* |
|
* Contributors |
|
* ~~~~~~~~~~~~ |
|
* Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) |
|
* |
|
* References |
|
* ~~~~~~~~~~ |
|
* [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling. |
|
* SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174. |
|
* [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the |
|
* singular value decomposition on a vector computer. |
|
* SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371. |
|
* [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR. |
|
* [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular |
|
* value computation in floating point arithmetic. |
|
* SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222. |
|
* [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I. |
|
* SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342. |
|
* LAPACK Working note 169. |
|
* [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II. |
|
* SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362. |
|
* LAPACK Working note 170. |
|
* [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, |
|
* QSVD, (H,K)-SVD computations. |
|
* Department of Mathematics, University of Zagreb, 2008. |
|
* |
|
* Bugs, Examples and Comments |
|
* ~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|
* Please report all bugs and send interesting test examples and comments to |
|
* drmac@math.hr. Thank you. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* JOBA (input) CHARACTER* 1 |
|
* Specifies the structure of A. |
|
* = 'L': The input matrix A is lower triangular; |
|
* = 'U': The input matrix A is upper triangular; |
|
* = 'G': The input matrix A is general M-by-N matrix, M >= N. |
|
* |
|
* JOBU (input) CHARACTER*1 |
|
* Specifies whether to compute the left singular vectors |
|
* (columns of U): |
|
* = 'U': The left singular vectors corresponding to the nonzero |
|
* singular values are computed and returned in the leading |
|
* columns of A. See more details in the description of A. |
|
* The default numerical orthogonality threshold is set to |
|
* approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E'). |
|
* = 'C': Analogous to JOBU='U', except that user can control the |
|
* level of numerical orthogonality of the computed left |
|
* singular vectors. TOL can be set to TOL = CTOL*EPS, where |
|
* CTOL is given on input in the array WORK. |
|
* No CTOL smaller than ONE is allowed. CTOL greater |
|
* than 1 / EPS is meaningless. The option 'C' |
|
* can be used if M*EPS is satisfactory orthogonality |
|
* of the computed left singular vectors, so CTOL=M could |
|
* save few sweeps of Jacobi rotations. |
|
* See the descriptions of A and WORK(1). |
|
* = 'N': The matrix U is not computed. However, see the |
|
* description of A. |
|
* |
|
* JOBV (input) CHARACTER*1 |
|
* Specifies whether to compute the right singular vectors, that |
|
* is, the matrix V: |
|
* = 'V' : the matrix V is computed and returned in the array V |
|
* = 'A' : the Jacobi rotations are applied to the MV-by-N |
|
* array V. In other words, the right singular vector |
|
* matrix V is not computed explicitly, instead it is |
|
* applied to an MV-by-N matrix initially stored in the |
|
* first MV rows of V. |
|
* = 'N' : the matrix V is not computed and the array V is not |
|
* referenced |
|
* |
|
* M (input) INTEGER |
|
* The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0. |
|
* |
|
* N (input) INTEGER |
|
* The number of columns of the input matrix A. |
|
* M >= N >= 0. |
|
* |
|
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) |
|
* On entry, the M-by-N matrix A. |
|
* On exit : |
|
* If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' : |
|
* If INFO .EQ. 0 : |
|
* RANKA orthonormal columns of U are returned in the |
|
* leading RANKA columns of the array A. Here RANKA <= N |
|
* is the number of computed singular values of A that are |
|
* above the underflow threshold DLAMCH('S'). The singular |
|
* vectors corresponding to underflowed or zero singular |
|
* values are not computed. The value of RANKA is returned |
|
* in the array WORK as RANKA=NINT(WORK(2)). Also see the |
|
* descriptions of SVA and WORK. The computed columns of U |
|
* are mutually numerically orthogonal up to approximately |
|
* TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'), |
|
* see the description of JOBU. |
|
* If INFO .GT. 0 : |
|
* the procedure DGESVJ did not converge in the given number |
|
* of iterations (sweeps). In that case, the computed |
|
* columns of U may not be orthogonal up to TOL. The output |
|
* U (stored in A), SIGMA (given by the computed singular |
|
* values in SVA(1:N)) and V is still a decomposition of the |
|
* input matrix A in the sense that the residual |
|
* ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small. |
|
* |
|
* If JOBU .EQ. 'N' : |
|
* If INFO .EQ. 0 : |
|
* Note that the left singular vectors are 'for free' in the |
|
* one-sided Jacobi SVD algorithm. However, if only the |
|
* singular values are needed, the level of numerical |
|
* orthogonality of U is not an issue and iterations are |
|
* stopped when the columns of the iterated matrix are |
|
* numerically orthogonal up to approximately M*EPS. Thus, |
|
* on exit, A contains the columns of U scaled with the |
|
* corresponding singular values. |
|
* If INFO .GT. 0 : |
|
* the procedure DGESVJ did not converge in the given number |
|
* of iterations (sweeps). |
|
* |
|
* LDA (input) INTEGER |
|
* The leading dimension of the array A. LDA >= max(1,M). |
|
* |
|
* SVA (workspace/output) DOUBLE PRECISION array, dimension (N) |
|
* On exit : |
|
* If INFO .EQ. 0 : |
|
* depending on the value SCALE = WORK(1), we have: |
|
* If SCALE .EQ. ONE : |
|
* SVA(1:N) contains the computed singular values of A. |
|
* During the computation SVA contains the Euclidean column |
|
* norms of the iterated matrices in the array A. |
|
* If SCALE .NE. ONE : |
|
* The singular values of A are SCALE*SVA(1:N), and this |
|
* factored representation is due to the fact that some of the |
|
* singular values of A might underflow or overflow. |
|
* If INFO .GT. 0 : |
|
* the procedure DGESVJ did not converge in the given number of |
|
* iterations (sweeps) and SCALE*SVA(1:N) may not be accurate. |
|
* |
|
* MV (input) INTEGER |
|
* If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ |
|
* is applied to the first MV rows of V. See the description of JOBV. |
|
* |
|
* V (input/output) DOUBLE PRECISION array, dimension (LDV,N) |
|
* If JOBV = 'V', then V contains on exit the N-by-N matrix of |
|
* the right singular vectors; |
|
* If JOBV = 'A', then V contains the product of the computed right |
|
* singular vector matrix and the initial matrix in |
|
* the array V. |
|
* If JOBV = 'N', then V is not referenced. |
|
* |
|
* LDV (input) INTEGER |
|
* The leading dimension of the array V, LDV .GE. 1. |
|
* If JOBV .EQ. 'V', then LDV .GE. max(1,N). |
|
* If JOBV .EQ. 'A', then LDV .GE. max(1,MV) . |
|
* |
|
* WORK (input/workspace/output) DOUBLE PRECISION array, dimension max(4,M+N). |
|
* On entry : |
|
* If JOBU .EQ. 'C' : |
|
* WORK(1) = CTOL, where CTOL defines the threshold for convergence. |
|
* The process stops if all columns of A are mutually |
|
* orthogonal up to CTOL*EPS, EPS=DLAMCH('E'). |
|
* It is required that CTOL >= ONE, i.e. it is not |
|
* allowed to force the routine to obtain orthogonality |
|
* below EPS. |
|
* On exit : |
|
* WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N) |
|
* are the computed singular values of A. |
|
* (See description of SVA().) |
|
* WORK(2) = NINT(WORK(2)) is the number of the computed nonzero |
|
* singular values. |
|
* WORK(3) = NINT(WORK(3)) is the number of the computed singular |
|
* values that are larger than the underflow threshold. |
|
* WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi |
|
* rotations needed for numerical convergence. |
|
* WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep. |
|
* This is useful information in cases when DGESVJ did |
|
* not converge, as it can be used to estimate whether |
|
* the output is stil useful and for post festum analysis. |
|
* WORK(6) = the largest absolute value over all sines of the |
|
* Jacobi rotation angles in the last sweep. It can be |
|
* useful for a post festum analysis. |
|
* |
|
* LWORK (input) INTEGER |
|
* length of WORK, WORK >= MAX(6,M+N) |
|
* |
|
* INFO (output) INTEGER |
|
* = 0 : successful exit. |
|
* < 0 : if INFO = -i, then the i-th argument had an illegal value |
|
* > 0 : DGESVJ did not converge in the maximal allowed number (30) |
|
* of sweeps. The output may still be useful. See the |
|
* description of WORK. |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Local Parameters .. |
* .. Local Parameters .. |