Annotation of rpl/lapack/lapack/zgesvdq.f, revision 1.1
1.1 ! bertrand 1: *> \brief <b> ZGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices</b>
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZGESVDQ + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvdq.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvdq.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvdq.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
! 22: * S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
! 23: * CWORK, LCWORK, RWORK, LRWORK, INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * IMPLICIT NONE
! 27: * CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV
! 28: * INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LCWORK, LRWORK,
! 29: * INFO
! 30: * ..
! 31: * .. Array Arguments ..
! 32: * COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * )
! 33: * DOUBLE PRECISION S( * ), RWORK( * )
! 34: * INTEGER IWORK( * )
! 35: * ..
! 36: *
! 37: *
! 38: *> \par Purpose:
! 39: * =============
! 40: *>
! 41: *> \verbatim
! 42: *>
! 43: * ZCGESVDQ computes the singular value decomposition (SVD) of a complex
! 44: *> M-by-N matrix A, where M >= N. The SVD of A is written as
! 45: *> [++] [xx] [x0] [xx]
! 46: *> A = U * SIGMA * V^*, [++] = [xx] * [ox] * [xx]
! 47: *> [++] [xx]
! 48: *> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
! 49: *> matrix, and V is an N-by-N unitary matrix. The diagonal elements
! 50: *> of SIGMA are the singular values of A. The columns of U and V are the
! 51: *> left and the right singular vectors of A, respectively.
! 52: *> \endverbatim
! 53: *
! 54: * Arguments
! 55: * =========
! 56: *
! 57: *> \param[in] JOBA
! 58: *> \verbatim
! 59: *> JOBA is CHARACTER*1
! 60: *> Specifies the level of accuracy in the computed SVD
! 61: *> = 'A' The requested accuracy corresponds to having the backward
! 62: *> error bounded by || delta A ||_F <= f(m,n) * EPS * || A ||_F,
! 63: *> where EPS = DLAMCH('Epsilon'). This authorises ZGESVDQ to
! 64: *> truncate the computed triangular factor in a rank revealing
! 65: *> QR factorization whenever the truncated part is below the
! 66: *> threshold of the order of EPS * ||A||_F. This is aggressive
! 67: *> truncation level.
! 68: *> = 'M' Similarly as with 'A', but the truncation is more gentle: it
! 69: *> is allowed only when there is a drop on the diagonal of the
! 70: *> triangular factor in the QR factorization. This is medium
! 71: *> truncation level.
! 72: *> = 'H' High accuracy requested. No numerical rank determination based
! 73: *> on the rank revealing QR factorization is attempted.
! 74: *> = 'E' Same as 'H', and in addition the condition number of column
! 75: *> scaled A is estimated and returned in RWORK(1).
! 76: *> N^(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <= N^(1/4)*RWORK(1)
! 77: *> \endverbatim
! 78: *>
! 79: *> \param[in] JOBP
! 80: *> \verbatim
! 81: *> JOBP is CHARACTER*1
! 82: *> = 'P' The rows of A are ordered in decreasing order with respect to
! 83: *> ||A(i,:)||_\infty. This enhances numerical accuracy at the cost
! 84: *> of extra data movement. Recommended for numerical robustness.
! 85: *> = 'N' No row pivoting.
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[in] JOBR
! 89: *> \verbatim
! 90: *> JOBR is CHARACTER*1
! 91: *> = 'T' After the initial pivoted QR factorization, ZGESVD is applied to
! 92: *> the adjoint R**H of the computed triangular factor R. This involves
! 93: *> some extra data movement (matrix transpositions). Useful for
! 94: *> experiments, research and development.
! 95: *> = 'N' The triangular factor R is given as input to CGESVD. This may be
! 96: *> preferred as it involves less data movement.
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] JOBU
! 100: *> \verbatim
! 101: *> JOBU is CHARACTER*1
! 102: *> = 'A' All M left singular vectors are computed and returned in the
! 103: *> matrix U. See the description of U.
! 104: *> = 'S' or 'U' N = min(M,N) left singular vectors are computed and returned
! 105: *> in the matrix U. See the description of U.
! 106: *> = 'R' Numerical rank NUMRANK is determined and only NUMRANK left singular
! 107: *> vectors are computed and returned in the matrix U.
! 108: *> = 'F' The N left singular vectors are returned in factored form as the
! 109: *> product of the Q factor from the initial QR factorization and the
! 110: *> N left singular vectors of (R**H , 0)**H. If row pivoting is used,
! 111: *> then the necessary information on the row pivoting is stored in
! 112: *> IWORK(N+1:N+M-1).
! 113: *> = 'N' The left singular vectors are not computed.
! 114: *> \endverbatim
! 115: *>
! 116: *> \param[in] JOBV
! 117: *> \verbatim
! 118: *> JOBV is CHARACTER*1
! 119: *> = 'A', 'V' All N right singular vectors are computed and returned in
! 120: *> the matrix V.
! 121: *> = 'R' Numerical rank NUMRANK is determined and only NUMRANK right singular
! 122: *> vectors are computed and returned in the matrix V. This option is
! 123: *> allowed only if JOBU = 'R' or JOBU = 'N'; otherwise it is illegal.
! 124: *> = 'N' The right singular vectors are not computed.
! 125: *> \endverbatim
! 126: *>
! 127: *> \param[in] M
! 128: *> \verbatim
! 129: *> M is INTEGER
! 130: *> The number of rows of the input matrix A. M >= 0.
! 131: *> \endverbatim
! 132: *>
! 133: *> \param[in] N
! 134: *> \verbatim
! 135: *> N is INTEGER
! 136: *> The number of columns of the input matrix A. M >= N >= 0.
! 137: *> \endverbatim
! 138: *>
! 139: *> \param[in,out] A
! 140: *> \verbatim
! 141: *> A is COMPLEX*16 array of dimensions LDA x N
! 142: *> On entry, the input matrix A.
! 143: *> On exit, if JOBU .NE. 'N' or JOBV .NE. 'N', the lower triangle of A contains
! 144: *> the Householder vectors as stored by ZGEQP3. If JOBU = 'F', these Householder
! 145: *> vectors together with CWORK(1:N) can be used to restore the Q factors from
! 146: *> the initial pivoted QR factorization of A. See the description of U.
! 147: *> \endverbatim
! 148: *>
! 149: *> \param[in] LDA
! 150: *> \verbatim
! 151: *> LDA is INTEGER.
! 152: *> The leading dimension of the array A. LDA >= max(1,M).
! 153: *> \endverbatim
! 154: *>
! 155: *> \param[out] S
! 156: *> \verbatim
! 157: *> S is DOUBLE PRECISION array of dimension N.
! 158: *> The singular values of A, ordered so that S(i) >= S(i+1).
! 159: *> \endverbatim
! 160: *>
! 161: *> \param[out] U
! 162: *> \verbatim
! 163: *> U is COMPLEX*16 array, dimension
! 164: *> LDU x M if JOBU = 'A'; see the description of LDU. In this case,
! 165: *> on exit, U contains the M left singular vectors.
! 166: *> LDU x N if JOBU = 'S', 'U', 'R' ; see the description of LDU. In this
! 167: *> case, U contains the leading N or the leading NUMRANK left singular vectors.
! 168: *> LDU x N if JOBU = 'F' ; see the description of LDU. In this case U
! 169: *> contains N x N unitary matrix that can be used to form the left
! 170: *> singular vectors.
! 171: *> If JOBU = 'N', U is not referenced.
! 172: *> \endverbatim
! 173: *>
! 174: *> \param[in] LDU
! 175: *> \verbatim
! 176: *> LDU is INTEGER.
! 177: *> The leading dimension of the array U.
! 178: *> If JOBU = 'A', 'S', 'U', 'R', LDU >= max(1,M).
! 179: *> If JOBU = 'F', LDU >= max(1,N).
! 180: *> Otherwise, LDU >= 1.
! 181: *> \endverbatim
! 182: *>
! 183: *> \param[out] V
! 184: *> \verbatim
! 185: *> V is COMPLEX*16 array, dimension
! 186: *> LDV x N if JOBV = 'A', 'V', 'R' or if JOBA = 'E' .
! 187: *> If JOBV = 'A', or 'V', V contains the N-by-N unitary matrix V**H;
! 188: *> If JOBV = 'R', V contains the first NUMRANK rows of V**H (the right
! 189: *> singular vectors, stored rowwise, of the NUMRANK largest singular values).
! 190: *> If JOBV = 'N' and JOBA = 'E', V is used as a workspace.
! 191: *> If JOBV = 'N', and JOBA.NE.'E', V is not referenced.
! 192: *> \endverbatim
! 193: *>
! 194: *> \param[in] LDV
! 195: *> \verbatim
! 196: *> LDV is INTEGER
! 197: *> The leading dimension of the array V.
! 198: *> If JOBV = 'A', 'V', 'R', or JOBA = 'E', LDV >= max(1,N).
! 199: *> Otherwise, LDV >= 1.
! 200: *> \endverbatim
! 201: *>
! 202: *> \param[out] NUMRANK
! 203: *> \verbatim
! 204: *> NUMRANK is INTEGER
! 205: *> NUMRANK is the numerical rank first determined after the rank
! 206: *> revealing QR factorization, following the strategy specified by the
! 207: *> value of JOBA. If JOBV = 'R' and JOBU = 'R', only NUMRANK
! 208: *> leading singular values and vectors are then requested in the call
! 209: *> of CGESVD. The final value of NUMRANK might be further reduced if
! 210: *> some singular values are computed as zeros.
! 211: *> \endverbatim
! 212: *>
! 213: *> \param[out] IWORK
! 214: *> \verbatim
! 215: *> IWORK is INTEGER array, dimension (max(1, LIWORK)).
! 216: *> On exit, IWORK(1:N) contains column pivoting permutation of the
! 217: *> rank revealing QR factorization.
! 218: *> If JOBP = 'P', IWORK(N+1:N+M-1) contains the indices of the sequence
! 219: *> of row swaps used in row pivoting. These can be used to restore the
! 220: *> left singular vectors in the case JOBU = 'F'.
! 221: *
! 222: *> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
! 223: *> LIWORK(1) returns the minimal LIWORK.
! 224: *> \endverbatim
! 225: *>
! 226: *> \param[in] LIWORK
! 227: *> \verbatim
! 228: *> LIWORK is INTEGER
! 229: *> The dimension of the array IWORK.
! 230: *> LIWORK >= N + M - 1, if JOBP = 'P';
! 231: *> LIWORK >= N if JOBP = 'N'.
! 232: *>
! 233: *> If LIWORK = -1, then a workspace query is assumed; the routine
! 234: *> only calculates and returns the optimal and minimal sizes
! 235: *> for the CWORK, IWORK, and RWORK arrays, and no error
! 236: *> message related to LCWORK is issued by XERBLA.
! 237: *> \endverbatim
! 238: *>
! 239: *> \param[out] CWORK
! 240: *> \verbatim
! 241: *> CWORK is COMPLEX*12 array, dimension (max(2, LCWORK)), used as a workspace.
! 242: *> On exit, if, on entry, LCWORK.NE.-1, CWORK(1:N) contains parameters
! 243: *> needed to recover the Q factor from the QR factorization computed by
! 244: *> ZGEQP3.
! 245: *
! 246: *> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
! 247: *> CWORK(1) returns the optimal LCWORK, and
! 248: *> CWORK(2) returns the minimal LCWORK.
! 249: *> \endverbatim
! 250: *>
! 251: *> \param[in,out] LCWORK
! 252: *> \verbatim
! 253: *> LCWORK is INTEGER
! 254: *> The dimension of the array CWORK. It is determined as follows:
! 255: *> Let LWQP3 = N+1, LWCON = 2*N, and let
! 256: *> LWUNQ = { MAX( N, 1 ), if JOBU = 'R', 'S', or 'U'
! 257: *> { MAX( M, 1 ), if JOBU = 'A'
! 258: *> LWSVD = MAX( 3*N, 1 )
! 259: *> LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 3*(N/2), 1 ), LWUNLQ = MAX( N, 1 ),
! 260: *> LWQRF = MAX( N/2, 1 ), LWUNQ2 = MAX( N, 1 )
! 261: *> Then the minimal value of LCWORK is:
! 262: *> = MAX( N + LWQP3, LWSVD ) if only the singular values are needed;
! 263: *> = MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values are needed,
! 264: *> and a scaled condition estimate requested;
! 265: *>
! 266: *> = N + MAX( LWQP3, LWSVD, LWUNQ ) if the singular values and the left
! 267: *> singular vectors are requested;
! 268: *> = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the singular values and the left
! 269: *> singular vectors are requested, and also
! 270: *> a scaled condition estimate requested;
! 271: *>
! 272: *> = N + MAX( LWQP3, LWSVD ) if the singular values and the right
! 273: *> singular vectors are requested;
! 274: *> = N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and the right
! 275: *> singular vectors are requested, and also
! 276: *> a scaled condition etimate requested;
! 277: *>
! 278: *> = N + MAX( LWQP3, LWSVD, LWUNQ ) if the full SVD is requested with JOBV = 'R';
! 279: *> independent of JOBR;
! 280: *> = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the full SVD is requested,
! 281: *> JOBV = 'R' and, also a scaled condition
! 282: *> estimate requested; independent of JOBR;
! 283: *> = MAX( N + MAX( LWQP3, LWSVD, LWUNQ ),
! 284: *> N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ) ) if the
! 285: *> full SVD is requested with JOBV = 'A' or 'V', and
! 286: *> JOBR ='N'
! 287: *> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
! 288: *> N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ ) )
! 289: *> if the full SVD is requested with JOBV = 'A' or 'V', and
! 290: *> JOBR ='N', and also a scaled condition number estimate
! 291: *> requested.
! 292: *> = MAX( N + MAX( LWQP3, LWSVD, LWUNQ ),
! 293: *> N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) ) if the
! 294: *> full SVD is requested with JOBV = 'A', 'V', and JOBR ='T'
! 295: *> = MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
! 296: *> N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) )
! 297: *> if the full SVD is requested with JOBV = 'A', 'V' and
! 298: *> JOBR ='T', and also a scaled condition number estimate
! 299: *> requested.
! 300: *> Finally, LCWORK must be at least two: LCWORK = MAX( 2, LCWORK ).
! 301: *>
! 302: *> If LCWORK = -1, then a workspace query is assumed; the routine
! 303: *> only calculates and returns the optimal and minimal sizes
! 304: *> for the CWORK, IWORK, and RWORK arrays, and no error
! 305: *> message related to LCWORK is issued by XERBLA.
! 306: *> \endverbatim
! 307: *>
! 308: *> \param[out] RWORK
! 309: *> \verbatim
! 310: *> RWORK is DOUBLE PRECISION array, dimension (max(1, LRWORK)).
! 311: *> On exit,
! 312: *> 1. If JOBA = 'E', RWORK(1) contains an estimate of the condition
! 313: *> number of column scaled A. If A = C * D where D is diagonal and C
! 314: *> has unit columns in the Euclidean norm, then, assuming full column rank,
! 315: *> N^(-1/4) * RWORK(1) <= ||pinv(C)||_2 <= N^(1/4) * RWORK(1).
! 316: *> Otherwise, RWORK(1) = -1.
! 317: *> 2. RWORK(2) contains the number of singular values computed as
! 318: *> exact zeros in ZGESVD applied to the upper triangular or trapeziodal
! 319: *> R (from the initial QR factorization). In case of early exit (no call to
! 320: *> ZGESVD, such as in the case of zero matrix) RWORK(2) = -1.
! 321: *
! 322: *> If LIWORK, LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
! 323: *> RWORK(1) returns the minimal LRWORK.
! 324: *> \endverbatim
! 325: *>
! 326: *> \param[in] LRWORK
! 327: *> \verbatim
! 328: *> LRWORK is INTEGER.
! 329: *> The dimension of the array RWORK.
! 330: *> If JOBP ='P', then LRWORK >= MAX(2, M, 5*N);
! 331: *> Otherwise, LRWORK >= MAX(2, 5*N).
! 332: *
! 333: *> If LRWORK = -1, then a workspace query is assumed; the routine
! 334: *> only calculates and returns the optimal and minimal sizes
! 335: *> for the CWORK, IWORK, and RWORK arrays, and no error
! 336: *> message related to LCWORK is issued by XERBLA.
! 337: *> \endverbatim
! 338: *>
! 339: *> \param[out] INFO
! 340: *> \verbatim
! 341: *> INFO is INTEGER
! 342: *> = 0: successful exit.
! 343: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 344: *> > 0: if ZBDSQR did not converge, INFO specifies how many superdiagonals
! 345: *> of an intermediate bidiagonal form B (computed in ZGESVD) did not
! 346: *> converge to zero.
! 347: *> \endverbatim
! 348: *
! 349: *> \par Further Details:
! 350: * ========================
! 351: *>
! 352: *> \verbatim
! 353: *>
! 354: *> 1. The data movement (matrix transpose) is coded using simple nested
! 355: *> DO-loops because BLAS and LAPACK do not provide corresponding subroutines.
! 356: *> Those DO-loops are easily identified in this source code - by the CONTINUE
! 357: *> statements labeled with 11**. In an optimized version of this code, the
! 358: *> nested DO loops should be replaced with calls to an optimized subroutine.
! 359: *> 2. This code scales A by 1/SQRT(M) if the largest ABS(A(i,j)) could cause
! 360: *> column norm overflow. This is the minial precaution and it is left to the
! 361: *> SVD routine (CGESVD) to do its own preemptive scaling if potential over-
! 362: *> or underflows are detected. To avoid repeated scanning of the array A,
! 363: *> an optimal implementation would do all necessary scaling before calling
! 364: *> CGESVD and the scaling in CGESVD can be switched off.
! 365: *> 3. Other comments related to code optimization are given in comments in the
! 366: *> code, enlosed in [[double brackets]].
! 367: *> \endverbatim
! 368: *
! 369: *> \par Bugs, examples and comments
! 370: * ===========================
! 371: *
! 372: *> \verbatim
! 373: *> Please report all bugs and send interesting examples and/or comments to
! 374: *> drmac@math.hr. Thank you.
! 375: *> \endverbatim
! 376: *
! 377: *> \par References
! 378: * ===============
! 379: *
! 380: *> \verbatim
! 381: *> [1] Zlatko Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
! 382: *> Computing the SVD with High Accuracy. ACM Trans. Math. Softw.
! 383: *> 44(1): 11:1-11:30 (2017)
! 384: *>
! 385: *> SIGMA library, xGESVDQ section updated February 2016.
! 386: *> Developed and coded by Zlatko Drmac, Department of Mathematics
! 387: *> University of Zagreb, Croatia, drmac@math.hr
! 388: *> \endverbatim
! 389: *
! 390: *
! 391: *> \par Contributors:
! 392: * ==================
! 393: *>
! 394: *> \verbatim
! 395: *> Developed and coded by Zlatko Drmac, Department of Mathematics
! 396: *> University of Zagreb, Croatia, drmac@math.hr
! 397: *> \endverbatim
! 398: *
! 399: * Authors:
! 400: * ========
! 401: *
! 402: *> \author Univ. of Tennessee
! 403: *> \author Univ. of California Berkeley
! 404: *> \author Univ. of Colorado Denver
! 405: *> \author NAG Ltd.
! 406: *
! 407: *> \date November 2018
! 408: *
! 409: *> \ingroup complex16GEsing
! 410: *
! 411: * =====================================================================
! 412: SUBROUTINE ZGESVDQ( JOBA, JOBP, JOBR, JOBU, JOBV, M, N, A, LDA,
! 413: $ S, U, LDU, V, LDV, NUMRANK, IWORK, LIWORK,
! 414: $ CWORK, LCWORK, RWORK, LRWORK, INFO )
! 415: * .. Scalar Arguments ..
! 416: IMPLICIT NONE
! 417: CHARACTER JOBA, JOBP, JOBR, JOBU, JOBV
! 418: INTEGER M, N, LDA, LDU, LDV, NUMRANK, LIWORK, LCWORK, LRWORK,
! 419: $ INFO
! 420: * ..
! 421: * .. Array Arguments ..
! 422: COMPLEX*16 A( LDA, * ), U( LDU, * ), V( LDV, * ), CWORK( * )
! 423: DOUBLE PRECISION S( * ), RWORK( * )
! 424: INTEGER IWORK( * )
! 425: *
! 426: * =====================================================================
! 427: *
! 428: * .. Parameters ..
! 429: DOUBLE PRECISION ZERO, ONE
! 430: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 431: COMPLEX*16 CZERO, CONE
! 432: PARAMETER ( CZERO = (0.0D0,0.0D0), CONE = (1.0D0,0.0D0) )
! 433: * ..
! 434: * .. Local Scalars ..
! 435: INTEGER IERR, NR, N1, OPTRATIO, p, q
! 436: INTEGER LWCON, LWQP3, LWRK_ZGELQF, LWRK_ZGESVD, LWRK_ZGESVD2,
! 437: $ LWRK_ZGEQP3, LWRK_ZGEQRF, LWRK_ZUNMLQ, LWRK_ZUNMQR,
! 438: $ LWRK_ZUNMQR2, LWLQF, LWQRF, LWSVD, LWSVD2, LWUNQ,
! 439: $ LWUNQ2, LWUNLQ, MINWRK, MINWRK2, OPTWRK, OPTWRK2,
! 440: $ IMINWRK, RMINWRK
! 441: LOGICAL ACCLA, ACCLM, ACCLH, ASCALED, CONDA, DNTWU, DNTWV,
! 442: $ LQUERY, LSVC0, LSVEC, ROWPRM, RSVEC, RTRANS, WNTUA,
! 443: $ WNTUF, WNTUR, WNTUS, WNTVA, WNTVR
! 444: DOUBLE PRECISION BIG, EPSLN, RTMP, SCONDA, SFMIN
! 445: COMPLEX*16 CTMP
! 446: * ..
! 447: * .. Local Arrays
! 448: COMPLEX*16 CDUMMY(1)
! 449: DOUBLE PRECISION RDUMMY(1)
! 450: * ..
! 451: * .. External Subroutines (BLAS, LAPACK)
! 452: EXTERNAL ZGELQF, ZGEQP3, ZGEQRF, ZGESVD, ZLACPY, ZLAPMT,
! 453: $ ZLASCL, ZLASET, ZLASWP, ZDSCAL, DLASET, DLASCL,
! 454: $ ZPOCON, ZUNMLQ, ZUNMQR, XERBLA
! 455: * ..
! 456: * .. External Functions (BLAS, LAPACK)
! 457: LOGICAL LSAME
! 458: INTEGER IDAMAX
! 459: DOUBLE PRECISION ZLANGE, DZNRM2, DLAMCH
! 460: EXTERNAL LSAME, ZLANGE, IDAMAX, DZNRM2, DLAMCH
! 461: * ..
! 462: * .. Intrinsic Functions ..
! 463: INTRINSIC ABS, CONJG, MAX, MIN, DBLE, SQRT
! 464: * ..
! 465: * .. Executable Statements ..
! 466: *
! 467: * Test the input arguments
! 468: *
! 469: WNTUS = LSAME( JOBU, 'S' ) .OR. LSAME( JOBU, 'U' )
! 470: WNTUR = LSAME( JOBU, 'R' )
! 471: WNTUA = LSAME( JOBU, 'A' )
! 472: WNTUF = LSAME( JOBU, 'F' )
! 473: LSVC0 = WNTUS .OR. WNTUR .OR. WNTUA
! 474: LSVEC = LSVC0 .OR. WNTUF
! 475: DNTWU = LSAME( JOBU, 'N' )
! 476: *
! 477: WNTVR = LSAME( JOBV, 'R' )
! 478: WNTVA = LSAME( JOBV, 'A' ) .OR. LSAME( JOBV, 'V' )
! 479: RSVEC = WNTVR .OR. WNTVA
! 480: DNTWV = LSAME( JOBV, 'N' )
! 481: *
! 482: ACCLA = LSAME( JOBA, 'A' )
! 483: ACCLM = LSAME( JOBA, 'M' )
! 484: CONDA = LSAME( JOBA, 'E' )
! 485: ACCLH = LSAME( JOBA, 'H' ) .OR. CONDA
! 486: *
! 487: ROWPRM = LSAME( JOBP, 'P' )
! 488: RTRANS = LSAME( JOBR, 'T' )
! 489: *
! 490: IF ( ROWPRM ) THEN
! 491: IMINWRK = MAX( 1, N + M - 1 )
! 492: RMINWRK = MAX( 2, M, 5*N )
! 493: ELSE
! 494: IMINWRK = MAX( 1, N )
! 495: RMINWRK = MAX( 2, 5*N )
! 496: END IF
! 497: LQUERY = (LIWORK .EQ. -1 .OR. LCWORK .EQ. -1 .OR. LRWORK .EQ. -1)
! 498: INFO = 0
! 499: IF ( .NOT. ( ACCLA .OR. ACCLM .OR. ACCLH ) ) THEN
! 500: INFO = -1
! 501: ELSE IF ( .NOT.( ROWPRM .OR. LSAME( JOBP, 'N' ) ) ) THEN
! 502: INFO = -2
! 503: ELSE IF ( .NOT.( RTRANS .OR. LSAME( JOBR, 'N' ) ) ) THEN
! 504: INFO = -3
! 505: ELSE IF ( .NOT.( LSVEC .OR. DNTWU ) ) THEN
! 506: INFO = -4
! 507: ELSE IF ( WNTUR .AND. WNTVA ) THEN
! 508: INFO = -5
! 509: ELSE IF ( .NOT.( RSVEC .OR. DNTWV )) THEN
! 510: INFO = -5
! 511: ELSE IF ( M.LT.0 ) THEN
! 512: INFO = -6
! 513: ELSE IF ( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
! 514: INFO = -7
! 515: ELSE IF ( LDA.LT.MAX( 1, M ) ) THEN
! 516: INFO = -9
! 517: ELSE IF ( LDU.LT.1 .OR. ( LSVC0 .AND. LDU.LT.M ) .OR.
! 518: $ ( WNTUF .AND. LDU.LT.N ) ) THEN
! 519: INFO = -12
! 520: ELSE IF ( LDV.LT.1 .OR. ( RSVEC .AND. LDV.LT.N ) .OR.
! 521: $ ( CONDA .AND. LDV.LT.N ) ) THEN
! 522: INFO = -14
! 523: ELSE IF ( LIWORK .LT. IMINWRK .AND. .NOT. LQUERY ) THEN
! 524: INFO = -17
! 525: END IF
! 526: *
! 527: *
! 528: IF ( INFO .EQ. 0 ) THEN
! 529: * .. compute the minimal and the optimal workspace lengths
! 530: * [[The expressions for computing the minimal and the optimal
! 531: * values of LCWORK are written with a lot of redundancy and
! 532: * can be simplified. However, this detailed form is easier for
! 533: * maintenance and modifications of the code.]]
! 534: *
! 535: * .. minimal workspace length for ZGEQP3 of an M x N matrix
! 536: LWQP3 = N+1
! 537: * .. minimal workspace length for ZUNMQR to build left singular vectors
! 538: IF ( WNTUS .OR. WNTUR ) THEN
! 539: LWUNQ = MAX( N , 1 )
! 540: ELSE IF ( WNTUA ) THEN
! 541: LWUNQ = MAX( M , 1 )
! 542: END IF
! 543: * .. minimal workspace length for ZPOCON of an N x N matrix
! 544: LWCON = 2 * N
! 545: * .. ZGESVD of an N x N matrix
! 546: LWSVD = MAX( 3 * N, 1 )
! 547: IF ( LQUERY ) THEN
! 548: CALL ZGEQP3( M, N, A, LDA, IWORK, CDUMMY, CDUMMY, -1,
! 549: $ RDUMMY, IERR )
! 550: LWRK_ZGEQP3 = INT( CDUMMY(1) )
! 551: IF ( WNTUS .OR. WNTUR ) THEN
! 552: CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
! 553: $ LDU, CDUMMY, -1, IERR )
! 554: LWRK_ZUNMQR = INT( CDUMMY(1) )
! 555: ELSE IF ( WNTUA ) THEN
! 556: CALL ZUNMQR( 'L', 'N', M, M, N, A, LDA, CDUMMY, U,
! 557: $ LDU, CDUMMY, -1, IERR )
! 558: LWRK_ZUNMQR = INT( CDUMMY(1) )
! 559: ELSE
! 560: LWRK_ZUNMQR = 0
! 561: END IF
! 562: END IF
! 563: MINWRK = 2
! 564: OPTWRK = 2
! 565: IF ( .NOT. (LSVEC .OR. RSVEC ) ) THEN
! 566: * .. minimal and optimal sizes of the complex workspace if
! 567: * only the singular values are requested
! 568: IF ( CONDA ) THEN
! 569: MINWRK = MAX( N+LWQP3, LWCON, LWSVD )
! 570: ELSE
! 571: MINWRK = MAX( N+LWQP3, LWSVD )
! 572: END IF
! 573: IF ( LQUERY ) THEN
! 574: CALL ZGESVD( 'N', 'N', N, N, A, LDA, S, U, LDU,
! 575: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
! 576: LWRK_ZGESVD = INT( CDUMMY(1) )
! 577: IF ( CONDA ) THEN
! 578: OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON, LWRK_ZGESVD )
! 579: ELSE
! 580: OPTWRK = MAX( N+LWRK_ZGEQP3, LWRK_ZGESVD )
! 581: END IF
! 582: END IF
! 583: ELSE IF ( LSVEC .AND. (.NOT.RSVEC) ) THEN
! 584: * .. minimal and optimal sizes of the complex workspace if the
! 585: * singular values and the left singular vectors are requested
! 586: IF ( CONDA ) THEN
! 587: MINWRK = N + MAX( LWQP3, LWCON, LWSVD, LWUNQ )
! 588: ELSE
! 589: MINWRK = N + MAX( LWQP3, LWSVD, LWUNQ )
! 590: END IF
! 591: IF ( LQUERY ) THEN
! 592: IF ( RTRANS ) THEN
! 593: CALL ZGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU,
! 594: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
! 595: ELSE
! 596: CALL ZGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU,
! 597: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
! 598: END IF
! 599: LWRK_ZGESVD = INT( CDUMMY(1) )
! 600: IF ( CONDA ) THEN
! 601: OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, LWRK_ZGESVD,
! 602: $ LWRK_ZUNMQR )
! 603: ELSE
! 604: OPTWRK = N + MAX( LWRK_ZGEQP3, LWRK_ZGESVD,
! 605: $ LWRK_ZUNMQR )
! 606: END IF
! 607: END IF
! 608: ELSE IF ( RSVEC .AND. (.NOT.LSVEC) ) THEN
! 609: * .. minimal and optimal sizes of the complex workspace if the
! 610: * singular values and the right singular vectors are requested
! 611: IF ( CONDA ) THEN
! 612: MINWRK = N + MAX( LWQP3, LWCON, LWSVD )
! 613: ELSE
! 614: MINWRK = N + MAX( LWQP3, LWSVD )
! 615: END IF
! 616: IF ( LQUERY ) THEN
! 617: IF ( RTRANS ) THEN
! 618: CALL ZGESVD( 'O', 'N', N, N, A, LDA, S, U, LDU,
! 619: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
! 620: ELSE
! 621: CALL ZGESVD( 'N', 'O', N, N, A, LDA, S, U, LDU,
! 622: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
! 623: END IF
! 624: LWRK_ZGESVD = INT( CDUMMY(1) )
! 625: IF ( CONDA ) THEN
! 626: OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, LWRK_ZGESVD )
! 627: ELSE
! 628: OPTWRK = N + MAX( LWRK_ZGEQP3, LWRK_ZGESVD )
! 629: END IF
! 630: END IF
! 631: ELSE
! 632: * .. minimal and optimal sizes of the complex workspace if the
! 633: * full SVD is requested
! 634: IF ( RTRANS ) THEN
! 635: MINWRK = MAX( LWQP3, LWSVD, LWUNQ )
! 636: IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON )
! 637: MINWRK = MINWRK + N
! 638: IF ( WNTVA ) THEN
! 639: * .. minimal workspace length for N x N/2 ZGEQRF
! 640: LWQRF = MAX( N/2, 1 )
! 641: * .. minimal workspace lengt for N/2 x N/2 ZGESVD
! 642: LWSVD2 = MAX( 3 * (N/2), 1 )
! 643: LWUNQ2 = MAX( N, 1 )
! 644: MINWRK2 = MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2,
! 645: $ N/2+LWUNQ2, LWUNQ )
! 646: IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON )
! 647: MINWRK2 = N + MINWRK2
! 648: MINWRK = MAX( MINWRK, MINWRK2 )
! 649: END IF
! 650: ELSE
! 651: MINWRK = MAX( LWQP3, LWSVD, LWUNQ )
! 652: IF ( CONDA ) MINWRK = MAX( MINWRK, LWCON )
! 653: MINWRK = MINWRK + N
! 654: IF ( WNTVA ) THEN
! 655: * .. minimal workspace length for N/2 x N ZGELQF
! 656: LWLQF = MAX( N/2, 1 )
! 657: LWSVD2 = MAX( 3 * (N/2), 1 )
! 658: LWUNLQ = MAX( N , 1 )
! 659: MINWRK2 = MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2,
! 660: $ N/2+LWUNLQ, LWUNQ )
! 661: IF ( CONDA ) MINWRK2 = MAX( MINWRK2, LWCON )
! 662: MINWRK2 = N + MINWRK2
! 663: MINWRK = MAX( MINWRK, MINWRK2 )
! 664: END IF
! 665: END IF
! 666: IF ( LQUERY ) THEN
! 667: IF ( RTRANS ) THEN
! 668: CALL ZGESVD( 'O', 'A', N, N, A, LDA, S, U, LDU,
! 669: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
! 670: LWRK_ZGESVD = INT( CDUMMY(1) )
! 671: OPTWRK = MAX(LWRK_ZGEQP3,LWRK_ZGESVD,LWRK_ZUNMQR)
! 672: IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON )
! 673: OPTWRK = N + OPTWRK
! 674: IF ( WNTVA ) THEN
! 675: CALL ZGEQRF(N,N/2,U,LDU,CDUMMY,CDUMMY,-1,IERR)
! 676: LWRK_ZGEQRF = INT( CDUMMY(1) )
! 677: CALL ZGESVD( 'S', 'O', N/2,N/2, V,LDV, S, U,LDU,
! 678: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
! 679: LWRK_ZGESVD2 = INT( CDUMMY(1) )
! 680: CALL ZUNMQR( 'R', 'C', N, N, N/2, U, LDU, CDUMMY,
! 681: $ V, LDV, CDUMMY, -1, IERR )
! 682: LWRK_ZUNMQR2 = INT( CDUMMY(1) )
! 683: OPTWRK2 = MAX( LWRK_ZGEQP3, N/2+LWRK_ZGEQRF,
! 684: $ N/2+LWRK_ZGESVD2, N/2+LWRK_ZUNMQR2 )
! 685: IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON )
! 686: OPTWRK2 = N + OPTWRK2
! 687: OPTWRK = MAX( OPTWRK, OPTWRK2 )
! 688: END IF
! 689: ELSE
! 690: CALL ZGESVD( 'S', 'O', N, N, A, LDA, S, U, LDU,
! 691: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
! 692: LWRK_ZGESVD = INT( CDUMMY(1) )
! 693: OPTWRK = MAX(LWRK_ZGEQP3,LWRK_ZGESVD,LWRK_ZUNMQR)
! 694: IF ( CONDA ) OPTWRK = MAX( OPTWRK, LWCON )
! 695: OPTWRK = N + OPTWRK
! 696: IF ( WNTVA ) THEN
! 697: CALL ZGELQF(N/2,N,U,LDU,CDUMMY,CDUMMY,-1,IERR)
! 698: LWRK_ZGELQF = INT( CDUMMY(1) )
! 699: CALL ZGESVD( 'S','O', N/2,N/2, V, LDV, S, U, LDU,
! 700: $ V, LDV, CDUMMY, -1, RDUMMY, IERR )
! 701: LWRK_ZGESVD2 = INT( CDUMMY(1) )
! 702: CALL ZUNMLQ( 'R', 'N', N, N, N/2, U, LDU, CDUMMY,
! 703: $ V, LDV, CDUMMY,-1,IERR )
! 704: LWRK_ZUNMLQ = INT( CDUMMY(1) )
! 705: OPTWRK2 = MAX( LWRK_ZGEQP3, N/2+LWRK_ZGELQF,
! 706: $ N/2+LWRK_ZGESVD2, N/2+LWRK_ZUNMLQ )
! 707: IF ( CONDA ) OPTWRK2 = MAX( OPTWRK2, LWCON )
! 708: OPTWRK2 = N + OPTWRK2
! 709: OPTWRK = MAX( OPTWRK, OPTWRK2 )
! 710: END IF
! 711: END IF
! 712: END IF
! 713: END IF
! 714: *
! 715: MINWRK = MAX( 2, MINWRK )
! 716: OPTWRK = MAX( 2, OPTWRK )
! 717: IF ( LCWORK .LT. MINWRK .AND. (.NOT.LQUERY) ) INFO = -19
! 718: *
! 719: END IF
! 720: *
! 721: IF (INFO .EQ. 0 .AND. LRWORK .LT. RMINWRK .AND. .NOT. LQUERY) THEN
! 722: INFO = -21
! 723: END IF
! 724: IF( INFO.NE.0 ) THEN
! 725: CALL XERBLA( 'ZGESVDQ', -INFO )
! 726: RETURN
! 727: ELSE IF ( LQUERY ) THEN
! 728: *
! 729: * Return optimal workspace
! 730: *
! 731: IWORK(1) = IMINWRK
! 732: CWORK(1) = OPTWRK
! 733: CWORK(2) = MINWRK
! 734: RWORK(1) = RMINWRK
! 735: RETURN
! 736: END IF
! 737: *
! 738: * Quick return if the matrix is void.
! 739: *
! 740: IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) ) THEN
! 741: * .. all output is void.
! 742: RETURN
! 743: END IF
! 744: *
! 745: BIG = DLAMCH('O')
! 746: ASCALED = .FALSE.
! 747: IF ( ROWPRM ) THEN
! 748: * .. reordering the rows in decreasing sequence in the
! 749: * ell-infinity norm - this enhances numerical robustness in
! 750: * the case of differently scaled rows.
! 751: DO 1904 p = 1, M
! 752: * RWORK(p) = ABS( A(p,IZAMAX(N,A(p,1),LDA)) )
! 753: * [[ZLANGE will return NaN if an entry of the p-th row is Nan]]
! 754: RWORK(p) = ZLANGE( 'M', 1, N, A(p,1), LDA, RDUMMY )
! 755: * .. check for NaN's and Inf's
! 756: IF ( ( RWORK(p) .NE. RWORK(p) ) .OR.
! 757: $ ( (RWORK(p)*ZERO) .NE. ZERO ) ) THEN
! 758: INFO = -8
! 759: CALL XERBLA( 'ZGESVDQ', -INFO )
! 760: RETURN
! 761: END IF
! 762: 1904 CONTINUE
! 763: DO 1952 p = 1, M - 1
! 764: q = IDAMAX( M-p+1, RWORK(p), 1 ) + p - 1
! 765: IWORK(N+p) = q
! 766: IF ( p .NE. q ) THEN
! 767: RTMP = RWORK(p)
! 768: RWORK(p) = RWORK(q)
! 769: RWORK(q) = RTMP
! 770: END IF
! 771: 1952 CONTINUE
! 772: *
! 773: IF ( RWORK(1) .EQ. ZERO ) THEN
! 774: * Quick return: A is the M x N zero matrix.
! 775: NUMRANK = 0
! 776: CALL DLASET( 'G', N, 1, ZERO, ZERO, S, N )
! 777: IF ( WNTUS ) CALL ZLASET('G', M, N, CZERO, CONE, U, LDU)
! 778: IF ( WNTUA ) CALL ZLASET('G', M, M, CZERO, CONE, U, LDU)
! 779: IF ( WNTVA ) CALL ZLASET('G', N, N, CZERO, CONE, V, LDV)
! 780: IF ( WNTUF ) THEN
! 781: CALL ZLASET( 'G', N, 1, CZERO, CZERO, CWORK, N )
! 782: CALL ZLASET( 'G', M, N, CZERO, CONE, U, LDU )
! 783: END IF
! 784: DO 5001 p = 1, N
! 785: IWORK(p) = p
! 786: 5001 CONTINUE
! 787: IF ( ROWPRM ) THEN
! 788: DO 5002 p = N + 1, N + M - 1
! 789: IWORK(p) = p - N
! 790: 5002 CONTINUE
! 791: END IF
! 792: IF ( CONDA ) RWORK(1) = -1
! 793: RWORK(2) = -1
! 794: RETURN
! 795: END IF
! 796: *
! 797: IF ( RWORK(1) .GT. BIG / SQRT(DBLE(M)) ) THEN
! 798: * .. to prevent overflow in the QR factorization, scale the
! 799: * matrix by 1/sqrt(M) if too large entry detected
! 800: CALL ZLASCL('G',0,0,SQRT(DBLE(M)),ONE, M,N, A,LDA, IERR)
! 801: ASCALED = .TRUE.
! 802: END IF
! 803: CALL ZLASWP( N, A, LDA, 1, M-1, IWORK(N+1), 1 )
! 804: END IF
! 805: *
! 806: * .. At this stage, preemptive scaling is done only to avoid column
! 807: * norms overflows during the QR factorization. The SVD procedure should
! 808: * have its own scaling to save the singular values from overflows and
! 809: * underflows. That depends on the SVD procedure.
! 810: *
! 811: IF ( .NOT.ROWPRM ) THEN
! 812: RTMP = ZLANGE( 'M', M, N, A, LDA, RWORK )
! 813: IF ( ( RTMP .NE. RTMP ) .OR.
! 814: $ ( (RTMP*ZERO) .NE. ZERO ) ) THEN
! 815: INFO = -8
! 816: CALL XERBLA( 'ZGESVDQ', -INFO )
! 817: RETURN
! 818: END IF
! 819: IF ( RTMP .GT. BIG / SQRT(DBLE(M)) ) THEN
! 820: * .. to prevent overflow in the QR factorization, scale the
! 821: * matrix by 1/sqrt(M) if too large entry detected
! 822: CALL ZLASCL('G',0,0, SQRT(DBLE(M)),ONE, M,N, A,LDA, IERR)
! 823: ASCALED = .TRUE.
! 824: END IF
! 825: END IF
! 826: *
! 827: * .. QR factorization with column pivoting
! 828: *
! 829: * A * P = Q * [ R ]
! 830: * [ 0 ]
! 831: *
! 832: DO 1963 p = 1, N
! 833: * .. all columns are free columns
! 834: IWORK(p) = 0
! 835: 1963 CONTINUE
! 836: CALL ZGEQP3( M, N, A, LDA, IWORK, CWORK, CWORK(N+1), LCWORK-N,
! 837: $ RWORK, IERR )
! 838: *
! 839: * If the user requested accuracy level allows truncation in the
! 840: * computed upper triangular factor, the matrix R is examined and,
! 841: * if possible, replaced with its leading upper trapezoidal part.
! 842: *
! 843: EPSLN = DLAMCH('E')
! 844: SFMIN = DLAMCH('S')
! 845: * SMALL = SFMIN / EPSLN
! 846: NR = N
! 847: *
! 848: IF ( ACCLA ) THEN
! 849: *
! 850: * Standard absolute error bound suffices. All sigma_i with
! 851: * sigma_i < N*EPS*||A||_F are flushed to zero. This is an
! 852: * aggressive enforcement of lower numerical rank by introducing a
! 853: * backward error of the order of N*EPS*||A||_F.
! 854: NR = 1
! 855: RTMP = SQRT(DBLE(N))*EPSLN
! 856: DO 3001 p = 2, N
! 857: IF ( ABS(A(p,p)) .LT. (RTMP*ABS(A(1,1))) ) GO TO 3002
! 858: NR = NR + 1
! 859: 3001 CONTINUE
! 860: 3002 CONTINUE
! 861: *
! 862: ELSEIF ( ACCLM ) THEN
! 863: * .. similarly as above, only slightly more gentle (less aggressive).
! 864: * Sudden drop on the diagonal of R is used as the criterion for being
! 865: * close-to-rank-deficient. The threshold is set to EPSLN=DLAMCH('E').
! 866: * [[This can be made more flexible by replacing this hard-coded value
! 867: * with a user specified threshold.]] Also, the values that underflow
! 868: * will be truncated.
! 869: NR = 1
! 870: DO 3401 p = 2, N
! 871: IF ( ( ABS(A(p,p)) .LT. (EPSLN*ABS(A(p-1,p-1))) ) .OR.
! 872: $ ( ABS(A(p,p)) .LT. SFMIN ) ) GO TO 3402
! 873: NR = NR + 1
! 874: 3401 CONTINUE
! 875: 3402 CONTINUE
! 876: *
! 877: ELSE
! 878: * .. RRQR not authorized to determine numerical rank except in the
! 879: * obvious case of zero pivots.
! 880: * .. inspect R for exact zeros on the diagonal;
! 881: * R(i,i)=0 => R(i:N,i:N)=0.
! 882: NR = 1
! 883: DO 3501 p = 2, N
! 884: IF ( ABS(A(p,p)) .EQ. ZERO ) GO TO 3502
! 885: NR = NR + 1
! 886: 3501 CONTINUE
! 887: 3502 CONTINUE
! 888: *
! 889: IF ( CONDA ) THEN
! 890: * Estimate the scaled condition number of A. Use the fact that it is
! 891: * the same as the scaled condition number of R.
! 892: * .. V is used as workspace
! 893: CALL ZLACPY( 'U', N, N, A, LDA, V, LDV )
! 894: * Only the leading NR x NR submatrix of the triangular factor
! 895: * is considered. Only if NR=N will this give a reliable error
! 896: * bound. However, even for NR < N, this can be used on an
! 897: * expert level and obtain useful information in the sense of
! 898: * perturbation theory.
! 899: DO 3053 p = 1, NR
! 900: RTMP = DZNRM2( p, V(1,p), 1 )
! 901: CALL ZDSCAL( p, ONE/RTMP, V(1,p), 1 )
! 902: 3053 CONTINUE
! 903: IF ( .NOT. ( LSVEC .OR. RSVEC ) ) THEN
! 904: CALL ZPOCON( 'U', NR, V, LDV, ONE, RTMP,
! 905: $ CWORK, RWORK, IERR )
! 906: ELSE
! 907: CALL ZPOCON( 'U', NR, V, LDV, ONE, RTMP,
! 908: $ CWORK(N+1), RWORK, IERR )
! 909: END IF
! 910: SCONDA = ONE / SQRT(RTMP)
! 911: * For NR=N, SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1),
! 912: * N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
! 913: * See the reference [1] for more details.
! 914: END IF
! 915: *
! 916: ENDIF
! 917: *
! 918: IF ( WNTUR ) THEN
! 919: N1 = NR
! 920: ELSE IF ( WNTUS .OR. WNTUF) THEN
! 921: N1 = N
! 922: ELSE IF ( WNTUA ) THEN
! 923: N1 = M
! 924: END IF
! 925: *
! 926: IF ( .NOT. ( RSVEC .OR. LSVEC ) ) THEN
! 927: *.......................................................................
! 928: * .. only the singular values are requested
! 929: *.......................................................................
! 930: IF ( RTRANS ) THEN
! 931: *
! 932: * .. compute the singular values of R**H = [A](1:NR,1:N)**H
! 933: * .. set the lower triangle of [A] to [A](1:NR,1:N)**H and
! 934: * the upper triangle of [A] to zero.
! 935: DO 1146 p = 1, MIN( N, NR )
! 936: A(p,p) = CONJG(A(p,p))
! 937: DO 1147 q = p + 1, N
! 938: A(q,p) = CONJG(A(p,q))
! 939: IF ( q .LE. NR ) A(p,q) = CZERO
! 940: 1147 CONTINUE
! 941: 1146 CONTINUE
! 942: *
! 943: CALL ZGESVD( 'N', 'N', N, NR, A, LDA, S, U, LDU,
! 944: $ V, LDV, CWORK, LCWORK, RWORK, INFO )
! 945: *
! 946: ELSE
! 947: *
! 948: * .. compute the singular values of R = [A](1:NR,1:N)
! 949: *
! 950: IF ( NR .GT. 1 )
! 951: $ CALL ZLASET( 'L', NR-1,NR-1, CZERO,CZERO, A(2,1), LDA )
! 952: CALL ZGESVD( 'N', 'N', NR, N, A, LDA, S, U, LDU,
! 953: $ V, LDV, CWORK, LCWORK, RWORK, INFO )
! 954: *
! 955: END IF
! 956: *
! 957: ELSE IF ( LSVEC .AND. ( .NOT. RSVEC) ) THEN
! 958: *.......................................................................
! 959: * .. the singular values and the left singular vectors requested
! 960: *.......................................................................""""""""
! 961: IF ( RTRANS ) THEN
! 962: * .. apply ZGESVD to R**H
! 963: * .. copy R**H into [U] and overwrite [U] with the right singular
! 964: * vectors of R
! 965: DO 1192 p = 1, NR
! 966: DO 1193 q = p, N
! 967: U(q,p) = CONJG(A(p,q))
! 968: 1193 CONTINUE
! 969: 1192 CONTINUE
! 970: IF ( NR .GT. 1 )
! 971: $ CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, U(1,2), LDU )
! 972: * .. the left singular vectors not computed, the NR right singular
! 973: * vectors overwrite [U](1:NR,1:NR) as conjugate transposed. These
! 974: * will be pre-multiplied by Q to build the left singular vectors of A.
! 975: CALL ZGESVD( 'N', 'O', N, NR, U, LDU, S, U, LDU,
! 976: $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO )
! 977: *
! 978: DO 1119 p = 1, NR
! 979: U(p,p) = CONJG(U(p,p))
! 980: DO 1120 q = p + 1, NR
! 981: CTMP = CONJG(U(q,p))
! 982: U(q,p) = CONJG(U(p,q))
! 983: U(p,q) = CTMP
! 984: 1120 CONTINUE
! 985: 1119 CONTINUE
! 986: *
! 987: ELSE
! 988: * .. apply ZGESVD to R
! 989: * .. copy R into [U] and overwrite [U] with the left singular vectors
! 990: CALL ZLACPY( 'U', NR, N, A, LDA, U, LDU )
! 991: IF ( NR .GT. 1 )
! 992: $ CALL ZLASET( 'L', NR-1, NR-1, CZERO, CZERO, U(2,1), LDU )
! 993: * .. the right singular vectors not computed, the NR left singular
! 994: * vectors overwrite [U](1:NR,1:NR)
! 995: CALL ZGESVD( 'O', 'N', NR, N, U, LDU, S, U, LDU,
! 996: $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO )
! 997: * .. now [U](1:NR,1:NR) contains the NR left singular vectors of
! 998: * R. These will be pre-multiplied by Q to build the left singular
! 999: * vectors of A.
! 1000: END IF
! 1001: *
! 1002: * .. assemble the left singular vector matrix U of dimensions
! 1003: * (M x NR) or (M x N) or (M x M).
! 1004: IF ( ( NR .LT. M ) .AND. ( .NOT.WNTUF ) ) THEN
! 1005: CALL ZLASET('A', M-NR, NR, CZERO, CZERO, U(NR+1,1), LDU)
! 1006: IF ( NR .LT. N1 ) THEN
! 1007: CALL ZLASET( 'A',NR,N1-NR,CZERO,CZERO,U(1,NR+1), LDU )
! 1008: CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,
! 1009: $ U(NR+1,NR+1), LDU )
! 1010: END IF
! 1011: END IF
! 1012: *
! 1013: * The Q matrix from the first QRF is built into the left singular
! 1014: * vectors matrix U.
! 1015: *
! 1016: IF ( .NOT.WNTUF )
! 1017: $ CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
! 1018: $ LDU, CWORK(N+1), LCWORK-N, IERR )
! 1019: IF ( ROWPRM .AND. .NOT.WNTUF )
! 1020: $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 )
! 1021: *
! 1022: ELSE IF ( RSVEC .AND. ( .NOT. LSVEC ) ) THEN
! 1023: *.......................................................................
! 1024: * .. the singular values and the right singular vectors requested
! 1025: *.......................................................................
! 1026: IF ( RTRANS ) THEN
! 1027: * .. apply ZGESVD to R**H
! 1028: * .. copy R**H into V and overwrite V with the left singular vectors
! 1029: DO 1165 p = 1, NR
! 1030: DO 1166 q = p, N
! 1031: V(q,p) = CONJG(A(p,q))
! 1032: 1166 CONTINUE
! 1033: 1165 CONTINUE
! 1034: IF ( NR .GT. 1 )
! 1035: $ CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV )
! 1036: * .. the left singular vectors of R**H overwrite V, the right singular
! 1037: * vectors not computed
! 1038: IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
! 1039: CALL ZGESVD( 'O', 'N', N, NR, V, LDV, S, U, LDU,
! 1040: $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO )
! 1041: *
! 1042: DO 1121 p = 1, NR
! 1043: V(p,p) = CONJG(V(p,p))
! 1044: DO 1122 q = p + 1, NR
! 1045: CTMP = CONJG(V(q,p))
! 1046: V(q,p) = CONJG(V(p,q))
! 1047: V(p,q) = CTMP
! 1048: 1122 CONTINUE
! 1049: 1121 CONTINUE
! 1050: *
! 1051: IF ( NR .LT. N ) THEN
! 1052: DO 1103 p = 1, NR
! 1053: DO 1104 q = NR + 1, N
! 1054: V(p,q) = CONJG(V(q,p))
! 1055: 1104 CONTINUE
! 1056: 1103 CONTINUE
! 1057: END IF
! 1058: CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK )
! 1059: ELSE
! 1060: * .. need all N right singular vectors and NR < N
! 1061: * [!] This is simple implementation that augments [V](1:N,1:NR)
! 1062: * by padding a zero block. In the case NR << N, a more efficient
! 1063: * way is to first use the QR factorization. For more details
! 1064: * how to implement this, see the " FULL SVD " branch.
! 1065: CALL ZLASET('G', N, N-NR, CZERO, CZERO, V(1,NR+1), LDV)
! 1066: CALL ZGESVD( 'O', 'N', N, N, V, LDV, S, U, LDU,
! 1067: $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO )
! 1068: *
! 1069: DO 1123 p = 1, N
! 1070: V(p,p) = CONJG(V(p,p))
! 1071: DO 1124 q = p + 1, N
! 1072: CTMP = CONJG(V(q,p))
! 1073: V(q,p) = CONJG(V(p,q))
! 1074: V(p,q) = CTMP
! 1075: 1124 CONTINUE
! 1076: 1123 CONTINUE
! 1077: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
! 1078: END IF
! 1079: *
! 1080: ELSE
! 1081: * .. aply ZGESVD to R
! 1082: * .. copy R into V and overwrite V with the right singular vectors
! 1083: CALL ZLACPY( 'U', NR, N, A, LDA, V, LDV )
! 1084: IF ( NR .GT. 1 )
! 1085: $ CALL ZLASET( 'L', NR-1, NR-1, CZERO, CZERO, V(2,1), LDV )
! 1086: * .. the right singular vectors overwrite V, the NR left singular
! 1087: * vectors stored in U(1:NR,1:NR)
! 1088: IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
! 1089: CALL ZGESVD( 'N', 'O', NR, N, V, LDV, S, U, LDU,
! 1090: $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO )
! 1091: CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK )
! 1092: * .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H
! 1093: ELSE
! 1094: * .. need all N right singular vectors and NR < N
! 1095: * [!] This is simple implementation that augments [V](1:NR,1:N)
! 1096: * by padding a zero block. In the case NR << N, a more efficient
! 1097: * way is to first use the LQ factorization. For more details
! 1098: * how to implement this, see the " FULL SVD " branch.
! 1099: CALL ZLASET('G', N-NR, N, CZERO,CZERO, V(NR+1,1), LDV)
! 1100: CALL ZGESVD( 'N', 'O', N, N, V, LDV, S, U, LDU,
! 1101: $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO )
! 1102: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
! 1103: END IF
! 1104: * .. now [V] contains the adjoint of the matrix of the right singular
! 1105: * vectors of A.
! 1106: END IF
! 1107: *
! 1108: ELSE
! 1109: *.......................................................................
! 1110: * .. FULL SVD requested
! 1111: *.......................................................................
! 1112: IF ( RTRANS ) THEN
! 1113: *
! 1114: * .. apply ZGESVD to R**H [[this option is left for R&D&T]]
! 1115: *
! 1116: IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
! 1117: * .. copy R**H into [V] and overwrite [V] with the left singular
! 1118: * vectors of R**H
! 1119: DO 1168 p = 1, NR
! 1120: DO 1169 q = p, N
! 1121: V(q,p) = CONJG(A(p,q))
! 1122: 1169 CONTINUE
! 1123: 1168 CONTINUE
! 1124: IF ( NR .GT. 1 )
! 1125: $ CALL ZLASET( 'U', NR-1,NR-1, CZERO,CZERO, V(1,2), LDV )
! 1126: *
! 1127: * .. the left singular vectors of R**H overwrite [V], the NR right
! 1128: * singular vectors of R**H stored in [U](1:NR,1:NR) as conjugate
! 1129: * transposed
! 1130: CALL ZGESVD( 'O', 'A', N, NR, V, LDV, S, V, LDV,
! 1131: $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO )
! 1132: * .. assemble V
! 1133: DO 1115 p = 1, NR
! 1134: V(p,p) = CONJG(V(p,p))
! 1135: DO 1116 q = p + 1, NR
! 1136: CTMP = CONJG(V(q,p))
! 1137: V(q,p) = CONJG(V(p,q))
! 1138: V(p,q) = CTMP
! 1139: 1116 CONTINUE
! 1140: 1115 CONTINUE
! 1141: IF ( NR .LT. N ) THEN
! 1142: DO 1101 p = 1, NR
! 1143: DO 1102 q = NR+1, N
! 1144: V(p,q) = CONJG(V(q,p))
! 1145: 1102 CONTINUE
! 1146: 1101 CONTINUE
! 1147: END IF
! 1148: CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK )
! 1149: *
! 1150: DO 1117 p = 1, NR
! 1151: U(p,p) = CONJG(U(p,p))
! 1152: DO 1118 q = p + 1, NR
! 1153: CTMP = CONJG(U(q,p))
! 1154: U(q,p) = CONJG(U(p,q))
! 1155: U(p,q) = CTMP
! 1156: 1118 CONTINUE
! 1157: 1117 CONTINUE
! 1158: *
! 1159: IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
! 1160: CALL ZLASET('A', M-NR,NR, CZERO,CZERO, U(NR+1,1), LDU)
! 1161: IF ( NR .LT. N1 ) THEN
! 1162: CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
! 1163: CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,
! 1164: $ U(NR+1,NR+1), LDU )
! 1165: END IF
! 1166: END IF
! 1167: *
! 1168: ELSE
! 1169: * .. need all N right singular vectors and NR < N
! 1170: * .. copy R**H into [V] and overwrite [V] with the left singular
! 1171: * vectors of R**H
! 1172: * [[The optimal ratio N/NR for using QRF instead of padding
! 1173: * with zeros. Here hard coded to 2; it must be at least
! 1174: * two due to work space constraints.]]
! 1175: * OPTRATIO = ILAENV(6, 'ZGESVD', 'S' // 'O', NR,N,0,0)
! 1176: * OPTRATIO = MAX( OPTRATIO, 2 )
! 1177: OPTRATIO = 2
! 1178: IF ( OPTRATIO*NR .GT. N ) THEN
! 1179: DO 1198 p = 1, NR
! 1180: DO 1199 q = p, N
! 1181: V(q,p) = CONJG(A(p,q))
! 1182: 1199 CONTINUE
! 1183: 1198 CONTINUE
! 1184: IF ( NR .GT. 1 )
! 1185: $ CALL ZLASET('U',NR-1,NR-1, CZERO,CZERO, V(1,2),LDV)
! 1186: *
! 1187: CALL ZLASET('A',N,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
! 1188: CALL ZGESVD( 'O', 'A', N, N, V, LDV, S, V, LDV,
! 1189: $ U, LDU, CWORK(N+1), LCWORK-N, RWORK, INFO )
! 1190: *
! 1191: DO 1113 p = 1, N
! 1192: V(p,p) = CONJG(V(p,p))
! 1193: DO 1114 q = p + 1, N
! 1194: CTMP = CONJG(V(q,p))
! 1195: V(q,p) = CONJG(V(p,q))
! 1196: V(p,q) = CTMP
! 1197: 1114 CONTINUE
! 1198: 1113 CONTINUE
! 1199: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
! 1200: * .. assemble the left singular vector matrix U of dimensions
! 1201: * (M x N1), i.e. (M x N) or (M x M).
! 1202: *
! 1203: DO 1111 p = 1, N
! 1204: U(p,p) = CONJG(U(p,p))
! 1205: DO 1112 q = p + 1, N
! 1206: CTMP = CONJG(U(q,p))
! 1207: U(q,p) = CONJG(U(p,q))
! 1208: U(p,q) = CTMP
! 1209: 1112 CONTINUE
! 1210: 1111 CONTINUE
! 1211: *
! 1212: IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN
! 1213: CALL ZLASET('A',M-N,N,CZERO,CZERO,U(N+1,1),LDU)
! 1214: IF ( N .LT. N1 ) THEN
! 1215: CALL ZLASET('A',N,N1-N,CZERO,CZERO,U(1,N+1),LDU)
! 1216: CALL ZLASET('A',M-N,N1-N,CZERO,CONE,
! 1217: $ U(N+1,N+1), LDU )
! 1218: END IF
! 1219: END IF
! 1220: ELSE
! 1221: * .. copy R**H into [U] and overwrite [U] with the right
! 1222: * singular vectors of R
! 1223: DO 1196 p = 1, NR
! 1224: DO 1197 q = p, N
! 1225: U(q,NR+p) = CONJG(A(p,q))
! 1226: 1197 CONTINUE
! 1227: 1196 CONTINUE
! 1228: IF ( NR .GT. 1 )
! 1229: $ CALL ZLASET('U',NR-1,NR-1,CZERO,CZERO,U(1,NR+2),LDU)
! 1230: CALL ZGEQRF( N, NR, U(1,NR+1), LDU, CWORK(N+1),
! 1231: $ CWORK(N+NR+1), LCWORK-N-NR, IERR )
! 1232: DO 1143 p = 1, NR
! 1233: DO 1144 q = 1, N
! 1234: V(q,p) = CONJG(U(p,NR+q))
! 1235: 1144 CONTINUE
! 1236: 1143 CONTINUE
! 1237: CALL ZLASET('U',NR-1,NR-1,CZERO,CZERO,V(1,2),LDV)
! 1238: CALL ZGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU,
! 1239: $ V,LDV, CWORK(N+NR+1),LCWORK-N-NR,RWORK, INFO )
! 1240: CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV)
! 1241: CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
! 1242: CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
! 1243: CALL ZUNMQR('R','C', N, N, NR, U(1,NR+1), LDU,
! 1244: $ CWORK(N+1),V,LDV,CWORK(N+NR+1),LCWORK-N-NR,IERR)
! 1245: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
! 1246: * .. assemble the left singular vector matrix U of dimensions
! 1247: * (M x NR) or (M x N) or (M x M).
! 1248: IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
! 1249: CALL ZLASET('A',M-NR,NR,CZERO,CZERO,U(NR+1,1),LDU)
! 1250: IF ( NR .LT. N1 ) THEN
! 1251: CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
! 1252: CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,
! 1253: $ U(NR+1,NR+1),LDU)
! 1254: END IF
! 1255: END IF
! 1256: END IF
! 1257: END IF
! 1258: *
! 1259: ELSE
! 1260: *
! 1261: * .. apply ZGESVD to R [[this is the recommended option]]
! 1262: *
! 1263: IF ( WNTVR .OR. ( NR .EQ. N ) ) THEN
! 1264: * .. copy R into [V] and overwrite V with the right singular vectors
! 1265: CALL ZLACPY( 'U', NR, N, A, LDA, V, LDV )
! 1266: IF ( NR .GT. 1 )
! 1267: $ CALL ZLASET( 'L', NR-1,NR-1, CZERO,CZERO, V(2,1), LDV )
! 1268: * .. the right singular vectors of R overwrite [V], the NR left
! 1269: * singular vectors of R stored in [U](1:NR,1:NR)
! 1270: CALL ZGESVD( 'S', 'O', NR, N, V, LDV, S, U, LDU,
! 1271: $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO )
! 1272: CALL ZLAPMT( .FALSE., NR, N, V, LDV, IWORK )
! 1273: * .. now [V](1:NR,1:N) contains V(1:N,1:NR)**H
! 1274: * .. assemble the left singular vector matrix U of dimensions
! 1275: * (M x NR) or (M x N) or (M x M).
! 1276: IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
! 1277: CALL ZLASET('A', M-NR,NR, CZERO,CZERO, U(NR+1,1), LDU)
! 1278: IF ( NR .LT. N1 ) THEN
! 1279: CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
! 1280: CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,
! 1281: $ U(NR+1,NR+1), LDU )
! 1282: END IF
! 1283: END IF
! 1284: *
! 1285: ELSE
! 1286: * .. need all N right singular vectors and NR < N
! 1287: * .. the requested number of the left singular vectors
! 1288: * is then N1 (N or M)
! 1289: * [[The optimal ratio N/NR for using LQ instead of padding
! 1290: * with zeros. Here hard coded to 2; it must be at least
! 1291: * two due to work space constraints.]]
! 1292: * OPTRATIO = ILAENV(6, 'ZGESVD', 'S' // 'O', NR,N,0,0)
! 1293: * OPTRATIO = MAX( OPTRATIO, 2 )
! 1294: OPTRATIO = 2
! 1295: IF ( OPTRATIO * NR .GT. N ) THEN
! 1296: CALL ZLACPY( 'U', NR, N, A, LDA, V, LDV )
! 1297: IF ( NR .GT. 1 )
! 1298: $ CALL ZLASET('L', NR-1,NR-1, CZERO,CZERO, V(2,1),LDV)
! 1299: * .. the right singular vectors of R overwrite [V], the NR left
! 1300: * singular vectors of R stored in [U](1:NR,1:NR)
! 1301: CALL ZLASET('A', N-NR,N, CZERO,CZERO, V(NR+1,1),LDV)
! 1302: CALL ZGESVD( 'S', 'O', N, N, V, LDV, S, U, LDU,
! 1303: $ V, LDV, CWORK(N+1), LCWORK-N, RWORK, INFO )
! 1304: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
! 1305: * .. now [V] contains the adjoint of the matrix of the right
! 1306: * singular vectors of A. The leading N left singular vectors
! 1307: * are in [U](1:N,1:N)
! 1308: * .. assemble the left singular vector matrix U of dimensions
! 1309: * (M x N1), i.e. (M x N) or (M x M).
! 1310: IF ( ( N .LT. M ) .AND. .NOT.(WNTUF)) THEN
! 1311: CALL ZLASET('A',M-N,N,CZERO,CZERO,U(N+1,1),LDU)
! 1312: IF ( N .LT. N1 ) THEN
! 1313: CALL ZLASET('A',N,N1-N,CZERO,CZERO,U(1,N+1),LDU)
! 1314: CALL ZLASET( 'A',M-N,N1-N,CZERO,CONE,
! 1315: $ U(N+1,N+1), LDU )
! 1316: END IF
! 1317: END IF
! 1318: ELSE
! 1319: CALL ZLACPY( 'U', NR, N, A, LDA, U(NR+1,1), LDU )
! 1320: IF ( NR .GT. 1 )
! 1321: $ CALL ZLASET('L',NR-1,NR-1,CZERO,CZERO,U(NR+2,1),LDU)
! 1322: CALL ZGELQF( NR, N, U(NR+1,1), LDU, CWORK(N+1),
! 1323: $ CWORK(N+NR+1), LCWORK-N-NR, IERR )
! 1324: CALL ZLACPY('L',NR,NR,U(NR+1,1),LDU,V,LDV)
! 1325: IF ( NR .GT. 1 )
! 1326: $ CALL ZLASET('U',NR-1,NR-1,CZERO,CZERO,V(1,2),LDV)
! 1327: CALL ZGESVD( 'S', 'O', NR, NR, V, LDV, S, U, LDU,
! 1328: $ V, LDV, CWORK(N+NR+1), LCWORK-N-NR, RWORK, INFO )
! 1329: CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV)
! 1330: CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
! 1331: CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
! 1332: CALL ZUNMLQ('R','N',N,N,NR,U(NR+1,1),LDU,CWORK(N+1),
! 1333: $ V, LDV, CWORK(N+NR+1),LCWORK-N-NR,IERR)
! 1334: CALL ZLAPMT( .FALSE., N, N, V, LDV, IWORK )
! 1335: * .. assemble the left singular vector matrix U of dimensions
! 1336: * (M x NR) or (M x N) or (M x M).
! 1337: IF ( ( NR .LT. M ) .AND. .NOT.(WNTUF)) THEN
! 1338: CALL ZLASET('A',M-NR,NR,CZERO,CZERO,U(NR+1,1),LDU)
! 1339: IF ( NR .LT. N1 ) THEN
! 1340: CALL ZLASET('A',NR,N1-NR,CZERO,CZERO,U(1,NR+1),LDU)
! 1341: CALL ZLASET( 'A',M-NR,N1-NR,CZERO,CONE,
! 1342: $ U(NR+1,NR+1), LDU )
! 1343: END IF
! 1344: END IF
! 1345: END IF
! 1346: END IF
! 1347: * .. end of the "R**H or R" branch
! 1348: END IF
! 1349: *
! 1350: * The Q matrix from the first QRF is built into the left singular
! 1351: * vectors matrix U.
! 1352: *
! 1353: IF ( .NOT. WNTUF )
! 1354: $ CALL ZUNMQR( 'L', 'N', M, N1, N, A, LDA, CWORK, U,
! 1355: $ LDU, CWORK(N+1), LCWORK-N, IERR )
! 1356: IF ( ROWPRM .AND. .NOT.WNTUF )
! 1357: $ CALL ZLASWP( N1, U, LDU, 1, M-1, IWORK(N+1), -1 )
! 1358: *
! 1359: * ... end of the "full SVD" branch
! 1360: END IF
! 1361: *
! 1362: * Check whether some singular values are returned as zeros, e.g.
! 1363: * due to underflow, and update the numerical rank.
! 1364: p = NR
! 1365: DO 4001 q = p, 1, -1
! 1366: IF ( S(q) .GT. ZERO ) GO TO 4002
! 1367: NR = NR - 1
! 1368: 4001 CONTINUE
! 1369: 4002 CONTINUE
! 1370: *
! 1371: * .. if numerical rank deficiency is detected, the truncated
! 1372: * singular values are set to zero.
! 1373: IF ( NR .LT. N ) CALL DLASET( 'G', N-NR,1, ZERO,ZERO, S(NR+1), N )
! 1374: * .. undo scaling; this may cause overflow in the largest singular
! 1375: * values.
! 1376: IF ( ASCALED )
! 1377: $ CALL DLASCL( 'G',0,0, ONE,SQRT(DBLE(M)), NR,1, S, N, IERR )
! 1378: IF ( CONDA ) RWORK(1) = SCONDA
! 1379: RWORK(2) = p - NR
! 1380: * .. p-NR is the number of singular values that are computed as
! 1381: * exact zeros in ZGESVD() applied to the (possibly truncated)
! 1382: * full row rank triangular (trapezoidal) factor of A.
! 1383: NUMRANK = NR
! 1384: *
! 1385: RETURN
! 1386: *
! 1387: * End of ZGESVDQ
! 1388: *
! 1389: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgesvdq.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack