diff --git a/lapack-netlib/zgejsv.f b/lapack-netlib/zgejsv.f deleted file mode 100644 index 5fe899e50..000000000 --- a/lapack-netlib/zgejsv.f +++ /dev/null @@ -1,2234 +0,0 @@ -*> \brief \b ZGEJSV -* -* =========== DOCUMENTATION =========== -* -* Online html documentation available at -* http://www.netlib.org/lapack/explore-html/ -* -*> \htmlonly -*> Download ZGEJSV + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -*> \endhtmlonly -* -* Definition: -* =========== -* -* SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, -* M, N, A, LDA, SVA, U, LDU, V, LDV, -* CWORK, LWORK, RWORK, LRWORK, IWORK, INFO ) -* -* .. Scalar Arguments .. -* IMPLICIT NONE -* INTEGER INFO, LDA, LDU, LDV, LWORK, M, N -* .. -* .. Array Arguments .. -* COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( LWORK ) -* DOUBLE PRECISION SVA( N ), RWORK( LRWORK ) -* INTEGER IWORK( * ) -* CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV -* .. -* -* -*> \par Purpose: -* ============= -*> -*> \verbatim -*> -*> ZGEJSV computes the singular value decomposition (SVD) of a complex M-by-N -*> matrix [A], where M >= N. The SVD of [A] is written as -*> -*> [A] = [U] * [SIGMA] * [V]^*, -*> -*> where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N -*> diagonal elements, [U] is an M-by-N (or M-by-M) unitary matrix, and -*> [V] is an N-by-N unitary 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. The matrices [U] and [V] -*> are computed and stored in the arrays U and V, respectively. The diagonal -*> of [SIGMA] is computed and stored in the array SVA. -*> \endverbatim -*> -*> Arguments: -*> ========== -*> -*> \param[in] JOBA -*> \verbatim -*> JOBA is CHARACTER*1 -*> Specifies the level of accuracy: -*> = 'C': This option works well (high relative accuracy) if A = B * D, -*> with well-conditioned B and arbitrary diagonal matrix D. -*> The accuracy cannot be spoiled by COLUMN scaling. The -*> accuracy of the computed output depends on the condition of -*> B, and the procedure aims at the best theoretical accuracy. -*> The relative error max_{i=1:N}|d sigma_i| / sigma_i is -*> bounded by f(M,N)*epsilon* cond(B), independent of D. -*> The input matrix is preprocessed with the QRF with column -*> pivoting. This initial preprocessing and preconditioning by -*> a rank revealing QR factorization is common for all values of -*> JOBA. Additional actions are specified as follows: -*> = 'E': Computation as with 'C' with an additional estimate of the -*> condition number of B. It provides a realistic error bound. -*> = 'F': If A = D1 * C * D2 with ill-conditioned diagonal scalings -*> D1, D2, and well-conditioned matrix C, this option gives -*> higher accuracy than the 'C' option. If the structure of the -*> input matrix is not known, and relative accuracy is -*> desirable, then this option is advisable. The input matrix A -*> is preprocessed with QR factorization with FULL (row and -*> column) pivoting. -*> = 'G': Computation as with 'F' with an additional estimate of the -*> condition number of B, where A=B*D. If A has heavily weighted -*> rows, then using this condition number gives too pessimistic -*> error bound. -*> = 'A': Small singular values are not well determined by the data -*> and are considered as noisy; the matrix is treated as -*> numerically rank deficient. The error in the computed -*> singular values is bounded by f(m,n)*epsilon*||A||. -*> The computed SVD A = U * S * V^* restores A up to -*> f(m,n)*epsilon*||A||. -*> This gives the procedure the licence to discard (set to zero) -*> all singular values below N*epsilon*||A||. -*> = 'R': Similar as in 'A'. Rank revealing property of the initial -*> QR factorization is used do reveal (using triangular factor) -*> a gap sigma_{r+1} < epsilon * sigma_r in which case the -*> numerical RANK is declared to be r. The SVD is computed with -*> absolute error bounds, but more accurately than with 'A'. -*> \endverbatim -*> -*> \param[in] JOBU -*> \verbatim -*> JOBU is CHARACTER*1 -*> Specifies whether to compute the columns of U: -*> = 'U': N columns of U are returned in the array U. -*> = 'F': full set of M left sing. vectors is returned in the array U. -*> = 'W': U may be used as workspace of length M*N. See the description -*> of U. -*> = 'N': U is not computed. -*> \endverbatim -*> -*> \param[in] JOBV -*> \verbatim -*> JOBV is CHARACTER*1 -*> Specifies whether to compute the matrix V: -*> = 'V': N columns of V are returned in the array V; Jacobi rotations -*> are not explicitly accumulated. -*> = 'J': N columns of V are returned in the array V, but they are -*> computed as the product of Jacobi rotations, if JOBT = 'N'. -*> = 'W': V may be used as workspace of length N*N. See the description -*> of V. -*> = 'N': V is not computed. -*> \endverbatim -*> -*> \param[in] JOBR -*> \verbatim -*> JOBR is CHARACTER*1 -*> Specifies the RANGE for the singular values. Issues the licence to -*> set to zero small positive singular values if they are outside -*> specified range. If A .NE. 0 is scaled so that the largest singular -*> value of c*A is around SQRT(BIG), BIG=DLAMCH('O'), then JOBR issues -*> the licence to kill columns of A whose norm in c*A is less than -*> SQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN, -*> where SFMIN=DLAMCH('S'), EPSLN=DLAMCH('E'). -*> = 'N': Do not kill small columns of c*A. This option assumes that -*> BLAS and QR factorizations and triangular solvers are -*> implemented to work in that range. If the condition of A -*> is greater than BIG, use ZGESVJ. -*> = 'R': RESTRICTED range for sigma(c*A) is [SQRT(SFMIN), SQRT(BIG)] -*> (roughly, as described above). This option is recommended. -*> =========================== -*> For computing the singular values in the FULL range [SFMIN,BIG] -*> use ZGESVJ. -*> \endverbatim -*> -*> \param[in] JOBT -*> \verbatim -*> JOBT is CHARACTER*1 -*> If the matrix is square then the procedure may determine to use -*> transposed A if A^* seems to be better with respect to convergence. -*> If the matrix is not square, JOBT is ignored. -*> The decision is based on two values of entropy over the adjoint -*> orbit of A^* * A. See the descriptions of RWORK(6) and RWORK(7). -*> = 'T': transpose if entropy test indicates possibly faster -*> convergence of Jacobi process if A^* is taken as input. If A is -*> replaced with A^*, then the row pivoting is included automatically. -*> = 'N': do not speculate. -*> The option 'T' can be used to compute only the singular values, or -*> the full SVD (U, SIGMA and V). For only one set of singular vectors -*> (U or V), the caller should provide both U and V, as one of the -*> matrices is used as workspace if the matrix A is transposed. -*> The implementer can easily remove this constraint and make the -*> code more complicated. See the descriptions of U and V. -*> In general, this option is considered experimental, and 'N'; should -*> be preferred. This is subject to changes in the future. -*> \endverbatim -*> -*> \param[in] JOBP -*> \verbatim -*> JOBP is CHARACTER*1 -*> Issues the licence to introduce structured perturbations to drown -*> denormalized numbers. This licence should be active if the -*> denormals are poorly implemented, causing slow computation, -*> especially in cases of fast convergence (!). For details see [1,2]. -*> For the sake of simplicity, this perturbations are included only -*> when the full SVD or only the singular values are requested. The -*> implementer/user can easily add the perturbation for the cases of -*> computing one set of singular vectors. -*> = 'P': introduce perturbation -*> = 'N': do not perturb -*> \endverbatim -*> -*> \param[in] M -*> \verbatim -*> M is INTEGER -*> The number of rows of the input matrix A. 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 COMPLEX*16 array, dimension (LDA,N) -*> On entry, the M-by-N matrix A. -*> \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, -*> - For RWORK(1)/RWORK(2) = ONE: The singular values of A. During -*> the computation SVA contains Euclidean column norms of the -*> iterated matrices in the array A. -*> - For RWORK(1) .NE. RWORK(2): The singular values of A are -*> (RWORK(1)/RWORK(2)) * SVA(1:N). This factored form is used if -*> sigma_max(A) overflows or if small singular values have been -*> saved from underflow by scaling the input matrix A. -*> - If JOBR='R' then some of the singular values may be returned -*> as exact zeros obtained by "set to zero" because they are -*> below the numerical rank threshold or are denormalized numbers. -*> \endverbatim -*> -*> \param[out] U -*> \verbatim -*> U is COMPLEX*16 array, dimension ( LDU, N ) -*> If JOBU = 'U', then U contains on exit the M-by-N matrix of -*> the left singular vectors. -*> If JOBU = 'F', then U contains on exit the M-by-M matrix of -*> the left singular vectors, including an ONB -*> of the orthogonal complement of the Range(A). -*> If JOBU = 'W' .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N), -*> then U is used as workspace if the procedure -*> replaces A with A^*. In that case, [V] is computed -*> in U as left singular vectors of A^* and then -*> copied back to the V array. This 'W' option is just -*> a reminder to the caller that in this case U is -*> reserved as workspace of length N*N. -*> If JOBU = 'N' U is not referenced, unless JOBT='T'. -*> \endverbatim -*> -*> \param[in] LDU -*> \verbatim -*> LDU is INTEGER -*> The leading dimension of the array U, LDU >= 1. -*> IF JOBU = 'U' or 'F' or 'W', then LDU >= M. -*> \endverbatim -*> -*> \param[out] V -*> \verbatim -*> V is COMPLEX*16 array, dimension ( LDV, N ) -*> If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of -*> the right singular vectors; -*> If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N), -*> then V is used as workspace if the pprocedure -*> replaces A with A^*. In that case, [U] is computed -*> in V as right singular vectors of A^* and then -*> copied back to the U array. This 'W' option is just -*> a reminder to the caller that in this case V is -*> reserved as workspace of length N*N. -*> If JOBV = 'N' V is not referenced, unless JOBT='T'. -*> \endverbatim -*> -*> \param[in] LDV -*> \verbatim -*> LDV is INTEGER -*> The leading dimension of the array V, LDV >= 1. -*> If JOBV = 'V' or 'J' or 'W', then LDV >= N. -*> \endverbatim -*> -*> \param[out] CWORK -*> \verbatim -*> CWORK is COMPLEX*16 array, dimension (MAX(2,LWORK)) -*> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or -*> LRWORK=-1), then on exit CWORK(1) contains the required length of -*> CWORK for the job parameters used in the call. -*> \endverbatim -*> -*> \param[in] LWORK -*> \verbatim -*> LWORK is INTEGER -*> Length of CWORK to confirm proper allocation of workspace. -*> LWORK depends on the job: -*> -*> 1. If only SIGMA is needed ( JOBU = 'N', JOBV = 'N' ) and -*> 1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'): -*> LWORK >= 2*N+1. This is the minimal requirement. -*> ->> For optimal performance (blocked code) the optimal value -*> is LWORK >= N + (N+1)*NB. Here NB is the optimal -*> block size for ZGEQP3 and ZGEQRF. -*> In general, optimal LWORK is computed as -*> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ)). -*> 1.2. .. an estimate of the scaled condition number of A is -*> required (JOBA='E', or 'G'). In this case, LWORK the minimal -*> requirement is LWORK >= N*N + 2*N. -*> ->> For optimal performance (blocked code) the optimal value -*> is LWORK >= max(N+(N+1)*NB, N*N+2*N)=N**2+2*N. -*> In general, the optimal length LWORK is computed as -*> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ), -*> N*N+LWORK(ZPOCON)). -*> 2. If SIGMA and the right singular vectors are needed (JOBV = 'V'), -*> (JOBU = 'N') -*> 2.1 .. no scaled condition estimate requested (JOBE = 'N'): -*> -> the minimal requirement is LWORK >= 3*N. -*> -> For optimal performance, -*> LWORK >= max(N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, -*> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQF, -*> ZUNMLQ. In general, the optimal length LWORK is computed as -*> LWORK >= max(N+LWORK(ZGEQP3), N+LWORK(ZGESVJ), -*> N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)). -*> 2.2 .. an estimate of the scaled condition number of A is -*> required (JOBA='E', or 'G'). -*> -> the minimal requirement is LWORK >= 3*N. -*> -> For optimal performance, -*> LWORK >= max(N+(N+1)*NB, 2*N,2*N+N*NB)=2*N+N*NB, -*> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQF, -*> ZUNMLQ. In general, the optimal length LWORK is computed as -*> LWORK >= max(N+LWORK(ZGEQP3), LWORK(ZPOCON), N+LWORK(ZGESVJ), -*> N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)). -*> 3. If SIGMA and the left singular vectors are needed -*> 3.1 .. no scaled condition estimate requested (JOBE = 'N'): -*> -> the minimal requirement is LWORK >= 3*N. -*> -> For optimal performance: -*> if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, -*> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR. -*> In general, the optimal length LWORK is computed as -*> LWORK >= max(N+LWORK(ZGEQP3), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)). -*> 3.2 .. an estimate of the scaled condition number of A is -*> required (JOBA='E', or 'G'). -*> -> the minimal requirement is LWORK >= 3*N. -*> -> For optimal performance: -*> if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB, -*> where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR. -*> In general, the optimal length LWORK is computed as -*> LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZPOCON), -*> 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)). -*> 4. If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and -*> 4.1. if JOBV = 'V' -*> the minimal requirement is LWORK >= 5*N+2*N*N. -*> 4.2. if JOBV = 'J' the minimal requirement is -*> LWORK >= 4*N+N*N. -*> In both cases, the allocated CWORK can accommodate blocked runs -*> of ZGEQP3, ZGEQRF, ZGELQF, SUNMQR, ZUNMLQ. -*> -*> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or -*> LRWORK=-1), then on exit CWORK(1) contains the optimal and CWORK(2) contains the -*> minimal length of CWORK for the job parameters used in the call. -*> \endverbatim -*> -*> \param[out] RWORK -*> \verbatim -*> RWORK is DOUBLE PRECISION array, dimension (MAX(7,LRWORK)) -*> On exit, -*> RWORK(1) = Determines the scaling factor SCALE = RWORK(2) / RWORK(1) -*> such that SCALE*SVA(1:N) are the computed singular values -*> of A. (See the description of SVA().) -*> RWORK(2) = See the description of RWORK(1). -*> RWORK(3) = SCONDA is an estimate for the condition number of -*> column equilibrated A. (If JOBA = 'E' or 'G') -*> SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1). -*> It is computed using ZPOCON. It holds -*> N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA -*> where R is the triangular factor from the QRF of A. -*> However, if R is truncated and the numerical rank is -*> determined to be strictly smaller than N, SCONDA is -*> returned as -1, thus indicating that the smallest -*> singular values might be lost. -*> -*> If full SVD is needed, the following two condition numbers are -*> useful for the analysis of the algorithm. They are provided for -*> a developer/implementer who is familiar with the details of -*> the method. -*> -*> RWORK(4) = an estimate of the scaled condition number of the -*> triangular factor in the first QR factorization. -*> RWORK(5) = an estimate of the scaled condition number of the -*> triangular factor in the second QR factorization. -*> The following two parameters are computed if JOBT = 'T'. -*> They are provided for a developer/implementer who is familiar -*> with the details of the method. -*> RWORK(6) = the entropy of A^* * A :: this is the Shannon entropy -*> of diag(A^* * A) / Trace(A^* * A) taken as point in the -*> probability simplex. -*> RWORK(7) = the entropy of A * A^*. (See the description of RWORK(6).) -*> If the call to ZGEJSV is a workspace query (indicated by LWORK=-1 or -*> LRWORK=-1), then on exit RWORK(1) contains the required length of -*> RWORK for the job parameters used in the call. -*> \endverbatim -*> -*> \param[in] LRWORK -*> \verbatim -*> LRWORK is INTEGER -*> Length of RWORK to confirm proper allocation of workspace. -*> LRWORK depends on the job: -*> -*> 1. If only the singular values are requested i.e. if -*> LSAME(JOBU,'N') .AND. LSAME(JOBV,'N') -*> then: -*> 1.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), -*> then: LRWORK = max( 7, 2 * M ). -*> 1.2. Otherwise, LRWORK = max( 7, N ). -*> 2. If singular values with the right singular vectors are requested -*> i.e. if -*> (LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) .AND. -*> .NOT.(LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) -*> then: -*> 2.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), -*> then LRWORK = max( 7, 2 * M ). -*> 2.2. Otherwise, LRWORK = max( 7, N ). -*> 3. If singular values with the left singular vectors are requested, i.e. if -*> (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND. -*> .NOT.(LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) -*> then: -*> 3.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), -*> then LRWORK = max( 7, 2 * M ). -*> 3.2. Otherwise, LRWORK = max( 7, N ). -*> 4. If singular values with both the left and the right singular vectors -*> are requested, i.e. if -*> (LSAME(JOBU,'U').OR.LSAME(JOBU,'F')) .AND. -*> (LSAME(JOBV,'V').OR.LSAME(JOBV,'J')) -*> then: -*> 4.1. If LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G'), -*> then LRWORK = max( 7, 2 * M ). -*> 4.2. Otherwise, LRWORK = max( 7, N ). -*> -*> If, on entry, LRWORK = -1 or LWORK=-1, a workspace query is assumed and -*> the length of RWORK is returned in RWORK(1). -*> \endverbatim -*> -*> \param[out] IWORK -*> \verbatim -*> IWORK is INTEGER array, of dimension at least 4, that further depends -*> on the job: -*> -*> 1. If only the singular values are requested then: -*> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) -*> then the length of IWORK is N+M; otherwise the length of IWORK is N. -*> 2. If the singular values and the right singular vectors are requested then: -*> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) -*> then the length of IWORK is N+M; otherwise the length of IWORK is N. -*> 3. If the singular values and the left singular vectors are requested then: -*> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) -*> then the length of IWORK is N+M; otherwise the length of IWORK is N. -*> 4. If the singular values with both the left and the right singular vectors -*> are requested, then: -*> 4.1. If LSAME(JOBV,'J') the length of IWORK is determined as follows: -*> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) -*> then the length of IWORK is N+M; otherwise the length of IWORK is N. -*> 4.2. If LSAME(JOBV,'V') the length of IWORK is determined as follows: -*> If ( LSAME(JOBT,'T') .OR. LSAME(JOBA,'F') .OR. LSAME(JOBA,'G') ) -*> then the length of IWORK is 2*N+M; otherwise the length of IWORK is 2*N. -*> -*> On exit, -*> IWORK(1) = the numerical rank determined after the initial -*> QR factorization with pivoting. See the descriptions -*> of JOBA and JOBR. -*> IWORK(2) = the number of the computed nonzero singular values -*> IWORK(3) = if nonzero, a warning message: -*> If IWORK(3) = 1 then some of the column norms of A -*> were denormalized floats. The requested high accuracy -*> is not warranted by the data. -*> IWORK(4) = 1 or -1. If IWORK(4) = 1, then the procedure used A^* to -*> do the job as specified by the JOB parameters. -*> If the call to ZGEJSV is a workspace query (indicated by LWORK = -1 or -*> LRWORK = -1), then on exit IWORK(1) contains the required length of -*> IWORK for the job parameters used in the call. -*> \endverbatim -*> -*> \param[out] INFO -*> \verbatim -*> INFO is INTEGER -*> < 0: if INFO = -i, then the i-th argument had an illegal value. -*> = 0: successful exit; -*> > 0: ZGEJSV did not converge in the maximal allowed number -*> of sweeps. The computed values may be inaccurate. -*> \endverbatim -* -* Authors: -* ======== -* -*> \author Univ. of Tennessee -*> \author Univ. of California Berkeley -*> \author Univ. of Colorado Denver -*> \author NAG Ltd. -* -*> \ingroup complex16GEsing -* -*> \par Further Details: -* ===================== -*> -*> \verbatim -*> -*> ZGEJSV implements a preconditioned Jacobi SVD algorithm. It uses ZGEQP3, -*> ZGEQRF, and ZGELQF as preprocessors and preconditioners. Optionally, an -*> 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 -*> structure A = D1 * C * D2, where D1, D2 are arbitrarily ill-conditioned -*> diagonal matrices and C is well-conditioned matrix. In that case, complete -*> pivoting in the first QR factorizations provides accuracy dependent on the -*> condition number of C, and independent of D1, D2. Such higher accuracy is -*> not completely understood theoretically, but it works well in practice. -*> Further, if A can be written as A = B*D, with well-conditioned B and some -*> diagonal D, then the high accuracy is guaranteed, both theoretically and -*> in software, independent of D. For more details see [1], [2]. -*> The computational range for the singular values can be the full range -*> ( UNDERFLOW,OVERFLOW ), provided that the machine arithmetic and the BLAS -*> & LAPACK routines called by ZGEJSV are implemented to work in that range. -*> If that is not the case, then the restriction for safe computation with -*> 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 -*> overflow. This code (ZGEJSV) is best used in this restricted range, -*> meaning that singular values of magnitude below ||A||_2 / DLAMCH('O') are -*> returned as zeros. See JOBR for details on this. -*> Further, this implementation is somewhat slower than the one described -*> in [1,2] due to replacement of some non-LAPACK components, and because -*> the choice of some tuning parameters in the iterative part (ZGESVJ) is -*> left to the implementer on a particular machine. -*> The rank revealing QR factorization (in this code: ZGEQP3) should be -*> implemented as in [3]. We have a new version of ZGEQP3 under development -*> that is more robust than the current one in LAPACK, with a cleaner cut in -*> rank deficient cases. It will be available in the SIGMA library [4]. -*> 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 -*> well known trick is not used in ZGEJSV because in some cases heavy row -*> weighting can be treated with complete pivoting. The overhead in cases -*> M much larger than N is then only due to pivoting, but the benefits in -*> terms of accuracy have prevailed. The implementer/user can incorporate -*> this extra QRF step easily. The implementer can also improve data movement -*> (matrix transpose, matrix copy, matrix transposed copy) - this -*> implementation of ZGEJSV uses only the simplest, naive data movement. -*> \endverbatim -* -*> \par Contributor: -* ================== -*> -*> Zlatko Drmac, Department of Mathematics, Faculty of Science, -*> University of Zagreb (Zagreb, Croatia); drmac@math.hr -* -*> \par References: -* ================ -*> -*> \verbatim -*> -*> [1] 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. -*> [2] 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. -*> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR -*> factorization software - a case study. -*> ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28. -*> LAPACK Working note 176. -*> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV, -*> QSVD, (H,K)-SVD computations. -*> Department of Mathematics, University of Zagreb, 2008, 2016. -*> \endverbatim -* -*> \par Bugs, examples and comments: -* ================================= -*> -*> Please report all bugs and send interesting examples and/or comments to -*> drmac@math.hr. Thank you. -*> -* ===================================================================== - SUBROUTINE ZGEJSV( JOBA, JOBU, JOBV, JOBR, JOBT, JOBP, - $ M, N, A, LDA, SVA, U, LDU, V, LDV, - $ CWORK, LWORK, RWORK, LRWORK, IWORK, INFO ) -* -* -- LAPACK computational routine -- -* -- LAPACK is a software package provided by Univ. of Tennessee, -- -* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* -* .. Scalar Arguments .. - IMPLICIT NONE - INTEGER INFO, LDA, LDU, LDV, LWORK, LRWORK, M, N -* .. -* .. Array Arguments .. - COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), - $ CWORK( LWORK ) - DOUBLE PRECISION SVA( N ), RWORK( LRWORK ) - INTEGER IWORK( * ) - CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV -* .. -* -* =========================================================================== -* -* .. Local Parameters .. - DOUBLE PRECISION ZERO, ONE - PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) - COMPLEX*16 CZERO, CONE - PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ), CONE = ( 1.0D0, 0.0D0 ) ) -* .. -* .. Local Scalars .. - COMPLEX*16 CTEMP - DOUBLE PRECISION AAPP, AAQQ, AATMAX, AATMIN, BIG, BIG1, - $ COND_OK, CONDR1, CONDR2, ENTRA, ENTRAT, EPSLN, - $ MAXPRJ, SCALEM, SCONDA, SFMIN, SMALL, TEMP1, - $ USCAL1, USCAL2, XSC - INTEGER IERR, N1, NR, NUMRANK, p, q, WARNING - LOGICAL ALMORT, DEFR, ERREST, GOSCAL, JRACC, KILL, LQUERY, - $ LSVEC, L2ABER, L2KILL, L2PERT, L2RANK, L2TRAN, NOSCAL, - $ ROWPIV, RSVEC, TRANSP -* - INTEGER OPTWRK, MINWRK, MINRWRK, MINIWRK - INTEGER LWCON, LWLQF, LWQP3, LWQRF, LWUNMLQ, LWUNMQR, LWUNMQRM, - $ LWSVDJ, LWSVDJV, LRWQP3, LRWCON, LRWSVDJ, IWOFF - INTEGER LWRK_ZGELQF, LWRK_ZGEQP3, LWRK_ZGEQP3N, LWRK_ZGEQRF, - $ LWRK_ZGESVJ, LWRK_ZGESVJV, LWRK_ZGESVJU, LWRK_ZUNMLQ, - $ LWRK_ZUNMQR, LWRK_ZUNMQRM -* .. -* .. Local Arrays - COMPLEX*16 CDUMMY(1) - DOUBLE PRECISION RDUMMY(1) -* -* .. Intrinsic Functions .. - INTRINSIC ABS, DCMPLX, CONJG, DLOG, MAX, MIN, DBLE, NINT, SQRT -* .. -* .. External Functions .. - DOUBLE PRECISION DLAMCH, DZNRM2 - INTEGER IDAMAX, IZAMAX - LOGICAL LSAME - EXTERNAL IDAMAX, IZAMAX, LSAME, DLAMCH, DZNRM2 -* .. -* .. External Subroutines .. - EXTERNAL DLASSQ, ZCOPY, ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLAPMR, - $ ZLASCL, DLASCL, ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ, - $ ZUNMQR, ZPOCON, DSCAL, ZDSCAL, ZSWAP, ZTRSM, ZLACGV, - $ XERBLA -* - EXTERNAL ZGESVJ -* .. -* -* Test the input arguments -* - LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' ) - JRACC = LSAME( JOBV, 'J' ) - RSVEC = LSAME( JOBV, 'V' ) .OR. JRACC - ROWPIV = LSAME( JOBA, 'F' ) .OR. LSAME( JOBA, 'G' ) - L2RANK = LSAME( JOBA, 'R' ) - L2ABER = LSAME( JOBA, 'A' ) - ERREST = LSAME( JOBA, 'E' ) .OR. LSAME( JOBA, 'G' ) - L2TRAN = LSAME( JOBT, 'T' ) .AND. ( M .EQ. N ) - L2KILL = LSAME( JOBR, 'R' ) - DEFR = LSAME( JOBR, 'N' ) - L2PERT = LSAME( JOBP, 'P' ) -* - LQUERY = ( LWORK .EQ. -1 ) .OR. ( LRWORK .EQ. -1 ) -* - IF ( .NOT.(ROWPIV .OR. L2RANK .OR. L2ABER .OR. - $ ERREST .OR. LSAME( JOBA, 'C' ) )) THEN - INFO = - 1 - ELSE IF ( .NOT.( LSVEC .OR. LSAME( JOBU, 'N' ) .OR. - $ ( LSAME( JOBU, 'W' ) .AND. RSVEC .AND. L2TRAN ) ) ) THEN - INFO = - 2 - ELSE IF ( .NOT.( RSVEC .OR. LSAME( JOBV, 'N' ) .OR. - $ ( LSAME( JOBV, 'W' ) .AND. LSVEC .AND. L2TRAN ) ) ) THEN - INFO = - 3 - ELSE IF ( .NOT. ( L2KILL .OR. DEFR ) ) THEN - INFO = - 4 - ELSE IF ( .NOT. ( LSAME(JOBT,'T') .OR. LSAME(JOBT,'N') ) ) THEN - INFO = - 5 - ELSE IF ( .NOT. ( L2PERT .OR. LSAME( JOBP, 'N' ) ) ) THEN - INFO = - 6 - ELSE IF ( M .LT. 0 ) THEN - INFO = - 7 - ELSE IF ( ( N .LT. 0 ) .OR. ( N .GT. M ) ) THEN - INFO = - 8 - ELSE IF ( LDA .LT. M ) THEN - INFO = - 10 - ELSE IF ( LSVEC .AND. ( LDU .LT. M ) ) THEN - INFO = - 13 - ELSE IF ( RSVEC .AND. ( LDV .LT. N ) ) THEN - INFO = - 15 - ELSE -* #:) - INFO = 0 - END IF -* - IF ( INFO .EQ. 0 ) THEN -* .. compute the minimal and the optimal workspace lengths -* [[The expressions for computing the minimal and the optimal -* values of LCWORK, LRWORK are written with a lot of redundancy and -* can be simplified. However, this verbose form is useful for -* maintenance and modifications of the code.]] -* -* .. minimal workspace length for ZGEQP3 of an M x N matrix, -* ZGEQRF of an N x N matrix, ZGELQF of an N x N matrix, -* ZUNMLQ for computing N x N matrix, ZUNMQR for computing N x N -* matrix, ZUNMQR for computing M x N matrix, respectively. - LWQP3 = N+1 - LWQRF = MAX( 1, N ) - LWLQF = MAX( 1, N ) - LWUNMLQ = MAX( 1, N ) - LWUNMQR = MAX( 1, N ) - LWUNMQRM = MAX( 1, M ) -* .. minimal workspace length for ZPOCON of an N x N matrix - LWCON = 2 * N -* .. minimal workspace length for ZGESVJ of an N x N matrix, -* without and with explicit accumulation of Jacobi rotations - LWSVDJ = MAX( 2 * N, 1 ) - LWSVDJV = MAX( 2 * N, 1 ) -* .. minimal REAL workspace length for ZGEQP3, ZPOCON, ZGESVJ - LRWQP3 = 2 * N - LRWCON = N - LRWSVDJ = N - IF ( LQUERY ) THEN - CALL ZGEQP3( M, N, A, LDA, IWORK, CDUMMY, CDUMMY, -1, - $ RDUMMY, IERR ) - LWRK_ZGEQP3 = INT( CDUMMY(1) ) - CALL ZGEQRF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR ) - LWRK_ZGEQRF = INT( CDUMMY(1) ) - CALL ZGELQF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR ) - LWRK_ZGELQF = INT( CDUMMY(1) ) - END IF - MINWRK = 2 - OPTWRK = 2 - MINIWRK = N - IF ( .NOT. (LSVEC .OR. RSVEC ) ) THEN -* .. minimal and optimal sizes of the complex workspace if -* only the singular values are requested - IF ( ERREST ) THEN - MINWRK = MAX( N+LWQP3, N**2+LWCON, N+LWQRF, LWSVDJ ) - ELSE - MINWRK = MAX( N+LWQP3, N+LWQRF, LWSVDJ ) - END IF - IF ( LQUERY ) THEN - CALL ZGESVJ( 'L', 'N', 'N', N, N, A, LDA, SVA, N, V, - $ LDV, CDUMMY, -1, RDUMMY, -1, IERR ) - LWRK_ZGESVJ = INT( CDUMMY(1) ) - IF ( ERREST ) THEN - OPTWRK = MAX( N+LWRK_ZGEQP3, N**2+LWCON, - $ N+LWRK_ZGEQRF, LWRK_ZGESVJ ) - ELSE - OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWRK_ZGEQRF, - $ LWRK_ZGESVJ ) - END IF - END IF - IF ( L2TRAN .OR. ROWPIV ) THEN - IF ( ERREST ) THEN - MINRWRK = MAX( 7, 2*M, LRWQP3, LRWCON, LRWSVDJ ) - ELSE - MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ ) - END IF - ELSE - IF ( ERREST ) THEN - MINRWRK = MAX( 7, LRWQP3, LRWCON, LRWSVDJ ) - ELSE - MINRWRK = MAX( 7, LRWQP3, LRWSVDJ ) - END IF - END IF - IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M - ELSE IF ( RSVEC .AND. (.NOT.LSVEC) ) THEN -* .. minimal and optimal sizes of the complex workspace if the -* singular values and the right singular vectors are requested - IF ( ERREST ) THEN - MINWRK = MAX( N+LWQP3, LWCON, LWSVDJ, N+LWLQF, - $ 2*N+LWQRF, N+LWSVDJ, N+LWUNMLQ ) - ELSE - MINWRK = MAX( N+LWQP3, LWSVDJ, N+LWLQF, 2*N+LWQRF, - $ N+LWSVDJ, N+LWUNMLQ ) - END IF - IF ( LQUERY ) THEN - CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A, - $ LDA, CDUMMY, -1, RDUMMY, -1, IERR ) - LWRK_ZGESVJ = INT( CDUMMY(1) ) - CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY, - $ V, LDV, CDUMMY, -1, IERR ) - LWRK_ZUNMLQ = INT( CDUMMY(1) ) - IF ( ERREST ) THEN - OPTWRK = MAX( N+LWRK_ZGEQP3, LWCON, LWRK_ZGESVJ, - $ N+LWRK_ZGELQF, 2*N+LWRK_ZGEQRF, - $ N+LWRK_ZGESVJ, N+LWRK_ZUNMLQ ) - ELSE - OPTWRK = MAX( N+LWRK_ZGEQP3, LWRK_ZGESVJ,N+LWRK_ZGELQF, - $ 2*N+LWRK_ZGEQRF, N+LWRK_ZGESVJ, - $ N+LWRK_ZUNMLQ ) - END IF - END IF - IF ( L2TRAN .OR. ROWPIV ) THEN - IF ( ERREST ) THEN - MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ, LRWCON ) - ELSE - MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ ) - END IF - ELSE - IF ( ERREST ) THEN - MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON ) - ELSE - MINRWRK = MAX( 7, LRWQP3, LRWSVDJ ) - END IF - END IF - IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M - ELSE IF ( LSVEC .AND. (.NOT.RSVEC) ) THEN -* .. minimal and optimal sizes of the complex workspace if the -* singular values and the left singular vectors are requested - IF ( ERREST ) THEN - MINWRK = N + MAX( LWQP3,LWCON,N+LWQRF,LWSVDJ,LWUNMQRM ) - ELSE - MINWRK = N + MAX( LWQP3, N+LWQRF, LWSVDJ, LWUNMQRM ) - END IF - IF ( LQUERY ) THEN - CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A, - $ LDA, CDUMMY, -1, RDUMMY, -1, IERR ) - LWRK_ZGESVJ = INT( CDUMMY(1) ) - CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U, - $ LDU, CDUMMY, -1, IERR ) - LWRK_ZUNMQRM = INT( CDUMMY(1) ) - IF ( ERREST ) THEN - OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, N+LWRK_ZGEQRF, - $ LWRK_ZGESVJ, LWRK_ZUNMQRM ) - ELSE - OPTWRK = N + MAX( LWRK_ZGEQP3, N+LWRK_ZGEQRF, - $ LWRK_ZGESVJ, LWRK_ZUNMQRM ) - END IF - END IF - IF ( L2TRAN .OR. ROWPIV ) THEN - IF ( ERREST ) THEN - MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ, LRWCON ) - ELSE - MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ ) - END IF - ELSE - IF ( ERREST ) THEN - MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON ) - ELSE - MINRWRK = MAX( 7, LRWQP3, LRWSVDJ ) - END IF - END IF - IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M - ELSE -* .. minimal and optimal sizes of the complex workspace if the -* full SVD is requested - IF ( .NOT. JRACC ) THEN - IF ( ERREST ) THEN - MINWRK = MAX( N+LWQP3, N+LWCON, 2*N+N**2+LWCON, - $ 2*N+LWQRF, 2*N+LWQP3, - $ 2*N+N**2+N+LWLQF, 2*N+N**2+N+N**2+LWCON, - $ 2*N+N**2+N+LWSVDJ, 2*N+N**2+N+LWSVDJV, - $ 2*N+N**2+N+LWUNMQR,2*N+N**2+N+LWUNMLQ, - $ N+N**2+LWSVDJ, N+LWUNMQRM ) - ELSE - MINWRK = MAX( N+LWQP3, 2*N+N**2+LWCON, - $ 2*N+LWQRF, 2*N+LWQP3, - $ 2*N+N**2+N+LWLQF, 2*N+N**2+N+N**2+LWCON, - $ 2*N+N**2+N+LWSVDJ, 2*N+N**2+N+LWSVDJV, - $ 2*N+N**2+N+LWUNMQR,2*N+N**2+N+LWUNMLQ, - $ N+N**2+LWSVDJ, N+LWUNMQRM ) - END IF - MINIWRK = MINIWRK + N - IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M - ELSE - IF ( ERREST ) THEN - MINWRK = MAX( N+LWQP3, N+LWCON, 2*N+LWQRF, - $ 2*N+N**2+LWSVDJV, 2*N+N**2+N+LWUNMQR, - $ N+LWUNMQRM ) - ELSE - MINWRK = MAX( N+LWQP3, 2*N+LWQRF, - $ 2*N+N**2+LWSVDJV, 2*N+N**2+N+LWUNMQR, - $ N+LWUNMQRM ) - END IF - IF ( ROWPIV .OR. L2TRAN ) MINIWRK = MINIWRK + M - END IF - IF ( LQUERY ) THEN - CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U, - $ LDU, CDUMMY, -1, IERR ) - LWRK_ZUNMQRM = INT( CDUMMY(1) ) - CALL ZUNMQR( 'L', 'N', N, N, N, A, LDA, CDUMMY, U, - $ LDU, CDUMMY, -1, IERR ) - LWRK_ZUNMQR = INT( CDUMMY(1) ) - IF ( .NOT. JRACC ) THEN - CALL ZGEQP3( N,N, A, LDA, IWORK, CDUMMY,CDUMMY, -1, - $ RDUMMY, IERR ) - LWRK_ZGEQP3N = INT( CDUMMY(1) ) - CALL ZGESVJ( 'L', 'U', 'N', N, N, U, LDU, SVA, - $ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR ) - LWRK_ZGESVJ = INT( CDUMMY(1) ) - CALL ZGESVJ( 'U', 'U', 'N', N, N, U, LDU, SVA, - $ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR ) - LWRK_ZGESVJU = INT( CDUMMY(1) ) - CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA, - $ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR ) - LWRK_ZGESVJV = INT( CDUMMY(1) ) - CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY, - $ V, LDV, CDUMMY, -1, IERR ) - LWRK_ZUNMLQ = INT( CDUMMY(1) ) - IF ( ERREST ) THEN - OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON, - $ 2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF, - $ 2*N+LWRK_ZGEQP3N, - $ 2*N+N**2+N+LWRK_ZGELQF, - $ 2*N+N**2+N+N**2+LWCON, - $ 2*N+N**2+N+LWRK_ZGESVJ, - $ 2*N+N**2+N+LWRK_ZGESVJV, - $ 2*N+N**2+N+LWRK_ZUNMQR, - $ 2*N+N**2+N+LWRK_ZUNMLQ, - $ N+N**2+LWRK_ZGESVJU, - $ N+LWRK_ZUNMQRM ) - ELSE - OPTWRK = MAX( N+LWRK_ZGEQP3, - $ 2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF, - $ 2*N+LWRK_ZGEQP3N, - $ 2*N+N**2+N+LWRK_ZGELQF, - $ 2*N+N**2+N+N**2+LWCON, - $ 2*N+N**2+N+LWRK_ZGESVJ, - $ 2*N+N**2+N+LWRK_ZGESVJV, - $ 2*N+N**2+N+LWRK_ZUNMQR, - $ 2*N+N**2+N+LWRK_ZUNMLQ, - $ N+N**2+LWRK_ZGESVJU, - $ N+LWRK_ZUNMQRM ) - END IF - ELSE - CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA, - $ N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR ) - LWRK_ZGESVJV = INT( CDUMMY(1) ) - CALL ZUNMQR( 'L', 'N', N, N, N, CDUMMY, N, CDUMMY, - $ V, LDV, CDUMMY, -1, IERR ) - LWRK_ZUNMQR = INT( CDUMMY(1) ) - CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U, - $ LDU, CDUMMY, -1, IERR ) - LWRK_ZUNMQRM = INT( CDUMMY(1) ) - IF ( ERREST ) THEN - OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON, - $ 2*N+LWRK_ZGEQRF, 2*N+N**2, - $ 2*N+N**2+LWRK_ZGESVJV, - $ 2*N+N**2+N+LWRK_ZUNMQR,N+LWRK_ZUNMQRM ) - ELSE - OPTWRK = MAX( N+LWRK_ZGEQP3, 2*N+LWRK_ZGEQRF, - $ 2*N+N**2, 2*N+N**2+LWRK_ZGESVJV, - $ 2*N+N**2+N+LWRK_ZUNMQR, - $ N+LWRK_ZUNMQRM ) - END IF - END IF - END IF - IF ( L2TRAN .OR. ROWPIV ) THEN - MINRWRK = MAX( 7, 2*M, LRWQP3, LRWSVDJ, LRWCON ) - ELSE - MINRWRK = MAX( 7, LRWQP3, LRWSVDJ, LRWCON ) - END IF - END IF - MINWRK = MAX( 2, MINWRK ) - OPTWRK = MAX( MINWRK, OPTWRK ) - IF ( LWORK .LT. MINWRK .AND. (.NOT.LQUERY) ) INFO = - 17 - IF ( LRWORK .LT. MINRWRK .AND. (.NOT.LQUERY) ) INFO = - 19 - END IF -* - IF ( INFO .NE. 0 ) THEN -* #:( - CALL XERBLA( 'ZGEJSV', - INFO ) - RETURN - ELSE IF ( LQUERY ) THEN - CWORK(1) = OPTWRK - CWORK(2) = MINWRK - RWORK(1) = MINRWRK - IWORK(1) = MAX( 4, MINIWRK ) - RETURN - END IF -* -* Quick return for void matrix (Y3K safe) -* #:) - IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN - IWORK(1:4) = 0 - RWORK(1:7) = 0 - RETURN - ENDIF -* -* Determine whether the matrix U should be M x N or M x M -* - IF ( LSVEC ) THEN - N1 = N - IF ( LSAME( JOBU, 'F' ) ) N1 = M - END IF -* -* Set numerical parameters -* -*! NOTE: Make sure DLAMCH() does not fail on the target architecture. -* - EPSLN = DLAMCH('Epsilon') - SFMIN = DLAMCH('SafeMinimum') - SMALL = SFMIN / EPSLN - BIG = DLAMCH('O') -* BIG = ONE / SFMIN -* -* Initialize SVA(1:N) = diag( ||A e_i||_2 )_1^N -* -*(!) If necessary, scale SVA() to protect the largest norm from -* overflow. It is possible that this scaling pushes the smallest -* column norm left from the underflow threshold (extreme case). -* - SCALEM = ONE / SQRT(DBLE(M)*DBLE(N)) - NOSCAL = .TRUE. - GOSCAL = .TRUE. - DO 1874 p = 1, N - AAPP = ZERO - AAQQ = ONE - CALL ZLASSQ( M, A(1,p), 1, AAPP, AAQQ ) - IF ( AAPP .GT. BIG ) THEN - INFO = - 9 - CALL XERBLA( 'ZGEJSV', -INFO ) - RETURN - END IF - AAQQ = SQRT(AAQQ) - IF ( ( AAPP .LT. (BIG / AAQQ) ) .AND. NOSCAL ) THEN - SVA(p) = AAPP * AAQQ - ELSE - NOSCAL = .FALSE. - SVA(p) = AAPP * ( AAQQ * SCALEM ) - IF ( GOSCAL ) THEN - GOSCAL = .FALSE. - CALL DSCAL( p-1, SCALEM, SVA, 1 ) - END IF - END IF - 1874 CONTINUE -* - IF ( NOSCAL ) SCALEM = ONE -* - AAPP = ZERO - AAQQ = BIG - DO 4781 p = 1, N - AAPP = MAX( AAPP, SVA(p) ) - IF ( SVA(p) .NE. ZERO ) AAQQ = MIN( AAQQ, SVA(p) ) - 4781 CONTINUE -* -* Quick return for zero M x N matrix -* #:) - IF ( AAPP .EQ. ZERO ) THEN - IF ( LSVEC ) CALL ZLASET( 'G', M, N1, CZERO, CONE, U, LDU ) - IF ( RSVEC ) CALL ZLASET( 'G', N, N, CZERO, CONE, V, LDV ) - RWORK(1) = ONE - RWORK(2) = ONE - IF ( ERREST ) RWORK(3) = ONE - IF ( LSVEC .AND. RSVEC ) THEN - RWORK(4) = ONE - RWORK(5) = ONE - END IF - IF ( L2TRAN ) THEN - RWORK(6) = ZERO - RWORK(7) = ZERO - END IF - IWORK(1) = 0 - IWORK(2) = 0 - IWORK(3) = 0 - IWORK(4) = -1 - RETURN - END IF -* -* Issue warning if denormalized column norms detected. Override the -* high relative accuracy request. Issue licence to kill nonzero columns -* (set them to zero) whose norm is less than sigma_max / BIG (roughly). -* #:( - WARNING = 0 - IF ( AAQQ .LE. SFMIN ) THEN - L2RANK = .TRUE. - L2KILL = .TRUE. - WARNING = 1 - END IF -* -* Quick return for one-column matrix -* #:) - IF ( N .EQ. 1 ) THEN -* - IF ( LSVEC ) THEN - CALL ZLASCL( 'G',0,0,SVA(1),SCALEM, M,1,A(1,1),LDA,IERR ) - CALL ZLACPY( 'A', M, 1, A, LDA, U, LDU ) -* computing all M left singular vectors of the M x 1 matrix - IF ( N1 .NE. N ) THEN - CALL ZGEQRF( M, N, U,LDU, CWORK, CWORK(N+1),LWORK-N,IERR ) - CALL ZUNGQR( M,N1,1, U,LDU,CWORK,CWORK(N+1),LWORK-N,IERR ) - CALL ZCOPY( M, A(1,1), 1, U(1,1), 1 ) - END IF - END IF - IF ( RSVEC ) THEN - V(1,1) = CONE - END IF - IF ( SVA(1) .LT. (BIG*SCALEM) ) THEN - SVA(1) = SVA(1) / SCALEM - SCALEM = ONE - END IF - RWORK(1) = ONE / SCALEM - RWORK(2) = ONE - IF ( SVA(1) .NE. ZERO ) THEN - IWORK(1) = 1 - IF ( ( SVA(1) / SCALEM) .GE. SFMIN ) THEN - IWORK(2) = 1 - ELSE - IWORK(2) = 0 - END IF - ELSE - IWORK(1) = 0 - IWORK(2) = 0 - END IF - IWORK(3) = 0 - IWORK(4) = -1 - IF ( ERREST ) RWORK(3) = ONE - IF ( LSVEC .AND. RSVEC ) THEN - RWORK(4) = ONE - RWORK(5) = ONE - END IF - IF ( L2TRAN ) THEN - RWORK(6) = ZERO - RWORK(7) = ZERO - END IF - RETURN -* - END IF -* - TRANSP = .FALSE. -* - AATMAX = -ONE - AATMIN = BIG - IF ( ROWPIV .OR. L2TRAN ) THEN -* -* Compute the row norms, needed to determine row pivoting sequence -* (in the case of heavily row weighted A, row pivoting is strongly -* advised) and to collect information needed to compare the -* structures of A * A^* and A^* * A (in the case L2TRAN.EQ..TRUE.). -* - IF ( L2TRAN ) THEN - DO 1950 p = 1, M - XSC = ZERO - TEMP1 = ONE - CALL ZLASSQ( N, A(p,1), LDA, XSC, TEMP1 ) -* ZLASSQ gets both the ell_2 and the ell_infinity norm -* in one pass through the vector - RWORK(M+p) = XSC * SCALEM - RWORK(p) = XSC * (SCALEM*SQRT(TEMP1)) - AATMAX = MAX( AATMAX, RWORK(p) ) - IF (RWORK(p) .NE. ZERO) - $ AATMIN = MIN(AATMIN,RWORK(p)) - 1950 CONTINUE - ELSE - DO 1904 p = 1, M - RWORK(M+p) = SCALEM*ABS( A(p,IZAMAX(N,A(p,1),LDA)) ) - AATMAX = MAX( AATMAX, RWORK(M+p) ) - AATMIN = MIN( AATMIN, RWORK(M+p) ) - 1904 CONTINUE - END IF -* - END IF -* -* For square matrix A try to determine whether A^* would be better -* input for the preconditioned Jacobi SVD, with faster convergence. -* The decision is based on an O(N) function of the vector of column -* and row norms of A, based on the Shannon entropy. This should give -* the right choice in most cases when the difference actually matters. -* It may fail and pick the slower converging side. -* - ENTRA = ZERO - ENTRAT = ZERO - IF ( L2TRAN ) THEN -* - XSC = ZERO - TEMP1 = ONE - CALL DLASSQ( N, SVA, 1, XSC, TEMP1 ) - TEMP1 = ONE / TEMP1 -* - ENTRA = ZERO - DO 1113 p = 1, N - BIG1 = ( ( SVA(p) / XSC )**2 ) * TEMP1 - IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1) - 1113 CONTINUE - ENTRA = - ENTRA / DLOG(DBLE(N)) -* -* Now, SVA().^2/Trace(A^* * A) is a point in the probability simplex. -* It is derived from the diagonal of A^* * A. Do the same with the -* diagonal of A * A^*, compute the entropy of the corresponding -* probability distribution. Note that A * A^* and A^* * A have the -* same trace. -* - ENTRAT = ZERO - DO 1114 p = 1, M - BIG1 = ( ( RWORK(p) / XSC )**2 ) * TEMP1 - IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1) - 1114 CONTINUE - ENTRAT = - ENTRAT / DLOG(DBLE(M)) -* -* Analyze the entropies and decide A or A^*. Smaller entropy -* usually means better input for the algorithm. -* - TRANSP = ( ENTRAT .LT. ENTRA ) -* -* If A^* is better than A, take the adjoint of A. This is allowed -* only for square matrices, M=N. - IF ( TRANSP ) THEN -* In an optimal implementation, this trivial transpose -* should be replaced with faster transpose. - DO 1115 p = 1, N - 1 - A(p,p) = CONJG(A(p,p)) - DO 1116 q = p + 1, N - CTEMP = CONJG(A(q,p)) - A(q,p) = CONJG(A(p,q)) - A(p,q) = CTEMP - 1116 CONTINUE - 1115 CONTINUE - A(N,N) = CONJG(A(N,N)) - DO 1117 p = 1, N - RWORK(M+p) = SVA(p) - SVA(p) = RWORK(p) -* previously computed row 2-norms are now column 2-norms -* of the transposed matrix - 1117 CONTINUE - TEMP1 = AAPP - AAPP = AATMAX - AATMAX = TEMP1 - TEMP1 = AAQQ - AAQQ = AATMIN - AATMIN = TEMP1 - KILL = LSVEC - LSVEC = RSVEC - RSVEC = KILL - IF ( LSVEC ) N1 = N -* - ROWPIV = .TRUE. - END IF -* - END IF -* END IF L2TRAN -* -* Scale the matrix so that its maximal singular value remains less -* than SQRT(BIG) -- the matrix is scaled so that its maximal column -* has Euclidean norm equal to SQRT(BIG/N). The only reason to keep -* SQRT(BIG) instead of BIG is the fact that ZGEJSV uses LAPACK and -* BLAS routines that, in some implementations, are not capable of -* working in the full interval [SFMIN,BIG] and that they may provoke -* overflows in the intermediate results. If the singular values spread -* from SFMIN to BIG, then ZGESVJ will compute them. So, in that case, -* one should use ZGESVJ instead of ZGEJSV. -* >> change in the April 2016 update: allow bigger range, i.e. the -* largest column is allowed up to BIG/N and ZGESVJ will do the rest. - BIG1 = SQRT( BIG ) - TEMP1 = SQRT( BIG / DBLE(N) ) -* TEMP1 = BIG/DBLE(N) -* - CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR ) - IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN - AAQQ = ( AAQQ / AAPP ) * TEMP1 - ELSE - AAQQ = ( AAQQ * TEMP1 ) / AAPP - END IF - TEMP1 = TEMP1 * SCALEM - CALL ZLASCL( 'G', 0, 0, AAPP, TEMP1, M, N, A, LDA, IERR ) -* -* To undo scaling at the end of this procedure, multiply the -* computed singular values with USCAL2 / USCAL1. -* - USCAL1 = TEMP1 - USCAL2 = AAPP -* - IF ( L2KILL ) THEN -* L2KILL enforces computation of nonzero singular values in -* the restricted range of condition number of the initial A, -* sigma_max(A) / sigma_min(A) approx. SQRT(BIG)/SQRT(SFMIN). - XSC = SQRT( SFMIN ) - ELSE - XSC = SMALL -* -* Now, if the condition number of A is too big, -* sigma_max(A) / sigma_min(A) .GT. SQRT(BIG/N) * EPSLN / SFMIN, -* as a precaution measure, the full SVD is computed using ZGESVJ -* with accumulated Jacobi rotations. This provides numerically -* more robust computation, at the cost of slightly increased run -* time. Depending on the concrete implementation of BLAS and LAPACK -* (i.e. how they behave in presence of extreme ill-conditioning) the -* implementor may decide to remove this switch. - IF ( ( AAQQ.LT.SQRT(SFMIN) ) .AND. LSVEC .AND. RSVEC ) THEN - JRACC = .TRUE. - END IF -* - END IF - IF ( AAQQ .LT. XSC ) THEN - DO 700 p = 1, N - IF ( SVA(p) .LT. XSC ) THEN - CALL ZLASET( 'A', M, 1, CZERO, CZERO, A(1,p), LDA ) - SVA(p) = ZERO - END IF - 700 CONTINUE - END IF -* -* Preconditioning using QR factorization with pivoting -* - IF ( ROWPIV ) THEN -* Optional row permutation (Bjoerck row pivoting): -* A result by Cox and Higham shows that the Bjoerck's -* row pivoting combined with standard column pivoting -* has similar effect as Powell-Reid complete pivoting. -* The ell-infinity norms of A are made nonincreasing. - IF ( ( LSVEC .AND. RSVEC ) .AND. .NOT.( JRACC ) ) THEN - IWOFF = 2*N - ELSE - IWOFF = N - END IF - DO 1952 p = 1, M - 1 - q = IDAMAX( M-p+1, RWORK(M+p), 1 ) + p - 1 - IWORK(IWOFF+p) = q - IF ( p .NE. q ) THEN - TEMP1 = RWORK(M+p) - RWORK(M+p) = RWORK(M+q) - RWORK(M+q) = TEMP1 - END IF - 1952 CONTINUE - CALL ZLASWP( N, A, LDA, 1, M-1, IWORK(IWOFF+1), 1 ) - END IF -* -* End of the preparation phase (scaling, optional sorting and -* transposing, optional flushing of small columns). -* -* Preconditioning -* -* If the full SVD is needed, the right singular vectors are computed -* from a matrix equation, and for that we need theoretical analysis -* of the Businger-Golub pivoting. So we use ZGEQP3 as the first RR QRF. -* In all other cases the first RR QRF can be chosen by other criteria -* (eg speed by replacing global with restricted window pivoting, such -* as in xGEQPX from TOMS # 782). Good results will be obtained using -* xGEQPX with properly (!) chosen numerical parameters. -* Any improvement of ZGEQP3 improves overall performance of ZGEJSV. -* -* A * P1 = Q1 * [ R1^* 0]^*: - DO 1963 p = 1, N -* .. all columns are free columns - IWORK(p) = 0 - 1963 CONTINUE - CALL ZGEQP3( M, N, A, LDA, IWORK, CWORK, CWORK(N+1), LWORK-N, - $ RWORK, IERR ) -* -* The upper triangular matrix R1 from the first QRF is inspected for -* rank deficiency and possibilities for deflation, or possible -* ill-conditioning. Depending on the user specified flag L2RANK, -* the procedure explores possibilities to reduce the numerical -* rank by inspecting the computed upper triangular factor. If -* L2RANK or L2ABER are up, then ZGEJSV will compute the SVD of -* A + dA, where ||dA|| <= f(M,N)*EPSLN. -* - NR = 1 - IF ( L2ABER ) THEN -* Standard absolute error bound suffices. All sigma_i with -* sigma_i < N*EPSLN*||A|| are flushed to zero. This is an -* aggressive enforcement of lower numerical rank by introducing a -* backward error of the order of N*EPSLN*||A||. - TEMP1 = SQRT(DBLE(N))*EPSLN - DO 3001 p = 2, N - IF ( ABS(A(p,p)) .GE. (TEMP1*ABS(A(1,1))) ) THEN - NR = NR + 1 - ELSE - GO TO 3002 - END IF - 3001 CONTINUE - 3002 CONTINUE - ELSE IF ( L2RANK ) THEN -* .. similarly as above, only slightly more gentle (less aggressive). -* Sudden drop on the diagonal of R1 is used as the criterion for -* close-to-rank-deficient. - TEMP1 = SQRT(SFMIN) - DO 3401 p = 2, N - IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR. - $ ( ABS(A(p,p)) .LT. SMALL ) .OR. - $ ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3402 - NR = NR + 1 - 3401 CONTINUE - 3402 CONTINUE -* - ELSE -* The goal is high relative accuracy. However, if the matrix -* has high scaled condition number the relative accuracy is in -* general not feasible. Later on, a condition number estimator -* will be deployed to estimate the scaled condition number. -* Here we just remove the underflowed part of the triangular -* factor. This prevents the situation in which the code is -* working hard to get the accuracy not warranted by the data. - TEMP1 = SQRT(SFMIN) - DO 3301 p = 2, N - IF ( ( ABS(A(p,p)) .LT. SMALL ) .OR. - $ ( L2KILL .AND. (ABS(A(p,p)) .LT. TEMP1) ) ) GO TO 3302 - NR = NR + 1 - 3301 CONTINUE - 3302 CONTINUE -* - END IF -* - ALMORT = .FALSE. - IF ( NR .EQ. N ) THEN - MAXPRJ = ONE - DO 3051 p = 2, N - TEMP1 = ABS(A(p,p)) / SVA(IWORK(p)) - MAXPRJ = MIN( MAXPRJ, TEMP1 ) - 3051 CONTINUE - IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE. - END IF -* -* - SCONDA = - ONE - CONDR1 = - ONE - CONDR2 = - ONE -* - IF ( ERREST ) THEN - IF ( N .EQ. NR ) THEN - IF ( RSVEC ) THEN -* .. V is available as workspace - CALL ZLACPY( 'U', N, N, A, LDA, V, LDV ) - DO 3053 p = 1, N - TEMP1 = SVA(IWORK(p)) - CALL ZDSCAL( p, ONE/TEMP1, V(1,p), 1 ) - 3053 CONTINUE - IF ( LSVEC )THEN - CALL ZPOCON( 'U', N, V, LDV, ONE, TEMP1, - $ CWORK(N+1), RWORK, IERR ) - ELSE - CALL ZPOCON( 'U', N, V, LDV, ONE, TEMP1, - $ CWORK, RWORK, IERR ) - END IF -* - ELSE IF ( LSVEC ) THEN -* .. U is available as workspace - CALL ZLACPY( 'U', N, N, A, LDA, U, LDU ) - DO 3054 p = 1, N - TEMP1 = SVA(IWORK(p)) - CALL ZDSCAL( p, ONE/TEMP1, U(1,p), 1 ) - 3054 CONTINUE - CALL ZPOCON( 'U', N, U, LDU, ONE, TEMP1, - $ CWORK(N+1), RWORK, IERR ) - ELSE - CALL ZLACPY( 'U', N, N, A, LDA, CWORK, N ) -*[] CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N ) -* Change: here index shifted by N to the left, CWORK(1:N) -* not needed for SIGMA only computation - DO 3052 p = 1, N - TEMP1 = SVA(IWORK(p)) -*[] CALL ZDSCAL( p, ONE/TEMP1, CWORK(N+(p-1)*N+1), 1 ) - CALL ZDSCAL( p, ONE/TEMP1, CWORK((p-1)*N+1), 1 ) - 3052 CONTINUE -* .. the columns of R are scaled to have unit Euclidean lengths. -*[] CALL ZPOCON( 'U', N, CWORK(N+1), N, ONE, TEMP1, -*[] $ CWORK(N+N*N+1), RWORK, IERR ) - CALL ZPOCON( 'U', N, CWORK, N, ONE, TEMP1, - $ CWORK(N*N+1), RWORK, IERR ) -* - END IF - IF ( TEMP1 .NE. ZERO ) THEN - SCONDA = ONE / SQRT(TEMP1) - ELSE - SCONDA = - ONE - END IF -* SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1). -* N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA - ELSE - SCONDA = - ONE - END IF - END IF -* - L2PERT = L2PERT .AND. ( ABS( A(1,1)/A(NR,NR) ) .GT. SQRT(BIG1) ) -* If there is no violent scaling, artificial perturbation is not needed. -* -* Phase 3: -* - IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN -* -* Singular Values only -* -* .. transpose A(1:NR,1:N) - DO 1946 p = 1, MIN( N-1, NR ) - CALL ZCOPY( N-p, A(p,p+1), LDA, A(p+1,p), 1 ) - CALL ZLACGV( N-p+1, A(p,p), 1 ) - 1946 CONTINUE - IF ( NR .EQ. N ) A(N,N) = CONJG(A(N,N)) -* -* The following two DO-loops introduce small relative perturbation -* into the strict upper triangle of the lower triangular matrix. -* Small entries below the main diagonal are also changed. -* This modification is useful if the computing environment does not -* provide/allow FLUSH TO ZERO underflow, for it prevents many -* annoying denormalized numbers in case of strongly scaled matrices. -* The perturbation is structured so that it does not introduce any -* new perturbation of the singular values, and it does not destroy -* the job done by the preconditioner. -* The licence for this perturbation is in the variable L2PERT, which -* should be .FALSE. if FLUSH TO ZERO underflow is active. -* - IF ( .NOT. ALMORT ) THEN -* - IF ( L2PERT ) THEN -* XSC = SQRT(SMALL) - XSC = EPSLN / DBLE(N) - DO 4947 q = 1, NR - CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO) - DO 4949 p = 1, N - IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) ) - $ .OR. ( p .LT. q ) ) -* $ A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) ) - $ A(p,q) = CTEMP - 4949 CONTINUE - 4947 CONTINUE - ELSE - CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, A(1,2),LDA ) - END IF -* -* .. second preconditioning using the QR factorization -* - CALL ZGEQRF( N,NR, A,LDA, CWORK, CWORK(N+1),LWORK-N, IERR ) -* -* .. and transpose upper to lower triangular - DO 1948 p = 1, NR - 1 - CALL ZCOPY( NR-p, A(p,p+1), LDA, A(p+1,p), 1 ) - CALL ZLACGV( NR-p+1, A(p,p), 1 ) - 1948 CONTINUE -* - END IF -* -* Row-cyclic Jacobi SVD algorithm with column pivoting -* -* .. again some perturbation (a "background noise") is added -* to drown denormals - IF ( L2PERT ) THEN -* XSC = SQRT(SMALL) - XSC = EPSLN / DBLE(N) - DO 1947 q = 1, NR - CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO) - DO 1949 p = 1, NR - IF ( ( (p.GT.q) .AND. (ABS(A(p,q)).LE.TEMP1) ) - $ .OR. ( p .LT. q ) ) -* $ A(p,q) = TEMP1 * ( A(p,q) / ABS(A(p,q)) ) - $ A(p,q) = CTEMP - 1949 CONTINUE - 1947 CONTINUE - ELSE - CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, A(1,2), LDA ) - END IF -* -* .. and one-sided Jacobi rotations are started on a lower -* triangular matrix (plus perturbation which is ignored in -* the part which destroys triangular form (confusing?!)) -* - CALL ZGESVJ( 'L', 'N', 'N', NR, NR, A, LDA, SVA, - $ N, V, LDV, CWORK, LWORK, RWORK, LRWORK, INFO ) -* - SCALEM = RWORK(1) - NUMRANK = NINT(RWORK(2)) -* -* - ELSE IF ( ( RSVEC .AND. ( .NOT. LSVEC ) .AND. ( .NOT. JRACC ) ) - $ .OR. - $ ( JRACC .AND. ( .NOT. LSVEC ) .AND. ( NR .NE. N ) ) ) THEN -* -* -> Singular Values and Right Singular Vectors <- -* - IF ( ALMORT ) THEN -* -* .. in this case NR equals N - DO 1998 p = 1, NR - CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) - CALL ZLACGV( N-p+1, V(p,p), 1 ) - 1998 CONTINUE - CALL ZLASET( 'U', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV ) -* - CALL ZGESVJ( 'L','U','N', N, NR, V, LDV, SVA, NR, A, LDA, - $ CWORK, LWORK, RWORK, LRWORK, INFO ) - SCALEM = RWORK(1) - NUMRANK = NINT(RWORK(2)) - - ELSE -* -* .. two more QR factorizations ( one QRF is not enough, two require -* accumulated product of Jacobi rotations, three are perfect ) -* - CALL ZLASET( 'L', NR-1,NR-1, CZERO, CZERO, A(2,1), LDA ) - CALL ZGELQF( NR,N, A, LDA, CWORK, CWORK(N+1), LWORK-N, IERR) - CALL ZLACPY( 'L', NR, NR, A, LDA, V, LDV ) - CALL ZLASET( 'U', NR-1,NR-1, CZERO, CZERO, V(1,2), LDV ) - CALL ZGEQRF( NR, NR, V, LDV, CWORK(N+1), CWORK(2*N+1), - $ LWORK-2*N, IERR ) - DO 8998 p = 1, NR - CALL ZCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 ) - CALL ZLACGV( NR-p+1, V(p,p), 1 ) - 8998 CONTINUE - CALL ZLASET('U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV) -* - CALL ZGESVJ( 'L', 'U','N', NR, NR, V,LDV, SVA, NR, U, - $ LDU, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO ) - SCALEM = RWORK(1) - NUMRANK = NINT(RWORK(2)) - IF ( NR .LT. N ) THEN - CALL ZLASET( 'A',N-NR, NR, CZERO,CZERO, V(NR+1,1), LDV ) - CALL ZLASET( 'A',NR, N-NR, CZERO,CZERO, V(1,NR+1), LDV ) - CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE, V(NR+1,NR+1),LDV ) - END IF -* - CALL ZUNMLQ( 'L', 'C', N, N, NR, A, LDA, CWORK, - $ V, LDV, CWORK(N+1), LWORK-N, IERR ) -* - END IF -* .. permute the rows of V -* DO 8991 p = 1, N -* CALL ZCOPY( N, V(p,1), LDV, A(IWORK(p),1), LDA ) -* 8991 CONTINUE -* CALL ZLACPY( 'All', N, N, A, LDA, V, LDV ) - CALL ZLAPMR( .FALSE., N, N, V, LDV, IWORK ) -* - IF ( TRANSP ) THEN - CALL ZLACPY( 'A', N, N, V, LDV, U, LDU ) - END IF -* - ELSE IF ( JRACC .AND. (.NOT. LSVEC) .AND. ( NR.EQ. N ) ) THEN -* - CALL ZLASET( 'L', N-1,N-1, CZERO, CZERO, A(2,1), LDA ) -* - CALL ZGESVJ( 'U','N','V', N, N, A, LDA, SVA, N, V, LDV, - $ CWORK, LWORK, RWORK, LRWORK, INFO ) - SCALEM = RWORK(1) - NUMRANK = NINT(RWORK(2)) - CALL ZLAPMR( .FALSE., N, N, V, LDV, IWORK ) -* - ELSE IF ( LSVEC .AND. ( .NOT. RSVEC ) ) THEN -* -* .. Singular Values and Left Singular Vectors .. -* -* .. second preconditioning step to avoid need to accumulate -* Jacobi rotations in the Jacobi iterations. - DO 1965 p = 1, NR - CALL ZCOPY( N-p+1, A(p,p), LDA, U(p,p), 1 ) - CALL ZLACGV( N-p+1, U(p,p), 1 ) - 1965 CONTINUE - CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU ) -* - CALL ZGEQRF( N, NR, U, LDU, CWORK(N+1), CWORK(2*N+1), - $ LWORK-2*N, IERR ) -* - DO 1967 p = 1, NR - 1 - CALL ZCOPY( NR-p, U(p,p+1), LDU, U(p+1,p), 1 ) - CALL ZLACGV( N-p+1, U(p,p), 1 ) - 1967 CONTINUE - CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU ) -* - CALL ZGESVJ( 'L', 'U', 'N', NR,NR, U, LDU, SVA, NR, A, - $ LDA, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO ) - SCALEM = RWORK(1) - NUMRANK = NINT(RWORK(2)) -* - IF ( NR .LT. M ) THEN - CALL ZLASET( 'A', M-NR, NR,CZERO, CZERO, U(NR+1,1), LDU ) - IF ( NR .LT. N1 ) THEN - CALL ZLASET( 'A',NR, N1-NR, CZERO, CZERO, U(1,NR+1),LDU ) - CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,U(NR+1,NR+1),LDU ) - END IF - END IF -* - CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U, - $ LDU, CWORK(N+1), LWORK-N, IERR ) -* - IF ( ROWPIV ) - $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 ) -* - DO 1974 p = 1, N1 - XSC = ONE / DZNRM2( M, U(1,p), 1 ) - CALL ZDSCAL( M, XSC, U(1,p), 1 ) - 1974 CONTINUE -* - IF ( TRANSP ) THEN - CALL ZLACPY( 'A', N, N, U, LDU, V, LDV ) - END IF -* - ELSE -* -* .. Full SVD .. -* - IF ( .NOT. JRACC ) THEN -* - IF ( .NOT. ALMORT ) THEN -* -* Second Preconditioning Step (QRF [with pivoting]) -* Note that the composition of TRANSPOSE, QRF and TRANSPOSE is -* equivalent to an LQF CALL. Since in many libraries the QRF -* seems to be better optimized than the LQF, we do explicit -* transpose and use the QRF. This is subject to changes in an -* optimized implementation of ZGEJSV. -* - DO 1968 p = 1, NR - CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) - CALL ZLACGV( N-p+1, V(p,p), 1 ) - 1968 CONTINUE -* -* .. the following two loops perturb small entries to avoid -* denormals in the second QR factorization, where they are -* as good as zeros. This is done to avoid painfully slow -* computation with denormals. The relative size of the perturbation -* is a parameter that can be changed by the implementer. -* This perturbation device will be obsolete on machines with -* properly implemented arithmetic. -* To switch it off, set L2PERT=.FALSE. To remove it from the -* code, remove the action under L2PERT=.TRUE., leave the ELSE part. -* The following two loops should be blocked and fused with the -* transposed copy above. -* - IF ( L2PERT ) THEN - XSC = SQRT(SMALL) - DO 2969 q = 1, NR - CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO) - DO 2968 p = 1, N - IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 ) - $ .OR. ( p .LT. q ) ) -* $ V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) ) - $ V(p,q) = CTEMP - IF ( p .LT. q ) V(p,q) = - V(p,q) - 2968 CONTINUE - 2969 CONTINUE - ELSE - CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV ) - END IF -* -* Estimate the row scaled condition number of R1 -* (If R1 is rectangular, N > NR, then the condition number -* of the leading NR x NR submatrix is estimated.) -* - CALL ZLACPY( 'L', NR, NR, V, LDV, CWORK(2*N+1), NR ) - DO 3950 p = 1, NR - TEMP1 = DZNRM2(NR-p+1,CWORK(2*N+(p-1)*NR+p),1) - CALL ZDSCAL(NR-p+1,ONE/TEMP1,CWORK(2*N+(p-1)*NR+p),1) - 3950 CONTINUE - CALL ZPOCON('L',NR,CWORK(2*N+1),NR,ONE,TEMP1, - $ CWORK(2*N+NR*NR+1),RWORK,IERR) - CONDR1 = ONE / SQRT(TEMP1) -* .. here need a second opinion on the condition number -* .. then assume worst case scenario -* R1 is OK for inverse <=> CONDR1 .LT. DBLE(N) -* more conservative <=> CONDR1 .LT. SQRT(DBLE(N)) -* - COND_OK = SQRT(SQRT(DBLE(NR))) -*[TP] COND_OK is a tuning parameter. -* - IF ( CONDR1 .LT. COND_OK ) THEN -* .. the second QRF without pivoting. Note: in an optimized -* implementation, this QRF should be implemented as the QRF -* of a lower triangular matrix. -* R1^* = Q2 * R2 - CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1), - $ LWORK-2*N, IERR ) -* - IF ( L2PERT ) THEN - XSC = SQRT(SMALL)/EPSLN - DO 3959 p = 2, NR - DO 3958 q = 1, p - 1 - CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))), - $ ZERO) - IF ( ABS(V(q,p)) .LE. TEMP1 ) -* $ V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) ) - $ V(q,p) = CTEMP - 3958 CONTINUE - 3959 CONTINUE - END IF -* - IF ( NR .NE. N ) - $ CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N ) -* .. save ... -* -* .. this transposed copy should be better than naive - DO 1969 p = 1, NR - 1 - CALL ZCOPY( NR-p, V(p,p+1), LDV, V(p+1,p), 1 ) - CALL ZLACGV(NR-p+1, V(p,p), 1 ) - 1969 CONTINUE - V(NR,NR)=CONJG(V(NR,NR)) -* - CONDR2 = CONDR1 -* - ELSE -* -* .. ill-conditioned case: second QRF with pivoting -* Note that windowed pivoting would be equally good -* numerically, and more run-time efficient. So, in -* an optimal implementation, the next call to ZGEQP3 -* should be replaced with eg. CALL ZGEQPX (ACM TOMS #782) -* with properly (carefully) chosen parameters. -* -* R1^* * P2 = Q2 * R2 - DO 3003 p = 1, NR - IWORK(N+p) = 0 - 3003 CONTINUE - CALL ZGEQP3( N, NR, V, LDV, IWORK(N+1), CWORK(N+1), - $ CWORK(2*N+1), LWORK-2*N, RWORK, IERR ) -** CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1), -** $ LWORK-2*N, IERR ) - IF ( L2PERT ) THEN - XSC = SQRT(SMALL) - DO 3969 p = 2, NR - DO 3968 q = 1, p - 1 - CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))), - $ ZERO) - IF ( ABS(V(q,p)) .LE. TEMP1 ) -* $ V(q,p) = TEMP1 * ( V(q,p) / ABS(V(q,p)) ) - $ V(q,p) = CTEMP - 3968 CONTINUE - 3969 CONTINUE - END IF -* - CALL ZLACPY( 'A', N, NR, V, LDV, CWORK(2*N+1), N ) -* - IF ( L2PERT ) THEN - XSC = SQRT(SMALL) - DO 8970 p = 2, NR - DO 8971 q = 1, p - 1 - CTEMP=DCMPLX(XSC*MIN(ABS(V(p,p)),ABS(V(q,q))), - $ ZERO) -* V(p,q) = - TEMP1*( V(q,p) / ABS(V(q,p)) ) - V(p,q) = - CTEMP - 8971 CONTINUE - 8970 CONTINUE - ELSE - CALL ZLASET( 'L',NR-1,NR-1,CZERO,CZERO,V(2,1),LDV ) - END IF -* Now, compute R2 = L3 * Q3, the LQ factorization. - CALL ZGELQF( NR, NR, V, LDV, CWORK(2*N+N*NR+1), - $ CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, IERR ) -* .. and estimate the condition number - CALL ZLACPY( 'L',NR,NR,V,LDV,CWORK(2*N+N*NR+NR+1),NR ) - DO 4950 p = 1, NR - TEMP1 = DZNRM2( p, CWORK(2*N+N*NR+NR+p), NR ) - CALL ZDSCAL( p, ONE/TEMP1, CWORK(2*N+N*NR+NR+p), NR ) - 4950 CONTINUE - CALL ZPOCON( 'L',NR,CWORK(2*N+N*NR+NR+1),NR,ONE,TEMP1, - $ CWORK(2*N+N*NR+NR+NR*NR+1),RWORK,IERR ) - CONDR2 = ONE / SQRT(TEMP1) -* -* - IF ( CONDR2 .GE. COND_OK ) THEN -* .. save the Householder vectors used for Q3 -* (this overwrites the copy of R2, as it will not be -* needed in this branch, but it does not overwritte the -* Huseholder vectors of Q2.). - CALL ZLACPY( 'U', NR, NR, V, LDV, CWORK(2*N+1), N ) -* .. and the rest of the information on Q3 is in -* WORK(2*N+N*NR+1:2*N+N*NR+N) - END IF -* - END IF -* - IF ( L2PERT ) THEN - XSC = SQRT(SMALL) - DO 4968 q = 2, NR - CTEMP = XSC * V(q,q) - DO 4969 p = 1, q - 1 -* V(p,q) = - TEMP1*( V(p,q) / ABS(V(p,q)) ) - V(p,q) = - CTEMP - 4969 CONTINUE - 4968 CONTINUE - ELSE - CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV ) - END IF -* -* Second preconditioning finished; continue with Jacobi SVD -* The input matrix is lower trinagular. -* -* Recover the right singular vectors as solution of a well -* conditioned triangular matrix equation. -* - IF ( CONDR1 .LT. COND_OK ) THEN -* - CALL ZGESVJ( 'L','U','N',NR,NR,V,LDV,SVA,NR,U, LDU, - $ CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,RWORK, - $ LRWORK, INFO ) - SCALEM = RWORK(1) - NUMRANK = NINT(RWORK(2)) - DO 3970 p = 1, NR - CALL ZCOPY( NR, V(1,p), 1, U(1,p), 1 ) - CALL ZDSCAL( NR, SVA(p), V(1,p), 1 ) - 3970 CONTINUE - -* .. pick the right matrix equation and solve it -* - IF ( NR .EQ. N ) THEN -* :)) .. best case, R1 is inverted. The solution of this matrix -* equation is Q2*V2 = the product of the Jacobi rotations -* used in ZGESVJ, premultiplied with the orthogonal matrix -* from the second QR factorization. - CALL ZTRSM('L','U','N','N', NR,NR,CONE, A,LDA, V,LDV) - ELSE -* .. R1 is well conditioned, but non-square. Adjoint of R2 -* is inverted to get the product of the Jacobi rotations -* used in ZGESVJ. The Q-factor from the second QR -* factorization is then built in explicitly. - CALL ZTRSM('L','U','C','N',NR,NR,CONE,CWORK(2*N+1), - $ N,V,LDV) - IF ( NR .LT. N ) THEN - CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV) - CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV) - CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV) - END IF - CALL ZUNMQR('L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1), - $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR) - END IF -* - ELSE IF ( CONDR2 .LT. COND_OK ) THEN -* -* The matrix R2 is inverted. The solution of the matrix equation -* is Q3^* * V3 = the product of the Jacobi rotations (appplied to -* the lower triangular L3 from the LQ factorization of -* R2=L3*Q3), pre-multiplied with the transposed Q3. - CALL ZGESVJ( 'L', 'U', 'N', NR, NR, V, LDV, SVA, NR, U, - $ LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, - $ RWORK, LRWORK, INFO ) - SCALEM = RWORK(1) - NUMRANK = NINT(RWORK(2)) - DO 3870 p = 1, NR - CALL ZCOPY( NR, V(1,p), 1, U(1,p), 1 ) - CALL ZDSCAL( NR, SVA(p), U(1,p), 1 ) - 3870 CONTINUE - CALL ZTRSM('L','U','N','N',NR,NR,CONE,CWORK(2*N+1),N, - $ U,LDU) -* .. apply the permutation from the second QR factorization - DO 873 q = 1, NR - DO 872 p = 1, NR - CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q) - 872 CONTINUE - DO 874 p = 1, NR - U(p,q) = CWORK(2*N+N*NR+NR+p) - 874 CONTINUE - 873 CONTINUE - IF ( NR .LT. N ) THEN - CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV ) - CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV ) - CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV) - END IF - CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1), - $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) - ELSE -* Last line of defense. -* #:( This is a rather pathological case: no scaled condition -* improvement after two pivoted QR factorizations. Other -* possibility is that the rank revealing QR factorization -* or the condition estimator has failed, or the COND_OK -* is set very close to ONE (which is unnecessary). Normally, -* this branch should never be executed, but in rare cases of -* failure of the RRQR or condition estimator, the last line of -* defense ensures that ZGEJSV completes the task. -* Compute the full SVD of L3 using ZGESVJ with explicit -* accumulation of Jacobi rotations. - CALL ZGESVJ( 'L', 'U', 'V', NR, NR, V, LDV, SVA, NR, U, - $ LDU, CWORK(2*N+N*NR+NR+1), LWORK-2*N-N*NR-NR, - $ RWORK, LRWORK, INFO ) - SCALEM = RWORK(1) - NUMRANK = NINT(RWORK(2)) - IF ( NR .LT. N ) THEN - CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV ) - CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV ) - CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV) - END IF - CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1), - $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) -* - CALL ZUNMLQ( 'L', 'C', NR, NR, NR, CWORK(2*N+1), N, - $ CWORK(2*N+N*NR+1), U, LDU, CWORK(2*N+N*NR+NR+1), - $ LWORK-2*N-N*NR-NR, IERR ) - DO 773 q = 1, NR - DO 772 p = 1, NR - CWORK(2*N+N*NR+NR+IWORK(N+p)) = U(p,q) - 772 CONTINUE - DO 774 p = 1, NR - U(p,q) = CWORK(2*N+N*NR+NR+p) - 774 CONTINUE - 773 CONTINUE -* - END IF -* -* Permute the rows of V using the (column) permutation from the -* first QRF. Also, scale the columns to make them unit in -* Euclidean norm. This applies to all cases. -* - TEMP1 = SQRT(DBLE(N)) * EPSLN - DO 1972 q = 1, N - DO 972 p = 1, N - CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q) - 972 CONTINUE - DO 973 p = 1, N - V(p,q) = CWORK(2*N+N*NR+NR+p) - 973 CONTINUE - XSC = ONE / DZNRM2( N, V(1,q), 1 ) - IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) - $ CALL ZDSCAL( N, XSC, V(1,q), 1 ) - 1972 CONTINUE -* At this moment, V contains the right singular vectors of A. -* Next, assemble the left singular vector matrix U (M x N). - IF ( NR .LT. M ) THEN - CALL ZLASET('A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU) - IF ( NR .LT. N1 ) THEN - CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU) - CALL ZLASET('A',M-NR,N1-NR,CZERO,CONE, - $ U(NR+1,NR+1),LDU) - END IF - END IF -* -* The Q matrix from the first QRF is built into the left singular -* matrix U. This applies to all cases. -* - CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U, - $ LDU, CWORK(N+1), LWORK-N, IERR ) - -* The columns of U are normalized. The cost is O(M*N) flops. - TEMP1 = SQRT(DBLE(M)) * EPSLN - DO 1973 p = 1, NR - XSC = ONE / DZNRM2( M, U(1,p), 1 ) - IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) - $ CALL ZDSCAL( M, XSC, U(1,p), 1 ) - 1973 CONTINUE -* -* If the initial QRF is computed with row pivoting, the left -* singular vectors must be adjusted. -* - IF ( ROWPIV ) - $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 ) -* - ELSE -* -* .. the initial matrix A has almost orthogonal columns and -* the second QRF is not needed -* - CALL ZLACPY( 'U', N, N, A, LDA, CWORK(N+1), N ) - IF ( L2PERT ) THEN - XSC = SQRT(SMALL) - DO 5970 p = 2, N - CTEMP = XSC * CWORK( N + (p-1)*N + p ) - DO 5971 q = 1, p - 1 -* CWORK(N+(q-1)*N+p)=-TEMP1 * ( CWORK(N+(p-1)*N+q) / -* $ ABS(CWORK(N+(p-1)*N+q)) ) - CWORK(N+(q-1)*N+p)=-CTEMP - 5971 CONTINUE - 5970 CONTINUE - ELSE - CALL ZLASET( 'L',N-1,N-1,CZERO,CZERO,CWORK(N+2),N ) - END IF -* - CALL ZGESVJ( 'U', 'U', 'N', N, N, CWORK(N+1), N, SVA, - $ N, U, LDU, CWORK(N+N*N+1), LWORK-N-N*N, RWORK, LRWORK, - $ INFO ) -* - SCALEM = RWORK(1) - NUMRANK = NINT(RWORK(2)) - DO 6970 p = 1, N - CALL ZCOPY( N, CWORK(N+(p-1)*N+1), 1, U(1,p), 1 ) - CALL ZDSCAL( N, SVA(p), CWORK(N+(p-1)*N+1), 1 ) - 6970 CONTINUE -* - CALL ZTRSM( 'L', 'U', 'N', 'N', N, N, - $ CONE, A, LDA, CWORK(N+1), N ) - DO 6972 p = 1, N - CALL ZCOPY( N, CWORK(N+p), N, V(IWORK(p),1), LDV ) - 6972 CONTINUE - TEMP1 = SQRT(DBLE(N))*EPSLN - DO 6971 p = 1, N - XSC = ONE / DZNRM2( N, V(1,p), 1 ) - IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) - $ CALL ZDSCAL( N, XSC, V(1,p), 1 ) - 6971 CONTINUE -* -* Assemble the left singular vector matrix U (M x N). -* - IF ( N .LT. M ) THEN - CALL ZLASET( 'A', M-N, N, CZERO, CZERO, U(N+1,1), LDU ) - IF ( N .LT. N1 ) THEN - CALL ZLASET('A',N, N1-N, CZERO, CZERO, U(1,N+1),LDU) - CALL ZLASET( 'A',M-N,N1-N, CZERO, CONE,U(N+1,N+1),LDU) - END IF - END IF - CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U, - $ LDU, CWORK(N+1), LWORK-N, IERR ) - TEMP1 = SQRT(DBLE(M))*EPSLN - DO 6973 p = 1, N1 - XSC = ONE / DZNRM2( M, U(1,p), 1 ) - IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) - $ CALL ZDSCAL( M, XSC, U(1,p), 1 ) - 6973 CONTINUE -* - IF ( ROWPIV ) - $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 ) -* - END IF -* -* end of the >> almost orthogonal case << in the full SVD -* - ELSE -* -* This branch deploys a preconditioned Jacobi SVD with explicitly -* accumulated rotations. It is included as optional, mainly for -* experimental purposes. It does perform well, and can also be used. -* In this implementation, this branch will be automatically activated -* if the condition number sigma_max(A) / sigma_min(A) is predicted -* to be greater than the overflow threshold. This is because the -* a posteriori computation of the singular vectors assumes robust -* implementation of BLAS and some LAPACK procedures, capable of working -* in presence of extreme values, e.g. when the singular values spread from -* the underflow to the overflow threshold. -* - DO 7968 p = 1, NR - CALL ZCOPY( N-p+1, A(p,p), LDA, V(p,p), 1 ) - CALL ZLACGV( N-p+1, V(p,p), 1 ) - 7968 CONTINUE -* - IF ( L2PERT ) THEN - XSC = SQRT(SMALL/EPSLN) - DO 5969 q = 1, NR - CTEMP = DCMPLX(XSC*ABS( V(q,q) ),ZERO) - DO 5968 p = 1, N - IF ( ( p .GT. q ) .AND. ( ABS(V(p,q)) .LE. TEMP1 ) - $ .OR. ( p .LT. q ) ) -* $ V(p,q) = TEMP1 * ( V(p,q) / ABS(V(p,q)) ) - $ V(p,q) = CTEMP - IF ( p .LT. q ) V(p,q) = - V(p,q) - 5968 CONTINUE - 5969 CONTINUE - ELSE - CALL ZLASET( 'U', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV ) - END IF - - CALL ZGEQRF( N, NR, V, LDV, CWORK(N+1), CWORK(2*N+1), - $ LWORK-2*N, IERR ) - CALL ZLACPY( 'L', N, NR, V, LDV, CWORK(2*N+1), N ) -* - DO 7969 p = 1, NR - CALL ZCOPY( NR-p+1, V(p,p), LDV, U(p,p), 1 ) - CALL ZLACGV( NR-p+1, U(p,p), 1 ) - 7969 CONTINUE - - IF ( L2PERT ) THEN - XSC = SQRT(SMALL/EPSLN) - DO 9970 q = 2, NR - DO 9971 p = 1, q - 1 - CTEMP = DCMPLX(XSC * MIN(ABS(U(p,p)),ABS(U(q,q))), - $ ZERO) -* U(p,q) = - TEMP1 * ( U(q,p) / ABS(U(q,p)) ) - U(p,q) = - CTEMP - 9971 CONTINUE - 9970 CONTINUE - ELSE - CALL ZLASET('U', NR-1, NR-1, CZERO, CZERO, U(1,2), LDU ) - END IF - - CALL ZGESVJ( 'L', 'U', 'V', NR, NR, U, LDU, SVA, - $ N, V, LDV, CWORK(2*N+N*NR+1), LWORK-2*N-N*NR, - $ RWORK, LRWORK, INFO ) - SCALEM = RWORK(1) - NUMRANK = NINT(RWORK(2)) - - IF ( NR .LT. N ) THEN - CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV ) - CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV ) - CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV ) - END IF - - CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1), - $ V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR ) -* -* Permute the rows of V using the (column) permutation from the -* first QRF. Also, scale the columns to make them unit in -* Euclidean norm. This applies to all cases. -* - TEMP1 = SQRT(DBLE(N)) * EPSLN - DO 7972 q = 1, N - DO 8972 p = 1, N - CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q) - 8972 CONTINUE - DO 8973 p = 1, N - V(p,q) = CWORK(2*N+N*NR+NR+p) - 8973 CONTINUE - XSC = ONE / DZNRM2( N, V(1,q), 1 ) - IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) ) - $ CALL ZDSCAL( N, XSC, V(1,q), 1 ) - 7972 CONTINUE -* -* At this moment, V contains the right singular vectors of A. -* Next, assemble the left singular vector matrix U (M x N). -* - IF ( NR .LT. M ) THEN - CALL ZLASET( 'A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU ) - IF ( NR .LT. N1 ) THEN - CALL ZLASET('A',NR, N1-NR, CZERO, CZERO, U(1,NR+1),LDU) - CALL ZLASET('A',M-NR,N1-NR, CZERO, CONE,U(NR+1,NR+1),LDU) - END IF - END IF -* - CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U, - $ LDU, CWORK(N+1), LWORK-N, IERR ) -* - IF ( ROWPIV ) - $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(IWOFF+1), -1 ) -* -* - END IF - IF ( TRANSP ) THEN -* .. swap U and V because the procedure worked on A^* - DO 6974 p = 1, N - CALL ZSWAP( N, U(1,p), 1, V(1,p), 1 ) - 6974 CONTINUE - END IF -* - END IF -* end of the full SVD -* -* Undo scaling, if necessary (and possible) -* - IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN - CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR ) - USCAL1 = ONE - USCAL2 = ONE - END IF -* - IF ( NR .LT. N ) THEN - DO 3004 p = NR+1, N - SVA(p) = ZERO - 3004 CONTINUE - END IF -* - RWORK(1) = USCAL2 * SCALEM - RWORK(2) = USCAL1 - IF ( ERREST ) RWORK(3) = SCONDA - IF ( LSVEC .AND. RSVEC ) THEN - RWORK(4) = CONDR1 - RWORK(5) = CONDR2 - END IF - IF ( L2TRAN ) THEN - RWORK(6) = ENTRA - RWORK(7) = ENTRAT - END IF -* - IWORK(1) = NR - IWORK(2) = NUMRANK - IWORK(3) = WARNING - IF ( TRANSP ) THEN - IWORK(4) = 1 - ELSE - IWORK(4) = -1 - END IF - -* - RETURN -* .. -* .. END OF ZGEJSV -* .. - END -*