*> \brief ZGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download ZGESVDQ + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE ZGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA, * S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK, * CWORK, LCWORK, RWORK, LRWORK, INFO ) * * .. Scalar Arguments .. * IMPLICIT NONE * CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV * INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LCWORK, LRWORK, * INFO * .. * .. Array Arguments .. * COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * ) * DOUBLE PRECISION S( * ), RWORK( * ) * INTEGER IWORK( * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> ZCGESVDQ computes the singular value decomposition (SVD) of a complex *> M-by-N matrix A, where M >= N. The SVD of A is written as *> [++] [xx] [x0] [xx] *> A = U * SIGMA * V^*, [++] = [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 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. *> \endverbatim * * Arguments * ========= * *> \param[in] JOBA *> \verbatim *> JOBA is CHARACTER*1 *> Specifies the level of accuracy in the computed SVD *> = 'A' The requested accuracy corresponds to having the backward *> error bounded by || delta A ||_F <= f(m,n) * EPS * || A ||_F, *> where EPS = DLAMCH('Epsilon'). This authorises ZGESVDQ to *> truncate the computed triangular factor in a rank revealing *> QR factorization whenever the truncated part is below the *> threshold of the order of EPS * ||A||_F. This is aggressive *> truncation level. *> = 'M' Similarly as with 'A', but the truncation is more gentle: it *> is allowed only when there is a drop on the diagonal of the *> triangular factor in the QR factorization. This is medium *> truncation level. *> = 'H' High accuracy requested. No numerical rank determination based *> on the rank revealing QR factorization is attempted. *> = 'E' Same as 'H', and in addition the condition number of column *> scaled A is estimated and returned in RWORK(1). *> N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1) *> \endverbatim *> *> \param[in] JOBP *> \verbatim *> JOBP is CHARACTER*1 *> = 'P' The rows of A are ordered in decreasing order with respect to *> ||A(i,:)||_\infty. This enhances numerical accuracy at the cost *> of extra data movement. Recommended for numerical robustness. *> = 'N' No row pivoting. *> \endverbatim *> *> \param[in] JOBR *> \verbatim *> JOBR is CHARACTER*1 *> = 'T' After the initial pivoted QR factorization, ZGESVD is applied to *> the adjoint R**H of the computed triangular factor R. This involves *> some extra data movement (matrix transpositions). Useful for *> experiments, research and development. *> = 'N' The triangular factor R is given as input to CGESVD. This may be *> preferred as it involves less data movement. *> \endverbatim *> *> \param[in] JOBU *> \verbatim *> JOBU is CHARACTER*1 *> = 'A' All M left singular vectors are computed and returned in the *> matrix U. See the description of U. *> = 'S' or 'U' N = min(M,N) left singular vectors are computed and returned *> in the matrix U. See the description of U. *> = 'R' Numerical rank NUMRANK is determined and only NUMRANK left singular *> vectors are computed and returned in the matrix U. *> = 'F' The N left singular vectors are returned in factored form as the *> product of the Q factor from the initial QR factorization and the *> N left singular vectors of (R**H , 0)**H. If row pivoting is used, *> then the necessary information on the row pivoting is stored in *> IWORK(N+1:N+M-1). *> = 'N' The left singular vectors are not computed. *> \endverbatim *> *> \param[in] JOBV *> \verbatim *> JOBV is CHARACTER*1 *> = 'A', 'V' All N right singular vectors are computed and returned in *> the matrix V. *> = 'R' Numerical rank NUMRANK is determined and only NUMRANK right singular *> vectors are computed and returned in the matrix V. This option is *> allowed only if JOBU = 'R' or JOBU = 'N'; otherwise it is illegal. *> = 'N' The right singular vectors are not computed. *> \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 of dimensions LDA x N *> On entry, the input matrix A. *> On exit, if JOBU .NE. 'N' or JOBV .NE. 'N', the lower triangle of A contains *> the Householder vectors as stored by ZGEQP3. If JOBU = 'F', these Householder *> vectors together with CWORK(1:N) can be used to restore the Q factors from *> the initial pivoted QR factorization of A. See the description of U. *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is INTEGER. *> The leading dimension of the array A. LDA >= max(1,M). *> \endverbatim *> *> \param[out] S *> \verbatim *> S is DOUBLE PRECISION array of dimension N. *> The singular values of A, ordered so that S(i) >= S(i+1). *> \endverbatim *> *> \param[out] U *> \verbatim *> U is COMPLEX*16 array, dimension *> LDU x M if JOBU = 'A'; see the description of LDU. In this case, *> on exit, U contains the M left singular vectors. *> LDU x N if JOBU = 'S', 'U', 'R' ; see the description of LDU. In this *> case, U contains the leading N or the leading NUMRANK left singular vectors. *> LDU x N if JOBU = 'F' ; see the description of LDU. In this case U *> contains N x N unitary matrix that can be used to form the left *> singular vectors. *> If JOBU = 'N', U is not referenced. *> \endverbatim *> *> \param[in] LDU *> \verbatim *> LDU is INTEGER. *> The leading dimension of the array U. *> If JOBU = 'A', 'S', 'U', 'R', LDU >= max(1,M). *> If JOBU = 'F', LDU >= max(1,N). *> Otherwise, LDU >= 1. *> \endverbatim *> *> \param[out] V *> \verbatim *> V is COMPLEX*16 array, dimension *> LDV x N if JOBV = 'A', 'V', 'R' or if JOBA = 'E' . *> If JOBV = 'A', or 'V', V contains the N-by-N unitary matrix V**H; *> If JOBV = 'R', V contains the first NUMRANK rows of V**H (the right *> singular vectors, stored rowwise, of the NUMRANK largest singular values). *> If JOBV = 'N' and JOBA = 'E', V is used as a workspace. *> If JOBV = 'N', and JOBA.NE.'E', V is not referenced. *> \endverbatim *> *> \param[in] LDV *> \verbatim *> LDV is INTEGER *> The leading dimension of the array V. *> If JOBV = 'A', 'V', 'R', or JOBA = 'E', LDV >= max(1,N). *> Otherwise, LDV >= 1. *> \endverbatim *> *> \param[out] NUMRANK *> \verbatim *> NUMRANK is INTEGER *> NUMRANK is the numerical rank first determined after the rank *> revealing QR factorization, following the strategy specified by the *> value of JOBA. If JOBV = 'R' and JOBU = 'R', only NUMRANK *> leading singular values and vectors are then requested in the call *> of CGESVD. The final value of NUMRANK might be further reduced if *> some singular values are computed as zeros. *> \endverbatim *> *> \param[out] IWORK *> \verbatim *> IWORK is INTEGER array, dimension (max(1, LIWORK)). *> On exit, IWORK(1:N) contains column pivoting permutation of the *> rank revealing QR factorization. *> If JOBP = 'P', IWORK(N+1:N+M-1) contains the indices of the sequence *> of row swaps used in row pivoting. These can be used to restore the *> left singular vectors in the case JOBU = 'F'. *> *> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0, *> IWORK(1) returns the minimal LIWORK. *> \endverbatim *> *> \param[in] LIWORK *> \verbatim *> LIWORK is INTEGER *> The dimension of the array IWORK. *> LIWORK >= N + M - 1, if JOBP = 'P'; *> LIWORK >= N if JOBP = 'N'. *> *> If LIWORK = -1, then a workspace query is assumed; the routine *> only calculates and returns the optimal and minimal sizes *> for the CWORK, IWORK, and RWORK arrays, and no error *> message related to LCWORK is issued by XERBLA. *> \endverbatim *> *> \param[out] CWORK *> \verbatim *> CWORK is COMPLEX*12 array, dimension (max(2, LCWORK)), used as a workspace. *> On exit, if, on entry, LCWORK.NE.-1, CWORK(1:N) contains parameters *> needed to recover the Q factor from the QR factorization computed by *> ZGEQP3. *> *> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0, *> CWORK(1) returns the optimal LCWORK, and *> CWORK(2) returns the minimal LCWORK. *> \endverbatim *> *> \param[in,out] LCWORK *> \verbatim *> LCWORK is INTEGER *> The dimension of the array CWORK. It is determined as follows: *> Let LWQP3 = N+1, LWCON = 2*N, and let *> LWUNQ = { MAX( N, 1 ), if JOBU = 'R', 'S', or 'U' *> { MAX( M, 1 ), if JOBU = 'A' *> LWSVD = MAX( 3*N, 1 ) *> LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 3*(N/2), 1 ), LWUNLQ = MAX( N, 1 ), *> LWQRF = MAX( N/2, 1 ), LWUNQ2 = MAX( N, 1 ) *> Then the minimal value of LCWORK is: *> = MAX( N + LWQP3, LWSVD ) if only the singular values are needed; *> = MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values are needed, *> and a scaled condition estimate requested; *> *> = N + MAX( LWQP3, LWSVD, LWUNQ ) if the singular values and the left *> singular vectors are requested; *> = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the singular values and the left *> singular vectors are requested, and also *> a scaled condition estimate requested; *> *> = N + MAX( LWQP3, LWSVD ) if the singular values and the right *> singular vectors are requested; *> = N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and the right *> singular vectors are requested, and also *> a scaled condition etimate requested; *> *> = N + MAX( LWQP3, LWSVD, LWUNQ ) if the full SVD is requested with JOBV = 'R'; *> independent of JOBR; *> = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the full SVD is requested, *> JOBV = 'R' and, also a scaled condition *> estimate requested; independent of JOBR; *> = MAX( N + MAX( LWQP3, LWSVD, LWUNQ ), *> N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ) ) if the *> full SVD is requested with JOBV = 'A' or 'V', and *> JOBR ='N' *> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ), *> N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ ) ) *> if the full SVD is requested with JOBV = 'A' or 'V', and *> JOBR ='N', and also a scaled condition number estimate *> requested. *> = MAX( N + MAX( LWQP3, LWSVD, LWUNQ ), *> N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) ) if the *> full SVD is requested with JOBV = 'A', 'V', and JOBR ='T' *> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ), *> N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) ) *> if the full SVD is requested with JOBV = 'A', 'V' and *> JOBR ='T', and also a scaled condition number estimate *> requested. *> Finally, LCWORK must be at least two: LCWORK = MAX( 2, LCWORK ). *> *> If LCWORK = -1, then a workspace query is assumed; the routine *> only calculates and returns the optimal and minimal sizes *> for the CWORK, IWORK, and RWORK arrays, and no error *> message related to LCWORK is issued by XERBLA. *> \endverbatim *> *> \param[out] RWORK *> \verbatim *> RWORK is DOUBLE PRECISION array, dimension (max(1, LRWORK)). *> On exit, *> 1. If JOBA = 'E', RWORK(1) contains an estimate of the condition *> number of column scaled A. If A = C * D where D is diagonal and C *> has unit columns in the Euclidean norm, then, assuming full column rank, *> N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1). *> Otherwise, RWORK(1) = -1. *> 2. RWORK(2) contains the number of singular values computed as *> exact zeros in ZGESVD applied to the upper triangular or trapezoidal *> R (from the initial QR factorization). In case of early exit (no call to *> ZGESVD, such as in the case of zero matrix) RWORK(2) = -1. *> *> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0, *> RWORK(1) returns the minimal LRWORK. *> \endverbatim *> *> \param[in] LRWORK *> \verbatim *> LRWORK is INTEGER. *> The dimension of the array RWORK. *> If JOBP ='P', then LRWORK >= MAX(2, M, 5*N); *> Otherwise, LRWORK >= MAX(2, 5*N). *> *> If LRWORK = -1, then a workspace query is assumed; the routine *> only calculates and returns the optimal and minimal sizes *> for the CWORK, IWORK, and RWORK arrays, and no error *> message related to LCWORK is issued by XERBLA. *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is INTEGER *> = 0: successful exit. *> < 0: if INFO = -i, the i-th argument had an illegal value. *> > 0: if ZBDSQR did not converge, INFO specifies how many superdiagonals *> of an intermediate bidiagonal form B (computed in ZGESVD) did not *> converge to zero. *> \endverbatim * *> \par Further Details: * ======================== *> *> \verbatim *> *> 1. The data movement (matrix transpose) is coded using simple nested *> DO-loops because BLAS and LAPACK do not provide corresponding subroutines. *> Those DO-loops are easily identified in this source code - by the CONTINUE *> statements labeled with 11**. In an optimized version of this code, the *> nested DO loops should be replaced with calls to an optimized subroutine. *> 2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause *> column norm overflow. This is the minial precaution and it is left to the *> SVD routine (CGESVD) to do its own preemptive scaling if potential over- *> or underflows are detected. To avoid repeated scanning of the array A, *> an optimal implementation would do all necessary scaling before calling *> CGESVD and the scaling in CGESVD can be switched off. *> 3. Other comments related to code optimization are given in comments in the *> code, enlosed in [[double brackets]]. *> \endverbatim * *> \par Bugs, examples and comments * =========================== * *> \verbatim *> Please report all bugs and send interesting examples and/or comments to *> drmac@math.hr. Thank you. *> \endverbatim * *> \par References * =============== * *> \verbatim *> [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for *> Computing the SVD with High Accuracy. ACM Trans. Math. Softw. *> 44(1): 11:1-11:30 (2017) *> *> SIGMA library, xGESVDQ section updated February 2016. *> Developed and coded by Zlatko Drmac, Department of Mathematics *> University of Zagreb, Croatia, drmac@math.hr *> \endverbatim * * *> \par Contributors: * ================== *> *> \verbatim *> Developed and coded by Zlatko Drmac, Department of Mathematics *> University of Zagreb, Croatia, drmac@math.hr *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \ingroup complex16GEsing * * ===================================================================== SUBROUTINE ZGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA, $ S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK, $ CWORK, LCWORK, RWORK, LRWORK, INFO ) * .. Scalar Arguments .. IMPLICIT NONE CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LCWORK, LRWORK, $ INFO * .. * .. Array Arguments .. COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * ) DOUBLE PRECISION S( * ), RWORK( * ) INTEGER IWORK( * ) * * ===================================================================== * * .. 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 .. INTEGER IERR, NR, N1, OPTRATIO, p, q INTEGER LWCON, LWQP3, LWRK_ZGELQF, LWRK_ZGESVD, LWRK_ZGESVD2, $ LWRK_ZGEQP3, LWRK_ZGEQRF, LWRK_ZUNMLQ, LWRK_ZUNMQR, $ LWRK_ZUNMQR2, LWLQF, LWQRF, LWSVD, LWSVD2, LWUNQ, $ LWUNQ2, LWUNLQ, MINWRK, MINWRK2, OPTWRK, OPTWRK2, $ IMINWRK, RMINWRK LOGICAL ACCLA, ACCLM, ACCLH, ASCALED, CONDA, DNTWU, DNTWV, $ LQUERY, LSVC0, LSVEC, ROWPRM, RSVEC, RTRANS, WNTUA, $ WNTUF, WNTUR, WNTUS, WNTVA, WNTVR DOUBLE PRECISION BIG, EPSLN, RTMP, SCONDA, SFMIN COMPLEX*16 CTMP * .. * .. Local Arrays COMPLEX*16 CDUMMY(1) DOUBLE PRECISION RDUMMY(1) * .. * .. External Subroutines (BLAS, LAPACK) EXTERNAL ZGELQF, ZGEQP3, ZGEQRF, ZGESVD, ZLACPY, ZLAPMT, $ ZLASCL, ZLASET, ZLASWP, ZDSCAL, DLASET, DLASCL, $ ZPOCON, ZUNMLQ, ZUNMQR, XERBLA * .. * .. External Functions (BLAS, LAPACK) LOGICAL LSAME INTEGER IDAMAX DOUBLE PRECISION ZLANGE, DZNRM2, DLAMCH EXTERNAL LSAME, ZLANGE, IDAMAX, DZNRM2, DLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, CONJG, MAX, MIN, DBLE, SQRT * .. * .. Executable Statements .. * * Test the input arguments * WNTUS = LSAME( JOBU, 'S' ) .OR. LSAME( JOBU, 'U' ) WNTUR = LSAME( JOBU, 'R' ) WNTUA = LSAME( JOBU, 'A' ) WNTUF = LSAME( JOBU, 'F' ) LSVC0 = WNTUS .OR. WNTUR .OR. WNTUA LSVEC = LSVC0 .OR. WNTUF DNTWU = LSAME( JOBU, 'N' ) * WNTVR = LSAME( JOBV, 'R' ) WNTVA = LSAME( JOBV, 'A' ) .OR. LSAME( JOBV, 'V' ) RSVEC = WNTVR .OR. WNTVA DNTWV = LSAME( JOBV, 'N' ) * ACCLA = LSAME( JOBA, 'A' ) ACCLM = LSAME( JOBA, 'M' ) CONDA = LSAME( JOBA, 'E' ) ACCLH = LSAME( JOBA, 'H' ) .OR. CONDA * ROWPRM = LSAME( JOBP, 'P' ) RTRANS = LSAME( JOBR, 'T' ) * IF ( ROWPRM ) THEN IMINWRK = MAX( 1, N + M - 1 ) RMINWRK = MAX( 2, M, 5*N ) ELSE IMINWRK = MAX( 1, N ) RMINWRK = MAX( 2, 5*N ) END IF LQUERY = (LIWORK .EQ. -1 .OR. LCWORK .EQ. -1 .OR. LRWORK .EQ. -1) INFO = 0 IF ( .NOT. ( ACCLA .OR. ACCLM .OR. ACCLH ) ) THEN INFO = -1 ELSE IF ( .NOT.( ROWPRM .OR. LSAME( JOBP, 'N' ) ) ) THEN INFO = -2 ELSE IF ( .NOT.( RTRANS .OR. LSAME( JOBR, 'N' ) ) ) THEN INFO = -3 ELSE IF ( .NOT.( LSVEC .OR. DNTWU ) ) THEN INFO = -4 ELSE IF ( WNTUR .AND. WNTVA ) THEN INFO = -5 ELSE IF ( .NOT.( RSVEC .OR. DNTWV )) THEN INFO = -5 ELSE IF ( M.LT.0 ) THEN INFO = -6 ELSE IF ( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN INFO = -7 ELSE IF ( LDA.LT.MAX( 1, M ) ) THEN INFO = -9 ELSE IF ( LDU.LT.1 .OR. ( LSVC0 .AND. LDU.LT.M ) .OR. $ ( WNTUF .AND. LDU.LT.N ) ) THEN INFO = -12 ELSE IF ( LDV.LT.1 .OR. ( RSVEC .AND. LDV.LT.N ) .OR. $ ( CONDA .AND. LDV.LT.N ) ) THEN INFO = -14 ELSE IF ( LIWORK .LT. IMINWRK .AND. .NOT. LQUERY ) THEN INFO = -17 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 are written with a lot of redundancy and * can be simplified. However, this detailed form is easier for * maintenance and modifications of the code.]] * * .. minimal workspace length for ZGEQP3 of an M x N matrix LWQP3 = N+1 * .. minimal workspace length for ZUNMQR to build left singular vectors IF ( WNTUS .OR. WNTUR ) THEN LWUNQ = MAX( N , 1 ) ELSE IF ( WNTUA ) THEN LWUNQ = MAX( M , 1 ) END IF * .. minimal workspace length for ZPOCON of an N x N matrix LWCON = 2 * N * .. ZGESVD of an N x N matrix LWSVD = MAX( 3 * N, 1 ) IF ( LQUERY ) THEN CALL ZGEQP3( M, N, A, LDA, IWORK, CDUMMY, CDUMMY, -1, $ RDUMMY, IERR ) LWRK_ZGEQP3 = INT( CDUMMY(1) ) IF ( WNTUS .OR. WNTUR ) THEN CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U, $ LDU, CDUMMY, -1, IERR ) LWRK_ZUNMQR = INT( CDUMMY(1) ) ELSE IF ( WNTUA ) THEN CALL ZUNMQR( 'L', 'N', M, M, N, A, LDA, CDUMMY, U, $ LDU, CDUMMY, -1, IERR ) LWRK_ZUNMQR = INT( CDUMMY(1) ) ELSE LWRK_ZUNMQR = 0 END IF END IF MINWRK = 2 OPTWRK = 2 IF ( .NOT. (LSVEC .OR. RSVEC ) ) THEN * .. minimal and optimal sizes of the complex workspace if * only the singular values are requested IF ( CONDA ) THEN MINWRK = MAX( N+LWQP3, LWCON, LWSVD ) ELSE MINWRK = MAX( N+LWQP3, LWSVD ) END IF IF ( LQUERY ) THEN CALL ZGESVD( 'N', 'N', N, N, A, LDA, S, U, LDU, $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) LWRK_ZGESVD = INT( CDUMMY(1) ) IF ( CONDA ) THEN OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON, LWRK_ZGESVD ) ELSE OPTWRK = MAX( N+LWRK_ZGEQP3, LWRK_ZGESVD ) END IF END IF 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 ( CONDA ) THEN MINWRK = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) ELSE MINWRK = N + MAX( LWQP3, LWSVD, LWUNQ ) END IF IF ( LQUERY ) THEN IF ( RTRANS ) THEN CALL ZGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU, $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) ELSE CALL ZGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU, $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) END IF LWRK_ZGESVD = INT( CDUMMY(1) ) IF ( CONDA ) THEN OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, LWRK_ZGESVD, $ LWRK_ZUNMQR ) ELSE OPTWRK = N + MAX( LWRK_ZGEQP3, LWRK_ZGESVD, $ LWRK_ZUNMQR ) END IF END IF 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 ( CONDA ) THEN MINWRK = N + MAX( LWQP3, LWCON, LWSVD ) ELSE MINWRK = N + MAX( LWQP3, LWSVD ) END IF IF ( LQUERY ) THEN IF ( RTRANS ) THEN CALL ZGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU, $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) ELSE CALL ZGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU, $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) END IF LWRK_ZGESVD = INT( CDUMMY(1) ) IF ( CONDA ) THEN OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, LWRK_ZGESVD ) ELSE OPTWRK = N + MAX( LWRK_ZGEQP3, LWRK_ZGESVD ) END IF END IF ELSE * .. minimal and optimal sizes of the complex workspace if the * full SVD is requested IF ( RTRANS ) THEN MINWRK = MAX( LWQP3, LWSVD, LWUNQ ) IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON ) MINWRK = MINWRK + N IF ( WNTVA ) THEN * .. minimal workspace length for N x N/2 ZGEQRF LWQRF = MAX( N/2, 1 ) * .. minimal workspace length for N/2 x N/2 ZGESVD LWSVD2 = MAX( 3 * (N/2), 1 ) LWUNQ2 = MAX( N, 1 ) MINWRK2 = MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, $ N/2+LWUNQ2, LWUNQ ) IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON ) MINWRK2 = N + MINWRK2 MINWRK = MAX( MINWRK, MINWRK2 ) END IF ELSE MINWRK = MAX( LWQP3, LWSVD, LWUNQ ) IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON ) MINWRK = MINWRK + N IF ( WNTVA ) THEN * .. minimal workspace length for N/2 x N ZGELQF LWLQF = MAX( N/2, 1 ) LWSVD2 = MAX( 3 * (N/2), 1 ) LWUNLQ = MAX( N , 1 ) MINWRK2 = MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, $ N/2+LWUNLQ, LWUNQ ) IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON ) MINWRK2 = N + MINWRK2 MINWRK = MAX( MINWRK, MINWRK2 ) END IF END IF IF ( LQUERY ) THEN IF ( RTRANS ) THEN CALL ZGESVD( 'O', 'A', N, N, A, LDA, S, U, LDU, $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) LWRK_ZGESVD = INT( CDUMMY(1) ) OPTWRK = MAX(LWRK_ZGEQP3,LWRK_ZGESVD,LWRK_ZUNMQR) IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON ) OPTWRK = N + OPTWRK IF ( WNTVA ) THEN CALL ZGEQRF(N,N/2,U,LDU,CDUMMY,CDUMMY,-1,IERR) LWRK_ZGEQRF = INT( CDUMMY(1) ) CALL ZGESVD( 'S', 'O', N/2,N/2, V,LDV, S, U,LDU, $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) LWRK_ZGESVD2 = INT( CDUMMY(1) ) CALL ZUNMQR( 'R', 'C', N, N, N/2, U, LDU, CDUMMY, $ V, LDV, CDUMMY, -1, IERR ) LWRK_ZUNMQR2 = INT( CDUMMY(1) ) OPTWRK2 = MAX( LWRK_ZGEQP3, N/2+LWRK_ZGEQRF, $ N/2+LWRK_ZGESVD2, N/2+LWRK_ZUNMQR2 ) IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON ) OPTWRK2 = N + OPTWRK2 OPTWRK = MAX( OPTWRK, OPTWRK2 ) END IF ELSE CALL ZGESVD( 'S', 'O', N, N, A, LDA, S, U, LDU, $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) LWRK_ZGESVD = INT( CDUMMY(1) ) OPTWRK = MAX(LWRK_ZGEQP3,LWRK_ZGESVD,LWRK_ZUNMQR) IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON ) OPTWRK = N + OPTWRK IF ( WNTVA ) THEN CALL ZGELQF(N/2,N,U,LDU,CDUMMY,CDUMMY,-1,IERR) LWRK_ZGELQF = INT( CDUMMY(1) ) CALL ZGESVD( 'S','O', N/2,N/2, V, LDV, S, U, LDU, $ V, LDV, CDUMMY, -1, RDUMMY, IERR ) LWRK_ZGESVD2 = INT( CDUMMY(1) ) CALL ZUNMLQ( 'R', 'N', N, N, N/2, U, LDU, CDUMMY, $ V, LDV, CDUMMY,-1,IERR ) LWRK_ZUNMLQ = INT( CDUMMY(1) ) OPTWRK2 = MAX( LWRK_ZGEQP3, N/2+LWRK_ZGELQF, $ N/2+LWRK_ZGESVD2, N/2+LWRK_ZUNMLQ ) IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON ) OPTWRK2 = N + OPTWRK2 OPTWRK = MAX( OPTWRK, OPTWRK2 ) END IF END IF END IF END IF * MINWRK = MAX( 2, MINWRK ) OPTWRK = MAX( 2, OPTWRK ) IF ( LCWORK .LT. MINWRK .AND. (.NOT.LQUERY) ) INFO = -19 * END IF * IF (INFO .EQ. 0 .AND. LRWORK .LT. RMINWRK .AND. .NOT. LQUERY) THEN INFO = -21 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'ZGESVDQ', -INFO ) RETURN ELSE IF ( LQUERY ) THEN * * Return optimal workspace * IWORK(1) = IMINWRK CWORK(1) = OPTWRK CWORK(2) = MINWRK RWORK(1) = RMINWRK RETURN END IF * * Quick return if the matrix is void. * IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) THEN * .. all output is void. RETURN END IF * BIG = DLAMCH('O') ASCALED = .FALSE. IF ( ROWPRM ) THEN * .. reordering the rows in decreasing sequence in the * ell-infinity norm - this enhances numerical robustness in * the case of differently scaled rows. DO 1904 p = 1, M * RWORK(p) = ABS( A(p,IZAMAX(N,A(p,1),LDA)) ) * [[ZLANGE will return NaN if an entry of the p-th row is Nan]] RWORK(p) = ZLANGE( 'M', 1, N, A(p,1), LDA, RDUMMY ) * .. check for NaN's and Inf's IF ( ( RWORK(p) .NE. RWORK(p) ) .OR. $ ( (RWORK(p)*ZERO) .NE. ZERO ) ) THEN INFO = -8 CALL XERBLA( 'ZGESVDQ', -INFO ) RETURN END IF 1904 CONTINUE DO 1952 p = 1, M - 1 q = IDAMAX( M-p+1, RWORK(p), 1 ) + p - 1 IWORK(N+p) = q IF ( p .NE. q ) THEN RTMP = RWORK(p) RWORK(p) = RWORK(q) RWORK(q) = RTMP END IF 1952 CONTINUE * IF ( RWORK(1) .EQ. ZERO ) THEN * Quick return: A is the M x N zero matrix. NUMRANK = 0 CALL DLASET( 'G', N, 1, ZERO, ZERO, S, N ) IF ( WNTUS ) CALL ZLASET('G', M, N, CZERO, CONE, U, LDU) IF ( WNTUA ) CALL ZLASET('G', M, M, CZERO, CONE, U, LDU) IF ( WNTVA ) CALL ZLASET('G', N, N, CZERO, CONE, V, LDV) IF ( WNTUF ) THEN CALL ZLASET( 'G', N, 1, CZERO, CZERO, CWORK, N ) CALL ZLASET( 'G', M, N, CZERO, CONE, U, LDU ) END IF DO 5001 p = 1, N IWORK(p) = p 5001 CONTINUE IF ( ROWPRM ) THEN DO 5002 p = N + 1, N + M - 1 IWORK(p) = p - N 5002 CONTINUE END IF IF ( CONDA ) RWORK(1) = -1 RWORK(2) = -1 RETURN END IF * IF ( RWORK(1) .GT. BIG / SQRT(DBLE(M)) ) THEN * .. to prevent overflow in the QR factorization, scale the * matrix by 1/sqrt(M) if too large entry detected CALL ZLASCL('G',0,0,SQRT(DBLE(M)),ONE, M,N, A,LDA, IERR) ASCALED = .TRUE. END IF CALL ZLASWP( N, A, LDA, 1, M-1, IWORK(N+1), 1 ) END IF * * .. At this stage, preemptive scaling is done only to avoid column * norms overflows during the QR factorization. The SVD procedure should * have its own scaling to save the singular values from overflows and * underflows. That depends on the SVD procedure. * IF ( .NOT.ROWPRM ) THEN RTMP = ZLANGE( 'M', M, N, A, LDA, RWORK ) IF ( ( RTMP .NE. RTMP ) .OR. $ ( (RTMP*ZERO) .NE. ZERO ) ) THEN INFO = -8 CALL XERBLA( 'ZGESVDQ', -INFO ) RETURN END IF IF ( RTMP .GT. BIG / SQRT(DBLE(M)) ) THEN * .. to prevent overflow in the QR factorization, scale the * matrix by 1/sqrt(M) if too large entry detected CALL ZLASCL('G',0,0, SQRT(DBLE(M)),ONE, M,N, A,LDA, IERR) ASCALED = .TRUE. END IF END IF * * .. QR factorization with column pivoting * * A * P = Q * [ R ] * [ 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), LCWORK-N, $ RWORK, IERR ) * * If the user requested accuracy level allows truncation in the * computed upper triangular factor, the matrix R is examined and, * if possible, replaced with its leading upper trapezoidal part. * EPSLN = DLAMCH('E') SFMIN = DLAMCH('S') * SMALL = SFMIN / EPSLN NR = N * IF ( ACCLA ) THEN * * Standard absolute error bound suffices. All sigma_i with * sigma_i < N*EPS*||A||_F are flushed to zero. This is an * aggressive enforcement of lower numerical rank by introducing a * backward error of the order of N*EPS*||A||_F. NR = 1 RTMP = SQRT(DBLE(N))*EPSLN DO 3001 p = 2, N IF ( ABS(A(p,p)) .LT. (RTMP*ABS(A(1,1))) ) GO TO 3002 NR = NR + 1 3001 CONTINUE 3002 CONTINUE * ELSEIF ( ACCLM ) THEN * .. similarly as above, only slightly more gentle (less aggressive). * Sudden drop on the diagonal of R is used as the criterion for being * close-to-rank-deficient. The threshold is set to EPSLN=DLAMCH('E'). * [[This can be made more flexible by replacing this hard-coded value * with a user specified threshold.]] Also, the values that underflow * will be truncated. NR = 1 DO 3401 p = 2, N IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR. $ ( ABS(A(p,p)) .LT. SFMIN ) ) GO TO 3402 NR = NR + 1 3401 CONTINUE 3402 CONTINUE * ELSE * .. RRQR not authorized to determine numerical rank except in the * obvious case of zero pivots. * .. inspect R for exact zeros on the diagonal; * R(i,i)=0 => R(i:N,i:N)=0. NR = 1 DO 3501 p = 2, N IF ( ABS(A(p,p)) .EQ. ZERO ) GO TO 3502 NR = NR + 1 3501 CONTINUE 3502 CONTINUE * IF ( CONDA ) THEN * Estimate the scaled condition number of A. Use the fact that it is * the same as the scaled condition number of R. * .. V is used as workspace CALL ZLACPY( 'U', N, N, A, LDA, V, LDV ) * Only the leading NR x NR submatrix of the triangular factor * is considered. Only if NR=N will this give a reliable error * bound. However, even for NR < N, this can be used on an * expert level and obtain useful information in the sense of * perturbation theory. DO 3053 p = 1, NR RTMP = DZNRM2( p, V(1,p), 1 ) CALL ZDSCAL( p, ONE/RTMP, V(1,p), 1 ) 3053 CONTINUE IF ( .NOT. ( LSVEC .OR. RSVEC ) ) THEN CALL ZPOCON( 'U', NR, V, LDV, ONE, RTMP, $ CWORK, RWORK, IERR ) ELSE CALL ZPOCON( 'U', NR, V, LDV, ONE, RTMP, $ CWORK(N+1), RWORK, IERR ) END IF SCONDA = ONE / SQRT(RTMP) * For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1), * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA * See the reference [1] for more details. END IF * ENDIF * IF ( WNTUR ) THEN N1 = NR ELSE IF ( WNTUS .OR. WNTUF) THEN N1 = N ELSE IF ( WNTUA ) THEN N1 = M END IF * IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN *....................................................................... * .. only the singular values are requested *....................................................................... IF ( RTRANS ) THEN * * .. compute the singular values of R**H = [A](1:NR,1:N)**H * .. set the lower triangle of [A] to [A](1:NR,1:N)**H and * the upper triangle of [A] to zero. DO 1146 p = 1, MIN( N, NR ) A(p,p) = CONJG(A(p,p)) DO 1147 q = p + 1, N A(q,p) = CONJG(A(p,q)) IF ( q .LE. NR ) A(p,q) = CZERO 1147 CONTINUE 1146 CONTINUE * CALL ZGESVD( 'N', 'N', N, NR, A, LDA, S, U, LDU, $ V, LDV, CWORK, LCWORK, RWORK, INFO ) * ELSE * * .. compute the singular values of R = [A](1:NR,1:N) * IF ( NR .GT. 1 ) $ CALL ZLASET( 'L', NR-1,NR-1, CZERO,CZERO, A(2,1), LDA ) CALL ZGESVD( 'N', 'N', NR, N, A, LDA, S, U, LDU, $ V, LDV, CWORK, LCWORK, RWORK, INFO ) * END IF * ELSE IF ( LSVEC .AND. ( .NOT. RSVEC) ) THEN *....................................................................... * .. the singular values and the left singular vectors requested *......................................................................."""""""" IF ( RTRANS ) THEN * .. apply ZGESVD to R**H * .. copy R**H into [U] and overwrite [U] with the right singular * vectors of R DO 1192 p = 1, NR DO 1193 q = p, N U(q,p) = CONJG(A(p,q)) 1193 CONTINUE 1192 CONTINUE IF ( NR .GT. 1 ) $ CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, U(1,2), LDU ) * .. the left singular vectors not computed, the NR right singular * vectors overwrite [U](1:NR,1:NR) as conjugate transposed. These * will be pre-multiplied by Q to build the left singular vectors of A. CALL ZGESVD( 'N', 'O', N, NR, U, LDU, S, U, LDU, $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO ) * DO 1119 p = 1, NR U(p,p) = CONJG(U(p,p)) DO 1120 q = p + 1, NR CTMP = CONJG(U(q,p)) U(q,p) = CONJG(U(p,q)) U(p,q) = CTMP 1120 CONTINUE 1119 CONTINUE * ELSE * .. apply ZGESVD to R * .. copy R into [U] and overwrite [U] with the left singular vectors CALL ZLACPY( 'U', NR, N, A, LDA, U, LDU ) IF ( NR .GT. 1 ) $ CALL ZLASET( 'L', NR-1, NR-1, CZERO, CZERO, U(2,1), LDU ) * .. the right singular vectors not computed, the NR left singular * vectors overwrite [U](1:NR,1:NR) CALL ZGESVD( 'O', 'N', NR, N, U, LDU, S, U, LDU, $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO ) * .. now [U](1:NR,1:NR) contains the NR left singular vectors of * R. These will be pre-multiplied by Q to build the left singular * vectors of A. END IF * * .. assemble the left singular vector matrix U of dimensions * (M x NR) or (M x N) or (M x M). IF ( ( NR .LT. M ) .AND. ( .NOT.WNTUF ) ) 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 * vectors matrix U. * IF ( .NOT.WNTUF ) $ CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U, $ LDU, CWORK(N+1), LCWORK-N, IERR ) IF ( ROWPRM .AND. .NOT.WNTUF ) $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 ) * ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN *....................................................................... * .. the singular values and the right singular vectors requested *....................................................................... IF ( RTRANS ) THEN * .. apply ZGESVD to R**H * .. copy R**H into V and overwrite V with the left singular vectors DO 1165 p = 1, NR DO 1166 q = p, N V(q,p) = CONJG(A(p,q)) 1166 CONTINUE 1165 CONTINUE IF ( NR .GT. 1 ) $ CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV ) * .. the left singular vectors of R**H overwrite V, the right singular * vectors not computed IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN CALL ZGESVD( 'O', 'N', N, NR, V, LDV, S, U, LDU, $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO ) * DO 1121 p = 1, NR V(p,p) = CONJG(V(p,p)) DO 1122 q = p + 1, NR CTMP = CONJG(V(q,p)) V(q,p) = CONJG(V(p,q)) V(p,q) = CTMP 1122 CONTINUE 1121 CONTINUE * IF ( NR .LT. N ) THEN DO 1103 p = 1, NR DO 1104 q = NR + 1, N V(p,q) = CONJG(V(q,p)) 1104 CONTINUE 1103 CONTINUE END IF CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK ) ELSE * .. need all N right singular vectors and NR < N * [!] This is simple implementation that augments [V](1:N,1:NR) * by padding a zero block. In the case NR << N, a more efficient * way is to first use the QR factorization. For more details * how to implement this, see the " FULL SVD " branch. CALL ZLASET('G', N, N-NR, CZERO, CZERO, V(1,NR+1), LDV) CALL ZGESVD( 'O', 'N', N, N, V, LDV, S, U, LDU, $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO ) * DO 1123 p = 1, N V(p,p) = CONJG(V(p,p)) DO 1124 q = p + 1, N CTMP = CONJG(V(q,p)) V(q,p) = CONJG(V(p,q)) V(p,q) = CTMP 1124 CONTINUE 1123 CONTINUE CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK ) END IF * ELSE * .. aply ZGESVD to R * .. copy R into V and overwrite V with the right singular vectors CALL ZLACPY( 'U', NR, N, A, LDA, V, LDV ) IF ( NR .GT. 1 ) $ CALL ZLASET( 'L', NR-1, NR-1, CZERO, CZERO, V(2,1), LDV ) * .. the right singular vectors overwrite V, the NR left singular * vectors stored in U(1:NR,1:NR) IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN CALL ZGESVD( 'N', 'O', NR, N, V, LDV, S, U, LDU, $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO ) CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK ) * .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H ELSE * .. need all N right singular vectors and NR < N * [!] This is simple implementation that augments [V](1:NR,1:N) * by padding a zero block. In the case NR << N, a more efficient * way is to first use the LQ factorization. For more details * how to implement this, see the " FULL SVD " branch. CALL ZLASET('G', N-NR, N, CZERO,CZERO, V(NR+1,1), LDV) CALL ZGESVD( 'N', 'O', N, N, V, LDV, S, U, LDU, $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO ) CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK ) END IF * .. now [V] contains the adjoint of the matrix of the right singular * vectors of A. END IF * ELSE *....................................................................... * .. FULL SVD requested *....................................................................... IF ( RTRANS ) THEN * * .. apply ZGESVD to R**H [[this option is left for R&D&T]] * IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN * .. copy R**H into [V] and overwrite [V] with the left singular * vectors of R**H DO 1168 p = 1, NR DO 1169 q = p, N V(q,p) = CONJG(A(p,q)) 1169 CONTINUE 1168 CONTINUE IF ( NR .GT. 1 ) $ CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV ) * * .. the left singular vectors of R**H overwrite [V], the NR right * singular vectors of R**H stored in [U](1:NR,1:NR) as conjugate * transposed CALL ZGESVD( 'O', 'A', N, NR, V, LDV, S, V, LDV, $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO ) * .. assemble V DO 1115 p = 1, NR V(p,p) = CONJG(V(p,p)) DO 1116 q = p + 1, NR CTMP = CONJG(V(q,p)) V(q,p) = CONJG(V(p,q)) V(p,q) = CTMP 1116 CONTINUE 1115 CONTINUE IF ( NR .LT. N ) THEN DO 1101 p = 1, NR DO 1102 q = NR+1, N V(p,q) = CONJG(V(q,p)) 1102 CONTINUE 1101 CONTINUE END IF CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK ) * DO 1117 p = 1, NR U(p,p) = CONJG(U(p,p)) DO 1118 q = p + 1, NR CTMP = CONJG(U(q,p)) U(q,p) = CONJG(U(p,q)) U(p,q) = CTMP 1118 CONTINUE 1117 CONTINUE * IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) 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 * ELSE * .. need all N right singular vectors and NR < N * .. copy R**H into [V] and overwrite [V] with the left singular * vectors of R**H * [[The optimal ratio N/NR for using QRF instead of padding * with zeros. Here hard coded to 2; it must be at least * two due to work space constraints.]] * OPTRATIO = ILAENV(6, 'ZGESVD', 'S' // 'O', NR,N,0,0) * OPTRATIO = MAX( OPTRATIO, 2 ) OPTRATIO = 2 IF ( OPTRATIO*NR .GT. N ) THEN DO 1198 p = 1, NR DO 1199 q = p, N V(q,p) = CONJG(A(p,q)) 1199 CONTINUE 1198 CONTINUE IF ( NR .GT. 1 ) $ CALL ZLASET('U',NR-1,NR-1, CZERO,CZERO, V(1,2),LDV) * CALL ZLASET('A',N,N-NR,CZERO,CZERO,V(1,NR+1),LDV) CALL ZGESVD( 'O', 'A', N, N, V, LDV, S, V, LDV, $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO ) * DO 1113 p = 1, N V(p,p) = CONJG(V(p,p)) DO 1114 q = p + 1, N CTMP = CONJG(V(q,p)) V(q,p) = CONJG(V(p,q)) V(p,q) = CTMP 1114 CONTINUE 1113 CONTINUE CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK ) * .. assemble the left singular vector matrix U of dimensions * (M x N1), i.e. (M x N) or (M x M). * DO 1111 p = 1, N U(p,p) = CONJG(U(p,p)) DO 1112 q = p + 1, N CTMP = CONJG(U(q,p)) U(q,p) = CONJG(U(p,q)) U(p,q) = CTMP 1112 CONTINUE 1111 CONTINUE * IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) 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 ELSE * .. copy R**H into [U] and overwrite [U] with the right * singular vectors of R DO 1196 p = 1, NR DO 1197 q = p, N U(q,NR+p) = CONJG(A(p,q)) 1197 CONTINUE 1196 CONTINUE IF ( NR .GT. 1 ) $ CALL ZLASET('U',NR-1,NR-1,CZERO,CZERO,U(1,NR+2),LDU) CALL ZGEQRF( N, NR, U(1,NR+1), LDU, CWORK(N+1), $ CWORK(N+NR+1), LCWORK-N-NR, IERR ) DO 1143 p = 1, NR DO 1144 q = 1, N V(q,p) = CONJG(U(p,NR+q)) 1144 CONTINUE 1143 CONTINUE CALL ZLASET('U',NR-1,NR-1,CZERO,CZERO,V(1,2),LDV) CALL ZGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU, $ V,LDV, CWORK(N+NR+1),LCWORK-N-NR,RWORK, INFO ) 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) CALL ZUNMQR('R','C', N, N, NR, U(1,NR+1), LDU, $ CWORK(N+1),V,LDV,CWORK(N+NR+1),LCWORK-N-NR,IERR) CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK ) * .. assemble the left singular vector matrix U of dimensions * (M x NR) or (M x N) or (M x M). IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) 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 END IF END IF * ELSE * * .. apply ZGESVD to R [[this is the recommended option]] * IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN * .. copy R into [V] and overwrite V with the right singular vectors CALL ZLACPY( 'U', NR, N, A, LDA, V, LDV ) IF ( NR .GT. 1 ) $ CALL ZLASET( 'L', NR-1,NR-1, CZERO,CZERO, V(2,1), LDV ) * .. the right singular vectors of R overwrite [V], the NR left * singular vectors of R stored in [U](1:NR,1:NR) CALL ZGESVD( 'S', 'O', NR, N, V, LDV, S, U, LDU, $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO ) CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK ) * .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H * .. assemble the left singular vector matrix U of dimensions * (M x NR) or (M x N) or (M x M). IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) 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 * ELSE * .. need all N right singular vectors and NR < N * .. the requested number of the left singular vectors * is then N1 (N or M) * [[The optimal ratio N/NR for using LQ instead of padding * with zeros. Here hard coded to 2; it must be at least * two due to work space constraints.]] * OPTRATIO = ILAENV(6, 'ZGESVD', 'S' // 'O', NR,N,0,0) * OPTRATIO = MAX( OPTRATIO, 2 ) OPTRATIO = 2 IF ( OPTRATIO * NR .GT. N ) THEN CALL ZLACPY( 'U', NR, N, A, LDA, V, LDV ) IF ( NR .GT. 1 ) $ CALL ZLASET('L', NR-1,NR-1, CZERO,CZERO, V(2,1),LDV) * .. the right singular vectors of R overwrite [V], the NR left * singular vectors of R stored in [U](1:NR,1:NR) CALL ZLASET('A', N-NR,N, CZERO,CZERO, V(NR+1,1),LDV) CALL ZGESVD( 'S', 'O', N, N, V, LDV, S, U, LDU, $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO ) CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK ) * .. now [V] contains the adjoint of the matrix of the right * singular vectors of A. The leading N left singular vectors * are in [U](1:N,1:N) * .. assemble the left singular vector matrix U of dimensions * (M x N1), i.e. (M x N) or (M x M). IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) 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 ELSE CALL ZLACPY( 'U', NR, N, A, LDA, U(NR+1,1), LDU ) IF ( NR .GT. 1 ) $ CALL ZLASET('L',NR-1,NR-1,CZERO,CZERO,U(NR+2,1),LDU) CALL ZGELQF( NR, N, U(NR+1,1), LDU, CWORK(N+1), $ CWORK(N+NR+1), LCWORK-N-NR, IERR ) CALL ZLACPY('L',NR,NR,U(NR+1,1),LDU,V,LDV) IF ( NR .GT. 1 ) $ CALL ZLASET('U',NR-1,NR-1,CZERO,CZERO,V(1,2),LDV) CALL ZGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU, $ V, LDV, CWORK(N+NR+1), LCWORK-N-NR, RWORK, INFO ) 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) CALL ZUNMLQ('R','N',N,N,NR,U(NR+1,1),LDU,CWORK(N+1), $ V, LDV, CWORK(N+NR+1),LCWORK-N-NR,IERR) CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK ) * .. assemble the left singular vector matrix U of dimensions * (M x NR) or (M x N) or (M x M). IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) 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 END IF END IF * .. end of the "R**H or R" branch END IF * * The Q matrix from the first QRF is built into the left singular * vectors matrix U. * IF ( .NOT. WNTUF ) $ CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U, $ LDU, CWORK(N+1), LCWORK-N, IERR ) IF ( ROWPRM .AND. .NOT.WNTUF ) $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 ) * * ... end of the "full SVD" branch END IF * * Check whether some singular values are returned as zeros, e.g. * due to underflow, and update the numerical rank. p = NR DO 4001 q = p, 1, -1 IF ( S(q) .GT. ZERO ) GO TO 4002 NR = NR - 1 4001 CONTINUE 4002 CONTINUE * * .. if numerical rank deficiency is detected, the truncated * singular values are set to zero. IF ( NR .LT. N ) CALL DLASET( 'G', N-NR,1, ZERO,ZERO, S(NR+1), N ) * .. undo scaling; this may cause overflow in the largest singular * values. IF ( ASCALED ) $ CALL DLASCL( 'G',0,0, ONE,SQRT(DBLE(M)), NR,1, S, N, IERR ) IF ( CONDA ) RWORK(1) = SCONDA RWORK(2) = p - NR * .. p-NR is the number of singular values that are computed as * exact zeros in ZGESVD() applied to the (possibly truncated) * full row rank triangular (trapezoidal) factor of A. NUMRANK = NR * RETURN * * End of ZGESVDQ * END