Annotation of rpl/lapack/lapack/zgesvd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
! 2: $ WORK, LWORK, RWORK, INFO )
! 3: *
! 4: * -- LAPACK driver routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER JOBU, JOBVT
! 11: INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
! 12: * ..
! 13: * .. Array Arguments ..
! 14: DOUBLE PRECISION RWORK( * ), S( * )
! 15: COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
! 16: $ WORK( * )
! 17: * ..
! 18: *
! 19: * Purpose
! 20: * =======
! 21: *
! 22: * ZGESVD computes the singular value decomposition (SVD) of a complex
! 23: * M-by-N matrix A, optionally computing the left and/or right singular
! 24: * vectors. The SVD is written
! 25: *
! 26: * A = U * SIGMA * conjugate-transpose(V)
! 27: *
! 28: * where SIGMA is an M-by-N matrix which is zero except for its
! 29: * min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
! 30: * V is an N-by-N unitary matrix. The diagonal elements of SIGMA
! 31: * are the singular values of A; they are real and non-negative, and
! 32: * are returned in descending order. The first min(m,n) columns of
! 33: * U and V are the left and right singular vectors of A.
! 34: *
! 35: * Note that the routine returns V**H, not V.
! 36: *
! 37: * Arguments
! 38: * =========
! 39: *
! 40: * JOBU (input) CHARACTER*1
! 41: * Specifies options for computing all or part of the matrix U:
! 42: * = 'A': all M columns of U are returned in array U:
! 43: * = 'S': the first min(m,n) columns of U (the left singular
! 44: * vectors) are returned in the array U;
! 45: * = 'O': the first min(m,n) columns of U (the left singular
! 46: * vectors) are overwritten on the array A;
! 47: * = 'N': no columns of U (no left singular vectors) are
! 48: * computed.
! 49: *
! 50: * JOBVT (input) CHARACTER*1
! 51: * Specifies options for computing all or part of the matrix
! 52: * V**H:
! 53: * = 'A': all N rows of V**H are returned in the array VT;
! 54: * = 'S': the first min(m,n) rows of V**H (the right singular
! 55: * vectors) are returned in the array VT;
! 56: * = 'O': the first min(m,n) rows of V**H (the right singular
! 57: * vectors) are overwritten on the array A;
! 58: * = 'N': no rows of V**H (no right singular vectors) are
! 59: * computed.
! 60: *
! 61: * JOBVT and JOBU cannot both be 'O'.
! 62: *
! 63: * M (input) INTEGER
! 64: * The number of rows of the input matrix A. M >= 0.
! 65: *
! 66: * N (input) INTEGER
! 67: * The number of columns of the input matrix A. N >= 0.
! 68: *
! 69: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
! 70: * On entry, the M-by-N matrix A.
! 71: * On exit,
! 72: * if JOBU = 'O', A is overwritten with the first min(m,n)
! 73: * columns of U (the left singular vectors,
! 74: * stored columnwise);
! 75: * if JOBVT = 'O', A is overwritten with the first min(m,n)
! 76: * rows of V**H (the right singular vectors,
! 77: * stored rowwise);
! 78: * if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
! 79: * are destroyed.
! 80: *
! 81: * LDA (input) INTEGER
! 82: * The leading dimension of the array A. LDA >= max(1,M).
! 83: *
! 84: * S (output) DOUBLE PRECISION array, dimension (min(M,N))
! 85: * The singular values of A, sorted so that S(i) >= S(i+1).
! 86: *
! 87: * U (output) COMPLEX*16 array, dimension (LDU,UCOL)
! 88: * (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
! 89: * If JOBU = 'A', U contains the M-by-M unitary matrix U;
! 90: * if JOBU = 'S', U contains the first min(m,n) columns of U
! 91: * (the left singular vectors, stored columnwise);
! 92: * if JOBU = 'N' or 'O', U is not referenced.
! 93: *
! 94: * LDU (input) INTEGER
! 95: * The leading dimension of the array U. LDU >= 1; if
! 96: * JOBU = 'S' or 'A', LDU >= M.
! 97: *
! 98: * VT (output) COMPLEX*16 array, dimension (LDVT,N)
! 99: * If JOBVT = 'A', VT contains the N-by-N unitary matrix
! 100: * V**H;
! 101: * if JOBVT = 'S', VT contains the first min(m,n) rows of
! 102: * V**H (the right singular vectors, stored rowwise);
! 103: * if JOBVT = 'N' or 'O', VT is not referenced.
! 104: *
! 105: * LDVT (input) INTEGER
! 106: * The leading dimension of the array VT. LDVT >= 1; if
! 107: * JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
! 108: *
! 109: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
! 110: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 111: *
! 112: * LWORK (input) INTEGER
! 113: * The dimension of the array WORK.
! 114: * LWORK >= MAX(1,2*MIN(M,N)+MAX(M,N)).
! 115: * For good performance, LWORK should generally be larger.
! 116: *
! 117: * If LWORK = -1, then a workspace query is assumed; the routine
! 118: * only calculates the optimal size of the WORK array, returns
! 119: * this value as the first entry of the WORK array, and no error
! 120: * message related to LWORK is issued by XERBLA.
! 121: *
! 122: * RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,N))
! 123: * On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the
! 124: * unconverged superdiagonal elements of an upper bidiagonal
! 125: * matrix B whose diagonal is in S (not necessarily sorted).
! 126: * B satisfies A = U * B * VT, so it has the same singular
! 127: * values as A, and singular vectors related by U and VT.
! 128: *
! 129: * INFO (output) INTEGER
! 130: * = 0: successful exit.
! 131: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 132: * > 0: if ZBDSQR did not converge, INFO specifies how many
! 133: * superdiagonals of an intermediate bidiagonal form B
! 134: * did not converge to zero. See the description of RWORK
! 135: * above for details.
! 136: *
! 137: * =====================================================================
! 138: *
! 139: * .. Parameters ..
! 140: COMPLEX*16 CZERO, CONE
! 141: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
! 142: $ CONE = ( 1.0D0, 0.0D0 ) )
! 143: DOUBLE PRECISION ZERO, ONE
! 144: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 145: * ..
! 146: * .. Local Scalars ..
! 147: LOGICAL LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS,
! 148: $ WNTVA, WNTVAS, WNTVN, WNTVO, WNTVS
! 149: INTEGER BLK, CHUNK, I, IE, IERR, IR, IRWORK, ISCL,
! 150: $ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,
! 151: $ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU,
! 152: $ NRVT, WRKBL
! 153: DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
! 154: * ..
! 155: * .. Local Arrays ..
! 156: DOUBLE PRECISION DUM( 1 )
! 157: COMPLEX*16 CDUM( 1 )
! 158: * ..
! 159: * .. External Subroutines ..
! 160: EXTERNAL DLASCL, XERBLA, ZBDSQR, ZGEBRD, ZGELQF, ZGEMM,
! 161: $ ZGEQRF, ZLACPY, ZLASCL, ZLASET, ZUNGBR, ZUNGLQ,
! 162: $ ZUNGQR, ZUNMBR
! 163: * ..
! 164: * .. External Functions ..
! 165: LOGICAL LSAME
! 166: INTEGER ILAENV
! 167: DOUBLE PRECISION DLAMCH, ZLANGE
! 168: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
! 169: * ..
! 170: * .. Intrinsic Functions ..
! 171: INTRINSIC MAX, MIN, SQRT
! 172: * ..
! 173: * .. Executable Statements ..
! 174: *
! 175: * Test the input arguments
! 176: *
! 177: INFO = 0
! 178: MINMN = MIN( M, N )
! 179: WNTUA = LSAME( JOBU, 'A' )
! 180: WNTUS = LSAME( JOBU, 'S' )
! 181: WNTUAS = WNTUA .OR. WNTUS
! 182: WNTUO = LSAME( JOBU, 'O' )
! 183: WNTUN = LSAME( JOBU, 'N' )
! 184: WNTVA = LSAME( JOBVT, 'A' )
! 185: WNTVS = LSAME( JOBVT, 'S' )
! 186: WNTVAS = WNTVA .OR. WNTVS
! 187: WNTVO = LSAME( JOBVT, 'O' )
! 188: WNTVN = LSAME( JOBVT, 'N' )
! 189: LQUERY = ( LWORK.EQ.-1 )
! 190: *
! 191: IF( .NOT.( WNTUA .OR. WNTUS .OR. WNTUO .OR. WNTUN ) ) THEN
! 192: INFO = -1
! 193: ELSE IF( .NOT.( WNTVA .OR. WNTVS .OR. WNTVO .OR. WNTVN ) .OR.
! 194: $ ( WNTVO .AND. WNTUO ) ) THEN
! 195: INFO = -2
! 196: ELSE IF( M.LT.0 ) THEN
! 197: INFO = -3
! 198: ELSE IF( N.LT.0 ) THEN
! 199: INFO = -4
! 200: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 201: INFO = -6
! 202: ELSE IF( LDU.LT.1 .OR. ( WNTUAS .AND. LDU.LT.M ) ) THEN
! 203: INFO = -9
! 204: ELSE IF( LDVT.LT.1 .OR. ( WNTVA .AND. LDVT.LT.N ) .OR.
! 205: $ ( WNTVS .AND. LDVT.LT.MINMN ) ) THEN
! 206: INFO = -11
! 207: END IF
! 208: *
! 209: * Compute workspace
! 210: * (Note: Comments in the code beginning "Workspace:" describe the
! 211: * minimal amount of workspace needed at that point in the code,
! 212: * as well as the preferred amount for good performance.
! 213: * CWorkspace refers to complex workspace, and RWorkspace to
! 214: * real workspace. NB refers to the optimal block size for the
! 215: * immediately following subroutine, as returned by ILAENV.)
! 216: *
! 217: IF( INFO.EQ.0 ) THEN
! 218: MINWRK = 1
! 219: MAXWRK = 1
! 220: IF( M.GE.N .AND. MINMN.GT.0 ) THEN
! 221: *
! 222: * Space needed for ZBDSQR is BDSPAC = 5*N
! 223: *
! 224: MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )
! 225: IF( M.GE.MNTHR ) THEN
! 226: IF( WNTUN ) THEN
! 227: *
! 228: * Path 1 (M much larger than N, JOBU='N')
! 229: *
! 230: MAXWRK = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1,
! 231: $ -1 )
! 232: MAXWRK = MAX( MAXWRK, 2*N+2*N*
! 233: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 234: IF( WNTVO .OR. WNTVAS )
! 235: $ MAXWRK = MAX( MAXWRK, 2*N+( N-1 )*
! 236: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
! 237: MINWRK = 3*N
! 238: ELSE IF( WNTUO .AND. WNTVN ) THEN
! 239: *
! 240: * Path 2 (M much larger than N, JOBU='O', JOBVT='N')
! 241: *
! 242: WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 243: WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
! 244: $ N, N, -1 ) )
! 245: WRKBL = MAX( WRKBL, 2*N+2*N*
! 246: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 247: WRKBL = MAX( WRKBL, 2*N+N*
! 248: $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
! 249: MAXWRK = MAX( N*N+WRKBL, N*N+M*N )
! 250: MINWRK = 2*N + M
! 251: ELSE IF( WNTUO .AND. WNTVAS ) THEN
! 252: *
! 253: * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
! 254: * 'A')
! 255: *
! 256: WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 257: WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
! 258: $ N, N, -1 ) )
! 259: WRKBL = MAX( WRKBL, 2*N+2*N*
! 260: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 261: WRKBL = MAX( WRKBL, 2*N+N*
! 262: $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
! 263: WRKBL = MAX( WRKBL, 2*N+( N-1 )*
! 264: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
! 265: MAXWRK = MAX( N*N+WRKBL, N*N+M*N )
! 266: MINWRK = 2*N + M
! 267: ELSE IF( WNTUS .AND. WNTVN ) THEN
! 268: *
! 269: * Path 4 (M much larger than N, JOBU='S', JOBVT='N')
! 270: *
! 271: WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 272: WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
! 273: $ N, N, -1 ) )
! 274: WRKBL = MAX( WRKBL, 2*N+2*N*
! 275: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 276: WRKBL = MAX( WRKBL, 2*N+N*
! 277: $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
! 278: MAXWRK = N*N + WRKBL
! 279: MINWRK = 2*N + M
! 280: ELSE IF( WNTUS .AND. WNTVO ) THEN
! 281: *
! 282: * Path 5 (M much larger than N, JOBU='S', JOBVT='O')
! 283: *
! 284: WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 285: WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
! 286: $ N, N, -1 ) )
! 287: WRKBL = MAX( WRKBL, 2*N+2*N*
! 288: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 289: WRKBL = MAX( WRKBL, 2*N+N*
! 290: $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
! 291: WRKBL = MAX( WRKBL, 2*N+( N-1 )*
! 292: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
! 293: MAXWRK = 2*N*N + WRKBL
! 294: MINWRK = 2*N + M
! 295: ELSE IF( WNTUS .AND. WNTVAS ) THEN
! 296: *
! 297: * Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
! 298: * 'A')
! 299: *
! 300: WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 301: WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
! 302: $ N, N, -1 ) )
! 303: WRKBL = MAX( WRKBL, 2*N+2*N*
! 304: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 305: WRKBL = MAX( WRKBL, 2*N+N*
! 306: $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
! 307: WRKBL = MAX( WRKBL, 2*N+( N-1 )*
! 308: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
! 309: MAXWRK = N*N + WRKBL
! 310: MINWRK = 2*N + M
! 311: ELSE IF( WNTUA .AND. WNTVN ) THEN
! 312: *
! 313: * Path 7 (M much larger than N, JOBU='A', JOBVT='N')
! 314: *
! 315: WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 316: WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M,
! 317: $ M, N, -1 ) )
! 318: WRKBL = MAX( WRKBL, 2*N+2*N*
! 319: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 320: WRKBL = MAX( WRKBL, 2*N+N*
! 321: $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
! 322: MAXWRK = N*N + WRKBL
! 323: MINWRK = 2*N + M
! 324: ELSE IF( WNTUA .AND. WNTVO ) THEN
! 325: *
! 326: * Path 8 (M much larger than N, JOBU='A', JOBVT='O')
! 327: *
! 328: WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 329: WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M,
! 330: $ M, N, -1 ) )
! 331: WRKBL = MAX( WRKBL, 2*N+2*N*
! 332: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 333: WRKBL = MAX( WRKBL, 2*N+N*
! 334: $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
! 335: WRKBL = MAX( WRKBL, 2*N+( N-1 )*
! 336: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
! 337: MAXWRK = 2*N*N + WRKBL
! 338: MINWRK = 2*N + M
! 339: ELSE IF( WNTUA .AND. WNTVAS ) THEN
! 340: *
! 341: * Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
! 342: * 'A')
! 343: *
! 344: WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 345: WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M,
! 346: $ M, N, -1 ) )
! 347: WRKBL = MAX( WRKBL, 2*N+2*N*
! 348: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 349: WRKBL = MAX( WRKBL, 2*N+N*
! 350: $ ILAENV( 1, 'ZUNGBR', 'Q', N, N, N, -1 ) )
! 351: WRKBL = MAX( WRKBL, 2*N+( N-1 )*
! 352: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
! 353: MAXWRK = N*N + WRKBL
! 354: MINWRK = 2*N + M
! 355: END IF
! 356: ELSE
! 357: *
! 358: * Path 10 (M at least N, but not much larger)
! 359: *
! 360: MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,
! 361: $ -1, -1 )
! 362: IF( WNTUS .OR. WNTUO )
! 363: $ MAXWRK = MAX( MAXWRK, 2*N+N*
! 364: $ ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) )
! 365: IF( WNTUA )
! 366: $ MAXWRK = MAX( MAXWRK, 2*N+M*
! 367: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )
! 368: IF( .NOT.WNTVN )
! 369: $ MAXWRK = MAX( MAXWRK, 2*N+( N-1 )*
! 370: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
! 371: MINWRK = 2*N + M
! 372: END IF
! 373: ELSE IF( MINMN.GT.0 ) THEN
! 374: *
! 375: * Space needed for ZBDSQR is BDSPAC = 5*M
! 376: *
! 377: MNTHR = ILAENV( 6, 'ZGESVD', JOBU // JOBVT, M, N, 0, 0 )
! 378: IF( N.GE.MNTHR ) THEN
! 379: IF( WNTVN ) THEN
! 380: *
! 381: * Path 1t(N much larger than M, JOBVT='N')
! 382: *
! 383: MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
! 384: $ -1 )
! 385: MAXWRK = MAX( MAXWRK, 2*M+2*M*
! 386: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 387: IF( WNTUO .OR. WNTUAS )
! 388: $ MAXWRK = MAX( MAXWRK, 2*M+M*
! 389: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
! 390: MINWRK = 3*M
! 391: ELSE IF( WNTVO .AND. WNTUN ) THEN
! 392: *
! 393: * Path 2t(N much larger than M, JOBU='N', JOBVT='O')
! 394: *
! 395: WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 396: WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
! 397: $ N, M, -1 ) )
! 398: WRKBL = MAX( WRKBL, 2*M+2*M*
! 399: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 400: WRKBL = MAX( WRKBL, 2*M+( M-1 )*
! 401: $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
! 402: MAXWRK = MAX( M*M+WRKBL, M*M+M*N )
! 403: MINWRK = 2*M + N
! 404: ELSE IF( WNTVO .AND. WNTUAS ) THEN
! 405: *
! 406: * Path 3t(N much larger than M, JOBU='S' or 'A',
! 407: * JOBVT='O')
! 408: *
! 409: WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 410: WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
! 411: $ N, M, -1 ) )
! 412: WRKBL = MAX( WRKBL, 2*M+2*M*
! 413: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 414: WRKBL = MAX( WRKBL, 2*M+( M-1 )*
! 415: $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
! 416: WRKBL = MAX( WRKBL, 2*M+M*
! 417: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
! 418: MAXWRK = MAX( M*M+WRKBL, M*M+M*N )
! 419: MINWRK = 2*M + N
! 420: ELSE IF( WNTVS .AND. WNTUN ) THEN
! 421: *
! 422: * Path 4t(N much larger than M, JOBU='N', JOBVT='S')
! 423: *
! 424: WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 425: WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
! 426: $ N, M, -1 ) )
! 427: WRKBL = MAX( WRKBL, 2*M+2*M*
! 428: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 429: WRKBL = MAX( WRKBL, 2*M+( M-1 )*
! 430: $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
! 431: MAXWRK = M*M + WRKBL
! 432: MINWRK = 2*M + N
! 433: ELSE IF( WNTVS .AND. WNTUO ) THEN
! 434: *
! 435: * Path 5t(N much larger than M, JOBU='O', JOBVT='S')
! 436: *
! 437: WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 438: WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
! 439: $ N, M, -1 ) )
! 440: WRKBL = MAX( WRKBL, 2*M+2*M*
! 441: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 442: WRKBL = MAX( WRKBL, 2*M+( M-1 )*
! 443: $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
! 444: WRKBL = MAX( WRKBL, 2*M+M*
! 445: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
! 446: MAXWRK = 2*M*M + WRKBL
! 447: MINWRK = 2*M + N
! 448: ELSE IF( WNTVS .AND. WNTUAS ) THEN
! 449: *
! 450: * Path 6t(N much larger than M, JOBU='S' or 'A',
! 451: * JOBVT='S')
! 452: *
! 453: WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 454: WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
! 455: $ N, M, -1 ) )
! 456: WRKBL = MAX( WRKBL, 2*M+2*M*
! 457: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 458: WRKBL = MAX( WRKBL, 2*M+( M-1 )*
! 459: $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
! 460: WRKBL = MAX( WRKBL, 2*M+M*
! 461: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
! 462: MAXWRK = M*M + WRKBL
! 463: MINWRK = 2*M + N
! 464: ELSE IF( WNTVA .AND. WNTUN ) THEN
! 465: *
! 466: * Path 7t(N much larger than M, JOBU='N', JOBVT='A')
! 467: *
! 468: WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 469: WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N,
! 470: $ N, M, -1 ) )
! 471: WRKBL = MAX( WRKBL, 2*M+2*M*
! 472: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 473: WRKBL = MAX( WRKBL, 2*M+( M-1 )*
! 474: $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
! 475: MAXWRK = M*M + WRKBL
! 476: MINWRK = 2*M + N
! 477: ELSE IF( WNTVA .AND. WNTUO ) THEN
! 478: *
! 479: * Path 8t(N much larger than M, JOBU='O', JOBVT='A')
! 480: *
! 481: WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 482: WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N,
! 483: $ N, M, -1 ) )
! 484: WRKBL = MAX( WRKBL, 2*M+2*M*
! 485: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 486: WRKBL = MAX( WRKBL, 2*M+( M-1 )*
! 487: $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
! 488: WRKBL = MAX( WRKBL, 2*M+M*
! 489: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
! 490: MAXWRK = 2*M*M + WRKBL
! 491: MINWRK = 2*M + N
! 492: ELSE IF( WNTVA .AND. WNTUAS ) THEN
! 493: *
! 494: * Path 9t(N much larger than M, JOBU='S' or 'A',
! 495: * JOBVT='A')
! 496: *
! 497: WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 498: WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N,
! 499: $ N, M, -1 ) )
! 500: WRKBL = MAX( WRKBL, 2*M+2*M*
! 501: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 502: WRKBL = MAX( WRKBL, 2*M+( M-1 )*
! 503: $ ILAENV( 1, 'ZUNGBR', 'P', M, M, M, -1 ) )
! 504: WRKBL = MAX( WRKBL, 2*M+M*
! 505: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
! 506: MAXWRK = M*M + WRKBL
! 507: MINWRK = 2*M + N
! 508: END IF
! 509: ELSE
! 510: *
! 511: * Path 10t(N greater than M, but not much larger)
! 512: *
! 513: MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,
! 514: $ -1, -1 )
! 515: IF( WNTVS .OR. WNTVO )
! 516: $ MAXWRK = MAX( MAXWRK, 2*M+M*
! 517: $ ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) )
! 518: IF( WNTVA )
! 519: $ MAXWRK = MAX( MAXWRK, 2*M+N*
! 520: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, M, -1 ) )
! 521: IF( .NOT.WNTUN )
! 522: $ MAXWRK = MAX( MAXWRK, 2*M+( M-1 )*
! 523: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, M, -1 ) )
! 524: MINWRK = 2*M + N
! 525: END IF
! 526: END IF
! 527: MAXWRK = MAX( MAXWRK, MINWRK )
! 528: WORK( 1 ) = MAXWRK
! 529: *
! 530: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
! 531: INFO = -13
! 532: END IF
! 533: END IF
! 534: *
! 535: IF( INFO.NE.0 ) THEN
! 536: CALL XERBLA( 'ZGESVD', -INFO )
! 537: RETURN
! 538: ELSE IF( LQUERY ) THEN
! 539: RETURN
! 540: END IF
! 541: *
! 542: * Quick return if possible
! 543: *
! 544: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
! 545: RETURN
! 546: END IF
! 547: *
! 548: * Get machine constants
! 549: *
! 550: EPS = DLAMCH( 'P' )
! 551: SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
! 552: BIGNUM = ONE / SMLNUM
! 553: *
! 554: * Scale A if max element outside range [SMLNUM,BIGNUM]
! 555: *
! 556: ANRM = ZLANGE( 'M', M, N, A, LDA, DUM )
! 557: ISCL = 0
! 558: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 559: ISCL = 1
! 560: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
! 561: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 562: ISCL = 1
! 563: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
! 564: END IF
! 565: *
! 566: IF( M.GE.N ) THEN
! 567: *
! 568: * A has at least as many rows as columns. If A has sufficiently
! 569: * more rows than columns, first reduce using the QR
! 570: * decomposition (if sufficient workspace available)
! 571: *
! 572: IF( M.GE.MNTHR ) THEN
! 573: *
! 574: IF( WNTUN ) THEN
! 575: *
! 576: * Path 1 (M much larger than N, JOBU='N')
! 577: * No left singular vectors to be computed
! 578: *
! 579: ITAU = 1
! 580: IWORK = ITAU + N
! 581: *
! 582: * Compute A=Q*R
! 583: * (CWorkspace: need 2*N, prefer N+N*NB)
! 584: * (RWorkspace: need 0)
! 585: *
! 586: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
! 587: $ LWORK-IWORK+1, IERR )
! 588: *
! 589: * Zero out below R
! 590: *
! 591: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
! 592: $ LDA )
! 593: IE = 1
! 594: ITAUQ = 1
! 595: ITAUP = ITAUQ + N
! 596: IWORK = ITAUP + N
! 597: *
! 598: * Bidiagonalize R in A
! 599: * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
! 600: * (RWorkspace: need N)
! 601: *
! 602: CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 603: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
! 604: $ IERR )
! 605: NCVT = 0
! 606: IF( WNTVO .OR. WNTVAS ) THEN
! 607: *
! 608: * If right singular vectors desired, generate P'.
! 609: * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
! 610: * (RWorkspace: 0)
! 611: *
! 612: CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
! 613: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 614: NCVT = N
! 615: END IF
! 616: IRWORK = IE + N
! 617: *
! 618: * Perform bidiagonal QR iteration, computing right
! 619: * singular vectors of A in A if desired
! 620: * (CWorkspace: 0)
! 621: * (RWorkspace: need BDSPAC)
! 622: *
! 623: CALL ZBDSQR( 'U', N, NCVT, 0, 0, S, RWORK( IE ), A, LDA,
! 624: $ CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO )
! 625: *
! 626: * If right singular vectors desired in VT, copy them there
! 627: *
! 628: IF( WNTVAS )
! 629: $ CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT )
! 630: *
! 631: ELSE IF( WNTUO .AND. WNTVN ) THEN
! 632: *
! 633: * Path 2 (M much larger than N, JOBU='O', JOBVT='N')
! 634: * N left singular vectors to be overwritten on A and
! 635: * no right singular vectors to be computed
! 636: *
! 637: IF( LWORK.GE.N*N+3*N ) THEN
! 638: *
! 639: * Sufficient workspace for a fast algorithm
! 640: *
! 641: IR = 1
! 642: IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN
! 643: *
! 644: * WORK(IU) is LDA by N, WORK(IR) is LDA by N
! 645: *
! 646: LDWRKU = LDA
! 647: LDWRKR = LDA
! 648: ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN
! 649: *
! 650: * WORK(IU) is LDA by N, WORK(IR) is N by N
! 651: *
! 652: LDWRKU = LDA
! 653: LDWRKR = N
! 654: ELSE
! 655: *
! 656: * WORK(IU) is LDWRKU by N, WORK(IR) is N by N
! 657: *
! 658: LDWRKU = ( LWORK-N*N ) / N
! 659: LDWRKR = N
! 660: END IF
! 661: ITAU = IR + LDWRKR*N
! 662: IWORK = ITAU + N
! 663: *
! 664: * Compute A=Q*R
! 665: * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
! 666: * (RWorkspace: 0)
! 667: *
! 668: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 669: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 670: *
! 671: * Copy R to WORK(IR) and zero out below it
! 672: *
! 673: CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
! 674: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 675: $ WORK( IR+1 ), LDWRKR )
! 676: *
! 677: * Generate Q in A
! 678: * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
! 679: * (RWorkspace: 0)
! 680: *
! 681: CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
! 682: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 683: IE = 1
! 684: ITAUQ = ITAU
! 685: ITAUP = ITAUQ + N
! 686: IWORK = ITAUP + N
! 687: *
! 688: * Bidiagonalize R in WORK(IR)
! 689: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
! 690: * (RWorkspace: need N)
! 691: *
! 692: CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
! 693: $ WORK( ITAUQ ), WORK( ITAUP ),
! 694: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 695: *
! 696: * Generate left vectors bidiagonalizing R
! 697: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
! 698: * (RWorkspace: need 0)
! 699: *
! 700: CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
! 701: $ WORK( ITAUQ ), WORK( IWORK ),
! 702: $ LWORK-IWORK+1, IERR )
! 703: IRWORK = IE + N
! 704: *
! 705: * Perform bidiagonal QR iteration, computing left
! 706: * singular vectors of R in WORK(IR)
! 707: * (CWorkspace: need N*N)
! 708: * (RWorkspace: need BDSPAC)
! 709: *
! 710: CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM, 1,
! 711: $ WORK( IR ), LDWRKR, CDUM, 1,
! 712: $ RWORK( IRWORK ), INFO )
! 713: IU = ITAUQ
! 714: *
! 715: * Multiply Q in A by left singular vectors of R in
! 716: * WORK(IR), storing result in WORK(IU) and copying to A
! 717: * (CWorkspace: need N*N+N, prefer N*N+M*N)
! 718: * (RWorkspace: 0)
! 719: *
! 720: DO 10 I = 1, M, LDWRKU
! 721: CHUNK = MIN( M-I+1, LDWRKU )
! 722: CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
! 723: $ LDA, WORK( IR ), LDWRKR, CZERO,
! 724: $ WORK( IU ), LDWRKU )
! 725: CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
! 726: $ A( I, 1 ), LDA )
! 727: 10 CONTINUE
! 728: *
! 729: ELSE
! 730: *
! 731: * Insufficient workspace for a fast algorithm
! 732: *
! 733: IE = 1
! 734: ITAUQ = 1
! 735: ITAUP = ITAUQ + N
! 736: IWORK = ITAUP + N
! 737: *
! 738: * Bidiagonalize A
! 739: * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
! 740: * (RWorkspace: N)
! 741: *
! 742: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ),
! 743: $ WORK( ITAUQ ), WORK( ITAUP ),
! 744: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 745: *
! 746: * Generate left vectors bidiagonalizing A
! 747: * (CWorkspace: need 3*N, prefer 2*N+N*NB)
! 748: * (RWorkspace: 0)
! 749: *
! 750: CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
! 751: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 752: IRWORK = IE + N
! 753: *
! 754: * Perform bidiagonal QR iteration, computing left
! 755: * singular vectors of A in A
! 756: * (CWorkspace: need 0)
! 757: * (RWorkspace: need BDSPAC)
! 758: *
! 759: CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM, 1,
! 760: $ A, LDA, CDUM, 1, RWORK( IRWORK ), INFO )
! 761: *
! 762: END IF
! 763: *
! 764: ELSE IF( WNTUO .AND. WNTVAS ) THEN
! 765: *
! 766: * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
! 767: * N left singular vectors to be overwritten on A and
! 768: * N right singular vectors to be computed in VT
! 769: *
! 770: IF( LWORK.GE.N*N+3*N ) THEN
! 771: *
! 772: * Sufficient workspace for a fast algorithm
! 773: *
! 774: IR = 1
! 775: IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*N ) THEN
! 776: *
! 777: * WORK(IU) is LDA by N and WORK(IR) is LDA by N
! 778: *
! 779: LDWRKU = LDA
! 780: LDWRKR = LDA
! 781: ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+N*N ) THEN
! 782: *
! 783: * WORK(IU) is LDA by N and WORK(IR) is N by N
! 784: *
! 785: LDWRKU = LDA
! 786: LDWRKR = N
! 787: ELSE
! 788: *
! 789: * WORK(IU) is LDWRKU by N and WORK(IR) is N by N
! 790: *
! 791: LDWRKU = ( LWORK-N*N ) / N
! 792: LDWRKR = N
! 793: END IF
! 794: ITAU = IR + LDWRKR*N
! 795: IWORK = ITAU + N
! 796: *
! 797: * Compute A=Q*R
! 798: * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
! 799: * (RWorkspace: 0)
! 800: *
! 801: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 802: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 803: *
! 804: * Copy R to VT, zeroing out below it
! 805: *
! 806: CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
! 807: IF( N.GT.1 )
! 808: $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 809: $ VT( 2, 1 ), LDVT )
! 810: *
! 811: * Generate Q in A
! 812: * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
! 813: * (RWorkspace: 0)
! 814: *
! 815: CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
! 816: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 817: IE = 1
! 818: ITAUQ = ITAU
! 819: ITAUP = ITAUQ + N
! 820: IWORK = ITAUP + N
! 821: *
! 822: * Bidiagonalize R in VT, copying result to WORK(IR)
! 823: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
! 824: * (RWorkspace: need N)
! 825: *
! 826: CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
! 827: $ WORK( ITAUQ ), WORK( ITAUP ),
! 828: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 829: CALL ZLACPY( 'L', N, N, VT, LDVT, WORK( IR ), LDWRKR )
! 830: *
! 831: * Generate left vectors bidiagonalizing R in WORK(IR)
! 832: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
! 833: * (RWorkspace: 0)
! 834: *
! 835: CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
! 836: $ WORK( ITAUQ ), WORK( IWORK ),
! 837: $ LWORK-IWORK+1, IERR )
! 838: *
! 839: * Generate right vectors bidiagonalizing R in VT
! 840: * (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB)
! 841: * (RWorkspace: 0)
! 842: *
! 843: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
! 844: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 845: IRWORK = IE + N
! 846: *
! 847: * Perform bidiagonal QR iteration, computing left
! 848: * singular vectors of R in WORK(IR) and computing right
! 849: * singular vectors of R in VT
! 850: * (CWorkspace: need N*N)
! 851: * (RWorkspace: need BDSPAC)
! 852: *
! 853: CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
! 854: $ LDVT, WORK( IR ), LDWRKR, CDUM, 1,
! 855: $ RWORK( IRWORK ), INFO )
! 856: IU = ITAUQ
! 857: *
! 858: * Multiply Q in A by left singular vectors of R in
! 859: * WORK(IR), storing result in WORK(IU) and copying to A
! 860: * (CWorkspace: need N*N+N, prefer N*N+M*N)
! 861: * (RWorkspace: 0)
! 862: *
! 863: DO 20 I = 1, M, LDWRKU
! 864: CHUNK = MIN( M-I+1, LDWRKU )
! 865: CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
! 866: $ LDA, WORK( IR ), LDWRKR, CZERO,
! 867: $ WORK( IU ), LDWRKU )
! 868: CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
! 869: $ A( I, 1 ), LDA )
! 870: 20 CONTINUE
! 871: *
! 872: ELSE
! 873: *
! 874: * Insufficient workspace for a fast algorithm
! 875: *
! 876: ITAU = 1
! 877: IWORK = ITAU + N
! 878: *
! 879: * Compute A=Q*R
! 880: * (CWorkspace: need 2*N, prefer N+N*NB)
! 881: * (RWorkspace: 0)
! 882: *
! 883: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 884: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 885: *
! 886: * Copy R to VT, zeroing out below it
! 887: *
! 888: CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
! 889: IF( N.GT.1 )
! 890: $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 891: $ VT( 2, 1 ), LDVT )
! 892: *
! 893: * Generate Q in A
! 894: * (CWorkspace: need 2*N, prefer N+N*NB)
! 895: * (RWorkspace: 0)
! 896: *
! 897: CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
! 898: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 899: IE = 1
! 900: ITAUQ = ITAU
! 901: ITAUP = ITAUQ + N
! 902: IWORK = ITAUP + N
! 903: *
! 904: * Bidiagonalize R in VT
! 905: * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
! 906: * (RWorkspace: N)
! 907: *
! 908: CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
! 909: $ WORK( ITAUQ ), WORK( ITAUP ),
! 910: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 911: *
! 912: * Multiply Q in A by left vectors bidiagonalizing R
! 913: * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
! 914: * (RWorkspace: 0)
! 915: *
! 916: CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
! 917: $ WORK( ITAUQ ), A, LDA, WORK( IWORK ),
! 918: $ LWORK-IWORK+1, IERR )
! 919: *
! 920: * Generate right vectors bidiagonalizing R in VT
! 921: * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
! 922: * (RWorkspace: 0)
! 923: *
! 924: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
! 925: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 926: IRWORK = IE + N
! 927: *
! 928: * Perform bidiagonal QR iteration, computing left
! 929: * singular vectors of A in A and computing right
! 930: * singular vectors of A in VT
! 931: * (CWorkspace: 0)
! 932: * (RWorkspace: need BDSPAC)
! 933: *
! 934: CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
! 935: $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
! 936: $ INFO )
! 937: *
! 938: END IF
! 939: *
! 940: ELSE IF( WNTUS ) THEN
! 941: *
! 942: IF( WNTVN ) THEN
! 943: *
! 944: * Path 4 (M much larger than N, JOBU='S', JOBVT='N')
! 945: * N left singular vectors to be computed in U and
! 946: * no right singular vectors to be computed
! 947: *
! 948: IF( LWORK.GE.N*N+3*N ) THEN
! 949: *
! 950: * Sufficient workspace for a fast algorithm
! 951: *
! 952: IR = 1
! 953: IF( LWORK.GE.WRKBL+LDA*N ) THEN
! 954: *
! 955: * WORK(IR) is LDA by N
! 956: *
! 957: LDWRKR = LDA
! 958: ELSE
! 959: *
! 960: * WORK(IR) is N by N
! 961: *
! 962: LDWRKR = N
! 963: END IF
! 964: ITAU = IR + LDWRKR*N
! 965: IWORK = ITAU + N
! 966: *
! 967: * Compute A=Q*R
! 968: * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
! 969: * (RWorkspace: 0)
! 970: *
! 971: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 972: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 973: *
! 974: * Copy R to WORK(IR), zeroing out below it
! 975: *
! 976: CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ),
! 977: $ LDWRKR )
! 978: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 979: $ WORK( IR+1 ), LDWRKR )
! 980: *
! 981: * Generate Q in A
! 982: * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
! 983: * (RWorkspace: 0)
! 984: *
! 985: CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
! 986: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 987: IE = 1
! 988: ITAUQ = ITAU
! 989: ITAUP = ITAUQ + N
! 990: IWORK = ITAUP + N
! 991: *
! 992: * Bidiagonalize R in WORK(IR)
! 993: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
! 994: * (RWorkspace: need N)
! 995: *
! 996: CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S,
! 997: $ RWORK( IE ), WORK( ITAUQ ),
! 998: $ WORK( ITAUP ), WORK( IWORK ),
! 999: $ LWORK-IWORK+1, IERR )
! 1000: *
! 1001: * Generate left vectors bidiagonalizing R in WORK(IR)
! 1002: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
! 1003: * (RWorkspace: 0)
! 1004: *
! 1005: CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
! 1006: $ WORK( ITAUQ ), WORK( IWORK ),
! 1007: $ LWORK-IWORK+1, IERR )
! 1008: IRWORK = IE + N
! 1009: *
! 1010: * Perform bidiagonal QR iteration, computing left
! 1011: * singular vectors of R in WORK(IR)
! 1012: * (CWorkspace: need N*N)
! 1013: * (RWorkspace: need BDSPAC)
! 1014: *
! 1015: CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM,
! 1016: $ 1, WORK( IR ), LDWRKR, CDUM, 1,
! 1017: $ RWORK( IRWORK ), INFO )
! 1018: *
! 1019: * Multiply Q in A by left singular vectors of R in
! 1020: * WORK(IR), storing result in U
! 1021: * (CWorkspace: need N*N)
! 1022: * (RWorkspace: 0)
! 1023: *
! 1024: CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
! 1025: $ WORK( IR ), LDWRKR, CZERO, U, LDU )
! 1026: *
! 1027: ELSE
! 1028: *
! 1029: * Insufficient workspace for a fast algorithm
! 1030: *
! 1031: ITAU = 1
! 1032: IWORK = ITAU + N
! 1033: *
! 1034: * Compute A=Q*R, copying result to U
! 1035: * (CWorkspace: need 2*N, prefer N+N*NB)
! 1036: * (RWorkspace: 0)
! 1037: *
! 1038: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 1039: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1040: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 1041: *
! 1042: * Generate Q in U
! 1043: * (CWorkspace: need 2*N, prefer N+N*NB)
! 1044: * (RWorkspace: 0)
! 1045: *
! 1046: CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ),
! 1047: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1048: IE = 1
! 1049: ITAUQ = ITAU
! 1050: ITAUP = ITAUQ + N
! 1051: IWORK = ITAUP + N
! 1052: *
! 1053: * Zero out below R in A
! 1054: *
! 1055: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 1056: $ A( 2, 1 ), LDA )
! 1057: *
! 1058: * Bidiagonalize R in A
! 1059: * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
! 1060: * (RWorkspace: need N)
! 1061: *
! 1062: CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
! 1063: $ WORK( ITAUQ ), WORK( ITAUP ),
! 1064: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1065: *
! 1066: * Multiply Q in U by left vectors bidiagonalizing R
! 1067: * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
! 1068: * (RWorkspace: 0)
! 1069: *
! 1070: CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
! 1071: $ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
! 1072: $ LWORK-IWORK+1, IERR )
! 1073: IRWORK = IE + N
! 1074: *
! 1075: * Perform bidiagonal QR iteration, computing left
! 1076: * singular vectors of A in U
! 1077: * (CWorkspace: 0)
! 1078: * (RWorkspace: need BDSPAC)
! 1079: *
! 1080: CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM,
! 1081: $ 1, U, LDU, CDUM, 1, RWORK( IRWORK ),
! 1082: $ INFO )
! 1083: *
! 1084: END IF
! 1085: *
! 1086: ELSE IF( WNTVO ) THEN
! 1087: *
! 1088: * Path 5 (M much larger than N, JOBU='S', JOBVT='O')
! 1089: * N left singular vectors to be computed in U and
! 1090: * N right singular vectors to be overwritten on A
! 1091: *
! 1092: IF( LWORK.GE.2*N*N+3*N ) THEN
! 1093: *
! 1094: * Sufficient workspace for a fast algorithm
! 1095: *
! 1096: IU = 1
! 1097: IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
! 1098: *
! 1099: * WORK(IU) is LDA by N and WORK(IR) is LDA by N
! 1100: *
! 1101: LDWRKU = LDA
! 1102: IR = IU + LDWRKU*N
! 1103: LDWRKR = LDA
! 1104: ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
! 1105: *
! 1106: * WORK(IU) is LDA by N and WORK(IR) is N by N
! 1107: *
! 1108: LDWRKU = LDA
! 1109: IR = IU + LDWRKU*N
! 1110: LDWRKR = N
! 1111: ELSE
! 1112: *
! 1113: * WORK(IU) is N by N and WORK(IR) is N by N
! 1114: *
! 1115: LDWRKU = N
! 1116: IR = IU + LDWRKU*N
! 1117: LDWRKR = N
! 1118: END IF
! 1119: ITAU = IR + LDWRKR*N
! 1120: IWORK = ITAU + N
! 1121: *
! 1122: * Compute A=Q*R
! 1123: * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
! 1124: * (RWorkspace: 0)
! 1125: *
! 1126: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 1127: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1128: *
! 1129: * Copy R to WORK(IU), zeroing out below it
! 1130: *
! 1131: CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
! 1132: $ LDWRKU )
! 1133: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 1134: $ WORK( IU+1 ), LDWRKU )
! 1135: *
! 1136: * Generate Q in A
! 1137: * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
! 1138: * (RWorkspace: 0)
! 1139: *
! 1140: CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
! 1141: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1142: IE = 1
! 1143: ITAUQ = ITAU
! 1144: ITAUP = ITAUQ + N
! 1145: IWORK = ITAUP + N
! 1146: *
! 1147: * Bidiagonalize R in WORK(IU), copying result to
! 1148: * WORK(IR)
! 1149: * (CWorkspace: need 2*N*N+3*N,
! 1150: * prefer 2*N*N+2*N+2*N*NB)
! 1151: * (RWorkspace: need N)
! 1152: *
! 1153: CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
! 1154: $ RWORK( IE ), WORK( ITAUQ ),
! 1155: $ WORK( ITAUP ), WORK( IWORK ),
! 1156: $ LWORK-IWORK+1, IERR )
! 1157: CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU,
! 1158: $ WORK( IR ), LDWRKR )
! 1159: *
! 1160: * Generate left bidiagonalizing vectors in WORK(IU)
! 1161: * (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
! 1162: * (RWorkspace: 0)
! 1163: *
! 1164: CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
! 1165: $ WORK( ITAUQ ), WORK( IWORK ),
! 1166: $ LWORK-IWORK+1, IERR )
! 1167: *
! 1168: * Generate right bidiagonalizing vectors in WORK(IR)
! 1169: * (CWorkspace: need 2*N*N+3*N-1,
! 1170: * prefer 2*N*N+2*N+(N-1)*NB)
! 1171: * (RWorkspace: 0)
! 1172: *
! 1173: CALL ZUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
! 1174: $ WORK( ITAUP ), WORK( IWORK ),
! 1175: $ LWORK-IWORK+1, IERR )
! 1176: IRWORK = IE + N
! 1177: *
! 1178: * Perform bidiagonal QR iteration, computing left
! 1179: * singular vectors of R in WORK(IU) and computing
! 1180: * right singular vectors of R in WORK(IR)
! 1181: * (CWorkspace: need 2*N*N)
! 1182: * (RWorkspace: need BDSPAC)
! 1183: *
! 1184: CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ),
! 1185: $ WORK( IR ), LDWRKR, WORK( IU ),
! 1186: $ LDWRKU, CDUM, 1, RWORK( IRWORK ),
! 1187: $ INFO )
! 1188: *
! 1189: * Multiply Q in A by left singular vectors of R in
! 1190: * WORK(IU), storing result in U
! 1191: * (CWorkspace: need N*N)
! 1192: * (RWorkspace: 0)
! 1193: *
! 1194: CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
! 1195: $ WORK( IU ), LDWRKU, CZERO, U, LDU )
! 1196: *
! 1197: * Copy right singular vectors of R to A
! 1198: * (CWorkspace: need N*N)
! 1199: * (RWorkspace: 0)
! 1200: *
! 1201: CALL ZLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
! 1202: $ LDA )
! 1203: *
! 1204: ELSE
! 1205: *
! 1206: * Insufficient workspace for a fast algorithm
! 1207: *
! 1208: ITAU = 1
! 1209: IWORK = ITAU + N
! 1210: *
! 1211: * Compute A=Q*R, copying result to U
! 1212: * (CWorkspace: need 2*N, prefer N+N*NB)
! 1213: * (RWorkspace: 0)
! 1214: *
! 1215: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 1216: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1217: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 1218: *
! 1219: * Generate Q in U
! 1220: * (CWorkspace: need 2*N, prefer N+N*NB)
! 1221: * (RWorkspace: 0)
! 1222: *
! 1223: CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ),
! 1224: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1225: IE = 1
! 1226: ITAUQ = ITAU
! 1227: ITAUP = ITAUQ + N
! 1228: IWORK = ITAUP + N
! 1229: *
! 1230: * Zero out below R in A
! 1231: *
! 1232: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 1233: $ A( 2, 1 ), LDA )
! 1234: *
! 1235: * Bidiagonalize R in A
! 1236: * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
! 1237: * (RWorkspace: need N)
! 1238: *
! 1239: CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
! 1240: $ WORK( ITAUQ ), WORK( ITAUP ),
! 1241: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1242: *
! 1243: * Multiply Q in U by left vectors bidiagonalizing R
! 1244: * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
! 1245: * (RWorkspace: 0)
! 1246: *
! 1247: CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
! 1248: $ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
! 1249: $ LWORK-IWORK+1, IERR )
! 1250: *
! 1251: * Generate right vectors bidiagonalizing R in A
! 1252: * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
! 1253: * (RWorkspace: 0)
! 1254: *
! 1255: CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
! 1256: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1257: IRWORK = IE + N
! 1258: *
! 1259: * Perform bidiagonal QR iteration, computing left
! 1260: * singular vectors of A in U and computing right
! 1261: * singular vectors of A in A
! 1262: * (CWorkspace: 0)
! 1263: * (RWorkspace: need BDSPAC)
! 1264: *
! 1265: CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A,
! 1266: $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
! 1267: $ INFO )
! 1268: *
! 1269: END IF
! 1270: *
! 1271: ELSE IF( WNTVAS ) THEN
! 1272: *
! 1273: * Path 6 (M much larger than N, JOBU='S', JOBVT='S'
! 1274: * or 'A')
! 1275: * N left singular vectors to be computed in U and
! 1276: * N right singular vectors to be computed in VT
! 1277: *
! 1278: IF( LWORK.GE.N*N+3*N ) THEN
! 1279: *
! 1280: * Sufficient workspace for a fast algorithm
! 1281: *
! 1282: IU = 1
! 1283: IF( LWORK.GE.WRKBL+LDA*N ) THEN
! 1284: *
! 1285: * WORK(IU) is LDA by N
! 1286: *
! 1287: LDWRKU = LDA
! 1288: ELSE
! 1289: *
! 1290: * WORK(IU) is N by N
! 1291: *
! 1292: LDWRKU = N
! 1293: END IF
! 1294: ITAU = IU + LDWRKU*N
! 1295: IWORK = ITAU + N
! 1296: *
! 1297: * Compute A=Q*R
! 1298: * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
! 1299: * (RWorkspace: 0)
! 1300: *
! 1301: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 1302: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1303: *
! 1304: * Copy R to WORK(IU), zeroing out below it
! 1305: *
! 1306: CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
! 1307: $ LDWRKU )
! 1308: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 1309: $ WORK( IU+1 ), LDWRKU )
! 1310: *
! 1311: * Generate Q in A
! 1312: * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
! 1313: * (RWorkspace: 0)
! 1314: *
! 1315: CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
! 1316: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1317: IE = 1
! 1318: ITAUQ = ITAU
! 1319: ITAUP = ITAUQ + N
! 1320: IWORK = ITAUP + N
! 1321: *
! 1322: * Bidiagonalize R in WORK(IU), copying result to VT
! 1323: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
! 1324: * (RWorkspace: need N)
! 1325: *
! 1326: CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
! 1327: $ RWORK( IE ), WORK( ITAUQ ),
! 1328: $ WORK( ITAUP ), WORK( IWORK ),
! 1329: $ LWORK-IWORK+1, IERR )
! 1330: CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
! 1331: $ LDVT )
! 1332: *
! 1333: * Generate left bidiagonalizing vectors in WORK(IU)
! 1334: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
! 1335: * (RWorkspace: 0)
! 1336: *
! 1337: CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
! 1338: $ WORK( ITAUQ ), WORK( IWORK ),
! 1339: $ LWORK-IWORK+1, IERR )
! 1340: *
! 1341: * Generate right bidiagonalizing vectors in VT
! 1342: * (CWorkspace: need N*N+3*N-1,
! 1343: * prefer N*N+2*N+(N-1)*NB)
! 1344: * (RWorkspace: 0)
! 1345: *
! 1346: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
! 1347: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1348: IRWORK = IE + N
! 1349: *
! 1350: * Perform bidiagonal QR iteration, computing left
! 1351: * singular vectors of R in WORK(IU) and computing
! 1352: * right singular vectors of R in VT
! 1353: * (CWorkspace: need N*N)
! 1354: * (RWorkspace: need BDSPAC)
! 1355: *
! 1356: CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
! 1357: $ LDVT, WORK( IU ), LDWRKU, CDUM, 1,
! 1358: $ RWORK( IRWORK ), INFO )
! 1359: *
! 1360: * Multiply Q in A by left singular vectors of R in
! 1361: * WORK(IU), storing result in U
! 1362: * (CWorkspace: need N*N)
! 1363: * (RWorkspace: 0)
! 1364: *
! 1365: CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA,
! 1366: $ WORK( IU ), LDWRKU, CZERO, U, LDU )
! 1367: *
! 1368: ELSE
! 1369: *
! 1370: * Insufficient workspace for a fast algorithm
! 1371: *
! 1372: ITAU = 1
! 1373: IWORK = ITAU + N
! 1374: *
! 1375: * Compute A=Q*R, copying result to U
! 1376: * (CWorkspace: need 2*N, prefer N+N*NB)
! 1377: * (RWorkspace: 0)
! 1378: *
! 1379: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 1380: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1381: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 1382: *
! 1383: * Generate Q in U
! 1384: * (CWorkspace: need 2*N, prefer N+N*NB)
! 1385: * (RWorkspace: 0)
! 1386: *
! 1387: CALL ZUNGQR( M, N, N, U, LDU, WORK( ITAU ),
! 1388: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1389: *
! 1390: * Copy R to VT, zeroing out below it
! 1391: *
! 1392: CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
! 1393: IF( N.GT.1 )
! 1394: $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 1395: $ VT( 2, 1 ), LDVT )
! 1396: IE = 1
! 1397: ITAUQ = ITAU
! 1398: ITAUP = ITAUQ + N
! 1399: IWORK = ITAUP + N
! 1400: *
! 1401: * Bidiagonalize R in VT
! 1402: * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
! 1403: * (RWorkspace: need N)
! 1404: *
! 1405: CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
! 1406: $ WORK( ITAUQ ), WORK( ITAUP ),
! 1407: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1408: *
! 1409: * Multiply Q in U by left bidiagonalizing vectors
! 1410: * in VT
! 1411: * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
! 1412: * (RWorkspace: 0)
! 1413: *
! 1414: CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
! 1415: $ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
! 1416: $ LWORK-IWORK+1, IERR )
! 1417: *
! 1418: * Generate right bidiagonalizing vectors in VT
! 1419: * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
! 1420: * (RWorkspace: 0)
! 1421: *
! 1422: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
! 1423: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1424: IRWORK = IE + N
! 1425: *
! 1426: * Perform bidiagonal QR iteration, computing left
! 1427: * singular vectors of A in U and computing right
! 1428: * singular vectors of A in VT
! 1429: * (CWorkspace: 0)
! 1430: * (RWorkspace: need BDSPAC)
! 1431: *
! 1432: CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
! 1433: $ LDVT, U, LDU, CDUM, 1,
! 1434: $ RWORK( IRWORK ), INFO )
! 1435: *
! 1436: END IF
! 1437: *
! 1438: END IF
! 1439: *
! 1440: ELSE IF( WNTUA ) THEN
! 1441: *
! 1442: IF( WNTVN ) THEN
! 1443: *
! 1444: * Path 7 (M much larger than N, JOBU='A', JOBVT='N')
! 1445: * M left singular vectors to be computed in U and
! 1446: * no right singular vectors to be computed
! 1447: *
! 1448: IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN
! 1449: *
! 1450: * Sufficient workspace for a fast algorithm
! 1451: *
! 1452: IR = 1
! 1453: IF( LWORK.GE.WRKBL+LDA*N ) THEN
! 1454: *
! 1455: * WORK(IR) is LDA by N
! 1456: *
! 1457: LDWRKR = LDA
! 1458: ELSE
! 1459: *
! 1460: * WORK(IR) is N by N
! 1461: *
! 1462: LDWRKR = N
! 1463: END IF
! 1464: ITAU = IR + LDWRKR*N
! 1465: IWORK = ITAU + N
! 1466: *
! 1467: * Compute A=Q*R, copying result to U
! 1468: * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
! 1469: * (RWorkspace: 0)
! 1470: *
! 1471: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 1472: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1473: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 1474: *
! 1475: * Copy R to WORK(IR), zeroing out below it
! 1476: *
! 1477: CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ),
! 1478: $ LDWRKR )
! 1479: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 1480: $ WORK( IR+1 ), LDWRKR )
! 1481: *
! 1482: * Generate Q in U
! 1483: * (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
! 1484: * (RWorkspace: 0)
! 1485: *
! 1486: CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
! 1487: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1488: IE = 1
! 1489: ITAUQ = ITAU
! 1490: ITAUP = ITAUQ + N
! 1491: IWORK = ITAUP + N
! 1492: *
! 1493: * Bidiagonalize R in WORK(IR)
! 1494: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
! 1495: * (RWorkspace: need N)
! 1496: *
! 1497: CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S,
! 1498: $ RWORK( IE ), WORK( ITAUQ ),
! 1499: $ WORK( ITAUP ), WORK( IWORK ),
! 1500: $ LWORK-IWORK+1, IERR )
! 1501: *
! 1502: * Generate left bidiagonalizing vectors in WORK(IR)
! 1503: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
! 1504: * (RWorkspace: 0)
! 1505: *
! 1506: CALL ZUNGBR( 'Q', N, N, N, WORK( IR ), LDWRKR,
! 1507: $ WORK( ITAUQ ), WORK( IWORK ),
! 1508: $ LWORK-IWORK+1, IERR )
! 1509: IRWORK = IE + N
! 1510: *
! 1511: * Perform bidiagonal QR iteration, computing left
! 1512: * singular vectors of R in WORK(IR)
! 1513: * (CWorkspace: need N*N)
! 1514: * (RWorkspace: need BDSPAC)
! 1515: *
! 1516: CALL ZBDSQR( 'U', N, 0, N, 0, S, RWORK( IE ), CDUM,
! 1517: $ 1, WORK( IR ), LDWRKR, CDUM, 1,
! 1518: $ RWORK( IRWORK ), INFO )
! 1519: *
! 1520: * Multiply Q in U by left singular vectors of R in
! 1521: * WORK(IR), storing result in A
! 1522: * (CWorkspace: need N*N)
! 1523: * (RWorkspace: 0)
! 1524: *
! 1525: CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
! 1526: $ WORK( IR ), LDWRKR, CZERO, A, LDA )
! 1527: *
! 1528: * Copy left singular vectors of A from A to U
! 1529: *
! 1530: CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
! 1531: *
! 1532: ELSE
! 1533: *
! 1534: * Insufficient workspace for a fast algorithm
! 1535: *
! 1536: ITAU = 1
! 1537: IWORK = ITAU + N
! 1538: *
! 1539: * Compute A=Q*R, copying result to U
! 1540: * (CWorkspace: need 2*N, prefer N+N*NB)
! 1541: * (RWorkspace: 0)
! 1542: *
! 1543: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 1544: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1545: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 1546: *
! 1547: * Generate Q in U
! 1548: * (CWorkspace: need N+M, prefer N+M*NB)
! 1549: * (RWorkspace: 0)
! 1550: *
! 1551: CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
! 1552: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1553: IE = 1
! 1554: ITAUQ = ITAU
! 1555: ITAUP = ITAUQ + N
! 1556: IWORK = ITAUP + N
! 1557: *
! 1558: * Zero out below R in A
! 1559: *
! 1560: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 1561: $ A( 2, 1 ), LDA )
! 1562: *
! 1563: * Bidiagonalize R in A
! 1564: * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
! 1565: * (RWorkspace: need N)
! 1566: *
! 1567: CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
! 1568: $ WORK( ITAUQ ), WORK( ITAUP ),
! 1569: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1570: *
! 1571: * Multiply Q in U by left bidiagonalizing vectors
! 1572: * in A
! 1573: * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
! 1574: * (RWorkspace: 0)
! 1575: *
! 1576: CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
! 1577: $ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
! 1578: $ LWORK-IWORK+1, IERR )
! 1579: IRWORK = IE + N
! 1580: *
! 1581: * Perform bidiagonal QR iteration, computing left
! 1582: * singular vectors of A in U
! 1583: * (CWorkspace: 0)
! 1584: * (RWorkspace: need BDSPAC)
! 1585: *
! 1586: CALL ZBDSQR( 'U', N, 0, M, 0, S, RWORK( IE ), CDUM,
! 1587: $ 1, U, LDU, CDUM, 1, RWORK( IRWORK ),
! 1588: $ INFO )
! 1589: *
! 1590: END IF
! 1591: *
! 1592: ELSE IF( WNTVO ) THEN
! 1593: *
! 1594: * Path 8 (M much larger than N, JOBU='A', JOBVT='O')
! 1595: * M left singular vectors to be computed in U and
! 1596: * N right singular vectors to be overwritten on A
! 1597: *
! 1598: IF( LWORK.GE.2*N*N+MAX( N+M, 3*N ) ) THEN
! 1599: *
! 1600: * Sufficient workspace for a fast algorithm
! 1601: *
! 1602: IU = 1
! 1603: IF( LWORK.GE.WRKBL+2*LDA*N ) THEN
! 1604: *
! 1605: * WORK(IU) is LDA by N and WORK(IR) is LDA by N
! 1606: *
! 1607: LDWRKU = LDA
! 1608: IR = IU + LDWRKU*N
! 1609: LDWRKR = LDA
! 1610: ELSE IF( LWORK.GE.WRKBL+( LDA+N )*N ) THEN
! 1611: *
! 1612: * WORK(IU) is LDA by N and WORK(IR) is N by N
! 1613: *
! 1614: LDWRKU = LDA
! 1615: IR = IU + LDWRKU*N
! 1616: LDWRKR = N
! 1617: ELSE
! 1618: *
! 1619: * WORK(IU) is N by N and WORK(IR) is N by N
! 1620: *
! 1621: LDWRKU = N
! 1622: IR = IU + LDWRKU*N
! 1623: LDWRKR = N
! 1624: END IF
! 1625: ITAU = IR + LDWRKR*N
! 1626: IWORK = ITAU + N
! 1627: *
! 1628: * Compute A=Q*R, copying result to U
! 1629: * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
! 1630: * (RWorkspace: 0)
! 1631: *
! 1632: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 1633: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1634: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 1635: *
! 1636: * Generate Q in U
! 1637: * (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
! 1638: * (RWorkspace: 0)
! 1639: *
! 1640: CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
! 1641: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1642: *
! 1643: * Copy R to WORK(IU), zeroing out below it
! 1644: *
! 1645: CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
! 1646: $ LDWRKU )
! 1647: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 1648: $ WORK( IU+1 ), LDWRKU )
! 1649: IE = 1
! 1650: ITAUQ = ITAU
! 1651: ITAUP = ITAUQ + N
! 1652: IWORK = ITAUP + N
! 1653: *
! 1654: * Bidiagonalize R in WORK(IU), copying result to
! 1655: * WORK(IR)
! 1656: * (CWorkspace: need 2*N*N+3*N,
! 1657: * prefer 2*N*N+2*N+2*N*NB)
! 1658: * (RWorkspace: need N)
! 1659: *
! 1660: CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
! 1661: $ RWORK( IE ), WORK( ITAUQ ),
! 1662: $ WORK( ITAUP ), WORK( IWORK ),
! 1663: $ LWORK-IWORK+1, IERR )
! 1664: CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU,
! 1665: $ WORK( IR ), LDWRKR )
! 1666: *
! 1667: * Generate left bidiagonalizing vectors in WORK(IU)
! 1668: * (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
! 1669: * (RWorkspace: 0)
! 1670: *
! 1671: CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
! 1672: $ WORK( ITAUQ ), WORK( IWORK ),
! 1673: $ LWORK-IWORK+1, IERR )
! 1674: *
! 1675: * Generate right bidiagonalizing vectors in WORK(IR)
! 1676: * (CWorkspace: need 2*N*N+3*N-1,
! 1677: * prefer 2*N*N+2*N+(N-1)*NB)
! 1678: * (RWorkspace: 0)
! 1679: *
! 1680: CALL ZUNGBR( 'P', N, N, N, WORK( IR ), LDWRKR,
! 1681: $ WORK( ITAUP ), WORK( IWORK ),
! 1682: $ LWORK-IWORK+1, IERR )
! 1683: IRWORK = IE + N
! 1684: *
! 1685: * Perform bidiagonal QR iteration, computing left
! 1686: * singular vectors of R in WORK(IU) and computing
! 1687: * right singular vectors of R in WORK(IR)
! 1688: * (CWorkspace: need 2*N*N)
! 1689: * (RWorkspace: need BDSPAC)
! 1690: *
! 1691: CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ),
! 1692: $ WORK( IR ), LDWRKR, WORK( IU ),
! 1693: $ LDWRKU, CDUM, 1, RWORK( IRWORK ),
! 1694: $ INFO )
! 1695: *
! 1696: * Multiply Q in U by left singular vectors of R in
! 1697: * WORK(IU), storing result in A
! 1698: * (CWorkspace: need N*N)
! 1699: * (RWorkspace: 0)
! 1700: *
! 1701: CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
! 1702: $ WORK( IU ), LDWRKU, CZERO, A, LDA )
! 1703: *
! 1704: * Copy left singular vectors of A from A to U
! 1705: *
! 1706: CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
! 1707: *
! 1708: * Copy right singular vectors of R from WORK(IR) to A
! 1709: *
! 1710: CALL ZLACPY( 'F', N, N, WORK( IR ), LDWRKR, A,
! 1711: $ LDA )
! 1712: *
! 1713: ELSE
! 1714: *
! 1715: * Insufficient workspace for a fast algorithm
! 1716: *
! 1717: ITAU = 1
! 1718: IWORK = ITAU + N
! 1719: *
! 1720: * Compute A=Q*R, copying result to U
! 1721: * (CWorkspace: need 2*N, prefer N+N*NB)
! 1722: * (RWorkspace: 0)
! 1723: *
! 1724: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 1725: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1726: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 1727: *
! 1728: * Generate Q in U
! 1729: * (CWorkspace: need N+M, prefer N+M*NB)
! 1730: * (RWorkspace: 0)
! 1731: *
! 1732: CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
! 1733: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1734: IE = 1
! 1735: ITAUQ = ITAU
! 1736: ITAUP = ITAUQ + N
! 1737: IWORK = ITAUP + N
! 1738: *
! 1739: * Zero out below R in A
! 1740: *
! 1741: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 1742: $ A( 2, 1 ), LDA )
! 1743: *
! 1744: * Bidiagonalize R in A
! 1745: * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
! 1746: * (RWorkspace: need N)
! 1747: *
! 1748: CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ),
! 1749: $ WORK( ITAUQ ), WORK( ITAUP ),
! 1750: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1751: *
! 1752: * Multiply Q in U by left bidiagonalizing vectors
! 1753: * in A
! 1754: * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
! 1755: * (RWorkspace: 0)
! 1756: *
! 1757: CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, A, LDA,
! 1758: $ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
! 1759: $ LWORK-IWORK+1, IERR )
! 1760: *
! 1761: * Generate right bidiagonalizing vectors in A
! 1762: * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
! 1763: * (RWorkspace: 0)
! 1764: *
! 1765: CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
! 1766: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1767: IRWORK = IE + N
! 1768: *
! 1769: * Perform bidiagonal QR iteration, computing left
! 1770: * singular vectors of A in U and computing right
! 1771: * singular vectors of A in A
! 1772: * (CWorkspace: 0)
! 1773: * (RWorkspace: need BDSPAC)
! 1774: *
! 1775: CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), A,
! 1776: $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
! 1777: $ INFO )
! 1778: *
! 1779: END IF
! 1780: *
! 1781: ELSE IF( WNTVAS ) THEN
! 1782: *
! 1783: * Path 9 (M much larger than N, JOBU='A', JOBVT='S'
! 1784: * or 'A')
! 1785: * M left singular vectors to be computed in U and
! 1786: * N right singular vectors to be computed in VT
! 1787: *
! 1788: IF( LWORK.GE.N*N+MAX( N+M, 3*N ) ) THEN
! 1789: *
! 1790: * Sufficient workspace for a fast algorithm
! 1791: *
! 1792: IU = 1
! 1793: IF( LWORK.GE.WRKBL+LDA*N ) THEN
! 1794: *
! 1795: * WORK(IU) is LDA by N
! 1796: *
! 1797: LDWRKU = LDA
! 1798: ELSE
! 1799: *
! 1800: * WORK(IU) is N by N
! 1801: *
! 1802: LDWRKU = N
! 1803: END IF
! 1804: ITAU = IU + LDWRKU*N
! 1805: IWORK = ITAU + N
! 1806: *
! 1807: * Compute A=Q*R, copying result to U
! 1808: * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
! 1809: * (RWorkspace: 0)
! 1810: *
! 1811: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 1812: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1813: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 1814: *
! 1815: * Generate Q in U
! 1816: * (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
! 1817: * (RWorkspace: 0)
! 1818: *
! 1819: CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
! 1820: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1821: *
! 1822: * Copy R to WORK(IU), zeroing out below it
! 1823: *
! 1824: CALL ZLACPY( 'U', N, N, A, LDA, WORK( IU ),
! 1825: $ LDWRKU )
! 1826: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 1827: $ WORK( IU+1 ), LDWRKU )
! 1828: IE = 1
! 1829: ITAUQ = ITAU
! 1830: ITAUP = ITAUQ + N
! 1831: IWORK = ITAUP + N
! 1832: *
! 1833: * Bidiagonalize R in WORK(IU), copying result to VT
! 1834: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
! 1835: * (RWorkspace: need N)
! 1836: *
! 1837: CALL ZGEBRD( N, N, WORK( IU ), LDWRKU, S,
! 1838: $ RWORK( IE ), WORK( ITAUQ ),
! 1839: $ WORK( ITAUP ), WORK( IWORK ),
! 1840: $ LWORK-IWORK+1, IERR )
! 1841: CALL ZLACPY( 'U', N, N, WORK( IU ), LDWRKU, VT,
! 1842: $ LDVT )
! 1843: *
! 1844: * Generate left bidiagonalizing vectors in WORK(IU)
! 1845: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
! 1846: * (RWorkspace: 0)
! 1847: *
! 1848: CALL ZUNGBR( 'Q', N, N, N, WORK( IU ), LDWRKU,
! 1849: $ WORK( ITAUQ ), WORK( IWORK ),
! 1850: $ LWORK-IWORK+1, IERR )
! 1851: *
! 1852: * Generate right bidiagonalizing vectors in VT
! 1853: * (CWorkspace: need N*N+3*N-1,
! 1854: * prefer N*N+2*N+(N-1)*NB)
! 1855: * (RWorkspace: need 0)
! 1856: *
! 1857: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
! 1858: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1859: IRWORK = IE + N
! 1860: *
! 1861: * Perform bidiagonal QR iteration, computing left
! 1862: * singular vectors of R in WORK(IU) and computing
! 1863: * right singular vectors of R in VT
! 1864: * (CWorkspace: need N*N)
! 1865: * (RWorkspace: need BDSPAC)
! 1866: *
! 1867: CALL ZBDSQR( 'U', N, N, N, 0, S, RWORK( IE ), VT,
! 1868: $ LDVT, WORK( IU ), LDWRKU, CDUM, 1,
! 1869: $ RWORK( IRWORK ), INFO )
! 1870: *
! 1871: * Multiply Q in U by left singular vectors of R in
! 1872: * WORK(IU), storing result in A
! 1873: * (CWorkspace: need N*N)
! 1874: * (RWorkspace: 0)
! 1875: *
! 1876: CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU,
! 1877: $ WORK( IU ), LDWRKU, CZERO, A, LDA )
! 1878: *
! 1879: * Copy left singular vectors of A from A to U
! 1880: *
! 1881: CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
! 1882: *
! 1883: ELSE
! 1884: *
! 1885: * Insufficient workspace for a fast algorithm
! 1886: *
! 1887: ITAU = 1
! 1888: IWORK = ITAU + N
! 1889: *
! 1890: * Compute A=Q*R, copying result to U
! 1891: * (CWorkspace: need 2*N, prefer N+N*NB)
! 1892: * (RWorkspace: 0)
! 1893: *
! 1894: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ),
! 1895: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1896: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 1897: *
! 1898: * Generate Q in U
! 1899: * (CWorkspace: need N+M, prefer N+M*NB)
! 1900: * (RWorkspace: 0)
! 1901: *
! 1902: CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
! 1903: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1904: *
! 1905: * Copy R from A to VT, zeroing out below it
! 1906: *
! 1907: CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
! 1908: IF( N.GT.1 )
! 1909: $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO,
! 1910: $ VT( 2, 1 ), LDVT )
! 1911: IE = 1
! 1912: ITAUQ = ITAU
! 1913: ITAUP = ITAUQ + N
! 1914: IWORK = ITAUP + N
! 1915: *
! 1916: * Bidiagonalize R in VT
! 1917: * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
! 1918: * (RWorkspace: need N)
! 1919: *
! 1920: CALL ZGEBRD( N, N, VT, LDVT, S, RWORK( IE ),
! 1921: $ WORK( ITAUQ ), WORK( ITAUP ),
! 1922: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1923: *
! 1924: * Multiply Q in U by left bidiagonalizing vectors
! 1925: * in VT
! 1926: * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
! 1927: * (RWorkspace: 0)
! 1928: *
! 1929: CALL ZUNMBR( 'Q', 'R', 'N', M, N, N, VT, LDVT,
! 1930: $ WORK( ITAUQ ), U, LDU, WORK( IWORK ),
! 1931: $ LWORK-IWORK+1, IERR )
! 1932: *
! 1933: * Generate right bidiagonalizing vectors in VT
! 1934: * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
! 1935: * (RWorkspace: 0)
! 1936: *
! 1937: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
! 1938: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1939: IRWORK = IE + N
! 1940: *
! 1941: * Perform bidiagonal QR iteration, computing left
! 1942: * singular vectors of A in U and computing right
! 1943: * singular vectors of A in VT
! 1944: * (CWorkspace: 0)
! 1945: * (RWorkspace: need BDSPAC)
! 1946: *
! 1947: CALL ZBDSQR( 'U', N, N, M, 0, S, RWORK( IE ), VT,
! 1948: $ LDVT, U, LDU, CDUM, 1,
! 1949: $ RWORK( IRWORK ), INFO )
! 1950: *
! 1951: END IF
! 1952: *
! 1953: END IF
! 1954: *
! 1955: END IF
! 1956: *
! 1957: ELSE
! 1958: *
! 1959: * M .LT. MNTHR
! 1960: *
! 1961: * Path 10 (M at least N, but not much larger)
! 1962: * Reduce to bidiagonal form without QR decomposition
! 1963: *
! 1964: IE = 1
! 1965: ITAUQ = 1
! 1966: ITAUP = ITAUQ + N
! 1967: IWORK = ITAUP + N
! 1968: *
! 1969: * Bidiagonalize A
! 1970: * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
! 1971: * (RWorkspace: need N)
! 1972: *
! 1973: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 1974: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
! 1975: $ IERR )
! 1976: IF( WNTUAS ) THEN
! 1977: *
! 1978: * If left singular vectors desired in U, copy result to U
! 1979: * and generate left bidiagonalizing vectors in U
! 1980: * (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB)
! 1981: * (RWorkspace: 0)
! 1982: *
! 1983: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 1984: IF( WNTUS )
! 1985: $ NCU = N
! 1986: IF( WNTUA )
! 1987: $ NCU = M
! 1988: CALL ZUNGBR( 'Q', M, NCU, N, U, LDU, WORK( ITAUQ ),
! 1989: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 1990: END IF
! 1991: IF( WNTVAS ) THEN
! 1992: *
! 1993: * If right singular vectors desired in VT, copy result to
! 1994: * VT and generate right bidiagonalizing vectors in VT
! 1995: * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
! 1996: * (RWorkspace: 0)
! 1997: *
! 1998: CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
! 1999: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
! 2000: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2001: END IF
! 2002: IF( WNTUO ) THEN
! 2003: *
! 2004: * If left singular vectors desired in A, generate left
! 2005: * bidiagonalizing vectors in A
! 2006: * (CWorkspace: need 3*N, prefer 2*N+N*NB)
! 2007: * (RWorkspace: 0)
! 2008: *
! 2009: CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
! 2010: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2011: END IF
! 2012: IF( WNTVO ) THEN
! 2013: *
! 2014: * If right singular vectors desired in A, generate right
! 2015: * bidiagonalizing vectors in A
! 2016: * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
! 2017: * (RWorkspace: 0)
! 2018: *
! 2019: CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
! 2020: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2021: END IF
! 2022: IRWORK = IE + N
! 2023: IF( WNTUAS .OR. WNTUO )
! 2024: $ NRU = M
! 2025: IF( WNTUN )
! 2026: $ NRU = 0
! 2027: IF( WNTVAS .OR. WNTVO )
! 2028: $ NCVT = N
! 2029: IF( WNTVN )
! 2030: $ NCVT = 0
! 2031: IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
! 2032: *
! 2033: * Perform bidiagonal QR iteration, if desired, computing
! 2034: * left singular vectors in U and computing right singular
! 2035: * vectors in VT
! 2036: * (CWorkspace: 0)
! 2037: * (RWorkspace: need BDSPAC)
! 2038: *
! 2039: CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT,
! 2040: $ LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ),
! 2041: $ INFO )
! 2042: ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
! 2043: *
! 2044: * Perform bidiagonal QR iteration, if desired, computing
! 2045: * left singular vectors in U and computing right singular
! 2046: * vectors in A
! 2047: * (CWorkspace: 0)
! 2048: * (RWorkspace: need BDSPAC)
! 2049: *
! 2050: CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), A,
! 2051: $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
! 2052: $ INFO )
! 2053: ELSE
! 2054: *
! 2055: * Perform bidiagonal QR iteration, if desired, computing
! 2056: * left singular vectors in A and computing right singular
! 2057: * vectors in VT
! 2058: * (CWorkspace: 0)
! 2059: * (RWorkspace: need BDSPAC)
! 2060: *
! 2061: CALL ZBDSQR( 'U', N, NCVT, NRU, 0, S, RWORK( IE ), VT,
! 2062: $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
! 2063: $ INFO )
! 2064: END IF
! 2065: *
! 2066: END IF
! 2067: *
! 2068: ELSE
! 2069: *
! 2070: * A has more columns than rows. If A has sufficiently more
! 2071: * columns than rows, first reduce using the LQ decomposition (if
! 2072: * sufficient workspace available)
! 2073: *
! 2074: IF( N.GE.MNTHR ) THEN
! 2075: *
! 2076: IF( WNTVN ) THEN
! 2077: *
! 2078: * Path 1t(N much larger than M, JOBVT='N')
! 2079: * No right singular vectors to be computed
! 2080: *
! 2081: ITAU = 1
! 2082: IWORK = ITAU + M
! 2083: *
! 2084: * Compute A=L*Q
! 2085: * (CWorkspace: need 2*M, prefer M+M*NB)
! 2086: * (RWorkspace: 0)
! 2087: *
! 2088: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
! 2089: $ LWORK-IWORK+1, IERR )
! 2090: *
! 2091: * Zero out above L
! 2092: *
! 2093: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
! 2094: $ LDA )
! 2095: IE = 1
! 2096: ITAUQ = 1
! 2097: ITAUP = ITAUQ + M
! 2098: IWORK = ITAUP + M
! 2099: *
! 2100: * Bidiagonalize L in A
! 2101: * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
! 2102: * (RWorkspace: need M)
! 2103: *
! 2104: CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 2105: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
! 2106: $ IERR )
! 2107: IF( WNTUO .OR. WNTUAS ) THEN
! 2108: *
! 2109: * If left singular vectors desired, generate Q
! 2110: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 2111: * (RWorkspace: 0)
! 2112: *
! 2113: CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
! 2114: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2115: END IF
! 2116: IRWORK = IE + M
! 2117: NRU = 0
! 2118: IF( WNTUO .OR. WNTUAS )
! 2119: $ NRU = M
! 2120: *
! 2121: * Perform bidiagonal QR iteration, computing left singular
! 2122: * vectors of A in A if desired
! 2123: * (CWorkspace: 0)
! 2124: * (RWorkspace: need BDSPAC)
! 2125: *
! 2126: CALL ZBDSQR( 'U', M, 0, NRU, 0, S, RWORK( IE ), CDUM, 1,
! 2127: $ A, LDA, CDUM, 1, RWORK( IRWORK ), INFO )
! 2128: *
! 2129: * If left singular vectors desired in U, copy them there
! 2130: *
! 2131: IF( WNTUAS )
! 2132: $ CALL ZLACPY( 'F', M, M, A, LDA, U, LDU )
! 2133: *
! 2134: ELSE IF( WNTVO .AND. WNTUN ) THEN
! 2135: *
! 2136: * Path 2t(N much larger than M, JOBU='N', JOBVT='O')
! 2137: * M right singular vectors to be overwritten on A and
! 2138: * no left singular vectors to be computed
! 2139: *
! 2140: IF( LWORK.GE.M*M+3*M ) THEN
! 2141: *
! 2142: * Sufficient workspace for a fast algorithm
! 2143: *
! 2144: IR = 1
! 2145: IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN
! 2146: *
! 2147: * WORK(IU) is LDA by N and WORK(IR) is LDA by M
! 2148: *
! 2149: LDWRKU = LDA
! 2150: CHUNK = N
! 2151: LDWRKR = LDA
! 2152: ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN
! 2153: *
! 2154: * WORK(IU) is LDA by N and WORK(IR) is M by M
! 2155: *
! 2156: LDWRKU = LDA
! 2157: CHUNK = N
! 2158: LDWRKR = M
! 2159: ELSE
! 2160: *
! 2161: * WORK(IU) is M by CHUNK and WORK(IR) is M by M
! 2162: *
! 2163: LDWRKU = M
! 2164: CHUNK = ( LWORK-M*M ) / M
! 2165: LDWRKR = M
! 2166: END IF
! 2167: ITAU = IR + LDWRKR*M
! 2168: IWORK = ITAU + M
! 2169: *
! 2170: * Compute A=L*Q
! 2171: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 2172: * (RWorkspace: 0)
! 2173: *
! 2174: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 2175: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2176: *
! 2177: * Copy L to WORK(IR) and zero out above it
! 2178: *
! 2179: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ), LDWRKR )
! 2180: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 2181: $ WORK( IR+LDWRKR ), LDWRKR )
! 2182: *
! 2183: * Generate Q in A
! 2184: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 2185: * (RWorkspace: 0)
! 2186: *
! 2187: CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
! 2188: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2189: IE = 1
! 2190: ITAUQ = ITAU
! 2191: ITAUP = ITAUQ + M
! 2192: IWORK = ITAUP + M
! 2193: *
! 2194: * Bidiagonalize L in WORK(IR)
! 2195: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
! 2196: * (RWorkspace: need M)
! 2197: *
! 2198: CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S, RWORK( IE ),
! 2199: $ WORK( ITAUQ ), WORK( ITAUP ),
! 2200: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2201: *
! 2202: * Generate right vectors bidiagonalizing L
! 2203: * (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
! 2204: * (RWorkspace: 0)
! 2205: *
! 2206: CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
! 2207: $ WORK( ITAUP ), WORK( IWORK ),
! 2208: $ LWORK-IWORK+1, IERR )
! 2209: IRWORK = IE + M
! 2210: *
! 2211: * Perform bidiagonal QR iteration, computing right
! 2212: * singular vectors of L in WORK(IR)
! 2213: * (CWorkspace: need M*M)
! 2214: * (RWorkspace: need BDSPAC)
! 2215: *
! 2216: CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
! 2217: $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
! 2218: $ RWORK( IRWORK ), INFO )
! 2219: IU = ITAUQ
! 2220: *
! 2221: * Multiply right singular vectors of L in WORK(IR) by Q
! 2222: * in A, storing result in WORK(IU) and copying to A
! 2223: * (CWorkspace: need M*M+M, prefer M*M+M*N)
! 2224: * (RWorkspace: 0)
! 2225: *
! 2226: DO 30 I = 1, N, CHUNK
! 2227: BLK = MIN( N-I+1, CHUNK )
! 2228: CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ),
! 2229: $ LDWRKR, A( 1, I ), LDA, CZERO,
! 2230: $ WORK( IU ), LDWRKU )
! 2231: CALL ZLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
! 2232: $ A( 1, I ), LDA )
! 2233: 30 CONTINUE
! 2234: *
! 2235: ELSE
! 2236: *
! 2237: * Insufficient workspace for a fast algorithm
! 2238: *
! 2239: IE = 1
! 2240: ITAUQ = 1
! 2241: ITAUP = ITAUQ + M
! 2242: IWORK = ITAUP + M
! 2243: *
! 2244: * Bidiagonalize A
! 2245: * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
! 2246: * (RWorkspace: need M)
! 2247: *
! 2248: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ),
! 2249: $ WORK( ITAUQ ), WORK( ITAUP ),
! 2250: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2251: *
! 2252: * Generate right vectors bidiagonalizing A
! 2253: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 2254: * (RWorkspace: 0)
! 2255: *
! 2256: CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
! 2257: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2258: IRWORK = IE + M
! 2259: *
! 2260: * Perform bidiagonal QR iteration, computing right
! 2261: * singular vectors of A in A
! 2262: * (CWorkspace: 0)
! 2263: * (RWorkspace: need BDSPAC)
! 2264: *
! 2265: CALL ZBDSQR( 'L', M, N, 0, 0, S, RWORK( IE ), A, LDA,
! 2266: $ CDUM, 1, CDUM, 1, RWORK( IRWORK ), INFO )
! 2267: *
! 2268: END IF
! 2269: *
! 2270: ELSE IF( WNTVO .AND. WNTUAS ) THEN
! 2271: *
! 2272: * Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
! 2273: * M right singular vectors to be overwritten on A and
! 2274: * M left singular vectors to be computed in U
! 2275: *
! 2276: IF( LWORK.GE.M*M+3*M ) THEN
! 2277: *
! 2278: * Sufficient workspace for a fast algorithm
! 2279: *
! 2280: IR = 1
! 2281: IF( LWORK.GE.MAX( WRKBL, LDA*N )+LDA*M ) THEN
! 2282: *
! 2283: * WORK(IU) is LDA by N and WORK(IR) is LDA by M
! 2284: *
! 2285: LDWRKU = LDA
! 2286: CHUNK = N
! 2287: LDWRKR = LDA
! 2288: ELSE IF( LWORK.GE.MAX( WRKBL, LDA*N )+M*M ) THEN
! 2289: *
! 2290: * WORK(IU) is LDA by N and WORK(IR) is M by M
! 2291: *
! 2292: LDWRKU = LDA
! 2293: CHUNK = N
! 2294: LDWRKR = M
! 2295: ELSE
! 2296: *
! 2297: * WORK(IU) is M by CHUNK and WORK(IR) is M by M
! 2298: *
! 2299: LDWRKU = M
! 2300: CHUNK = ( LWORK-M*M ) / M
! 2301: LDWRKR = M
! 2302: END IF
! 2303: ITAU = IR + LDWRKR*M
! 2304: IWORK = ITAU + M
! 2305: *
! 2306: * Compute A=L*Q
! 2307: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 2308: * (RWorkspace: 0)
! 2309: *
! 2310: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 2311: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2312: *
! 2313: * Copy L to U, zeroing about above it
! 2314: *
! 2315: CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
! 2316: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ),
! 2317: $ LDU )
! 2318: *
! 2319: * Generate Q in A
! 2320: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 2321: * (RWorkspace: 0)
! 2322: *
! 2323: CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
! 2324: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2325: IE = 1
! 2326: ITAUQ = ITAU
! 2327: ITAUP = ITAUQ + M
! 2328: IWORK = ITAUP + M
! 2329: *
! 2330: * Bidiagonalize L in U, copying result to WORK(IR)
! 2331: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
! 2332: * (RWorkspace: need M)
! 2333: *
! 2334: CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
! 2335: $ WORK( ITAUQ ), WORK( ITAUP ),
! 2336: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2337: CALL ZLACPY( 'U', M, M, U, LDU, WORK( IR ), LDWRKR )
! 2338: *
! 2339: * Generate right vectors bidiagonalizing L in WORK(IR)
! 2340: * (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
! 2341: * (RWorkspace: 0)
! 2342: *
! 2343: CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
! 2344: $ WORK( ITAUP ), WORK( IWORK ),
! 2345: $ LWORK-IWORK+1, IERR )
! 2346: *
! 2347: * Generate left vectors bidiagonalizing L in U
! 2348: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
! 2349: * (RWorkspace: 0)
! 2350: *
! 2351: CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
! 2352: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2353: IRWORK = IE + M
! 2354: *
! 2355: * Perform bidiagonal QR iteration, computing left
! 2356: * singular vectors of L in U, and computing right
! 2357: * singular vectors of L in WORK(IR)
! 2358: * (CWorkspace: need M*M)
! 2359: * (RWorkspace: need BDSPAC)
! 2360: *
! 2361: CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
! 2362: $ WORK( IR ), LDWRKR, U, LDU, CDUM, 1,
! 2363: $ RWORK( IRWORK ), INFO )
! 2364: IU = ITAUQ
! 2365: *
! 2366: * Multiply right singular vectors of L in WORK(IR) by Q
! 2367: * in A, storing result in WORK(IU) and copying to A
! 2368: * (CWorkspace: need M*M+M, prefer M*M+M*N))
! 2369: * (RWorkspace: 0)
! 2370: *
! 2371: DO 40 I = 1, N, CHUNK
! 2372: BLK = MIN( N-I+1, CHUNK )
! 2373: CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IR ),
! 2374: $ LDWRKR, A( 1, I ), LDA, CZERO,
! 2375: $ WORK( IU ), LDWRKU )
! 2376: CALL ZLACPY( 'F', M, BLK, WORK( IU ), LDWRKU,
! 2377: $ A( 1, I ), LDA )
! 2378: 40 CONTINUE
! 2379: *
! 2380: ELSE
! 2381: *
! 2382: * Insufficient workspace for a fast algorithm
! 2383: *
! 2384: ITAU = 1
! 2385: IWORK = ITAU + M
! 2386: *
! 2387: * Compute A=L*Q
! 2388: * (CWorkspace: need 2*M, prefer M+M*NB)
! 2389: * (RWorkspace: 0)
! 2390: *
! 2391: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 2392: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2393: *
! 2394: * Copy L to U, zeroing out above it
! 2395: *
! 2396: CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
! 2397: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, U( 1, 2 ),
! 2398: $ LDU )
! 2399: *
! 2400: * Generate Q in A
! 2401: * (CWorkspace: need 2*M, prefer M+M*NB)
! 2402: * (RWorkspace: 0)
! 2403: *
! 2404: CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
! 2405: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2406: IE = 1
! 2407: ITAUQ = ITAU
! 2408: ITAUP = ITAUQ + M
! 2409: IWORK = ITAUP + M
! 2410: *
! 2411: * Bidiagonalize L in U
! 2412: * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
! 2413: * (RWorkspace: need M)
! 2414: *
! 2415: CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
! 2416: $ WORK( ITAUQ ), WORK( ITAUP ),
! 2417: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2418: *
! 2419: * Multiply right vectors bidiagonalizing L by Q in A
! 2420: * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
! 2421: * (RWorkspace: 0)
! 2422: *
! 2423: CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
! 2424: $ WORK( ITAUP ), A, LDA, WORK( IWORK ),
! 2425: $ LWORK-IWORK+1, IERR )
! 2426: *
! 2427: * Generate left vectors bidiagonalizing L in U
! 2428: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 2429: * (RWorkspace: 0)
! 2430: *
! 2431: CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
! 2432: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2433: IRWORK = IE + M
! 2434: *
! 2435: * Perform bidiagonal QR iteration, computing left
! 2436: * singular vectors of A in U and computing right
! 2437: * singular vectors of A in A
! 2438: * (CWorkspace: 0)
! 2439: * (RWorkspace: need BDSPAC)
! 2440: *
! 2441: CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), A, LDA,
! 2442: $ U, LDU, CDUM, 1, RWORK( IRWORK ), INFO )
! 2443: *
! 2444: END IF
! 2445: *
! 2446: ELSE IF( WNTVS ) THEN
! 2447: *
! 2448: IF( WNTUN ) THEN
! 2449: *
! 2450: * Path 4t(N much larger than M, JOBU='N', JOBVT='S')
! 2451: * M right singular vectors to be computed in VT and
! 2452: * no left singular vectors to be computed
! 2453: *
! 2454: IF( LWORK.GE.M*M+3*M ) THEN
! 2455: *
! 2456: * Sufficient workspace for a fast algorithm
! 2457: *
! 2458: IR = 1
! 2459: IF( LWORK.GE.WRKBL+LDA*M ) THEN
! 2460: *
! 2461: * WORK(IR) is LDA by M
! 2462: *
! 2463: LDWRKR = LDA
! 2464: ELSE
! 2465: *
! 2466: * WORK(IR) is M by M
! 2467: *
! 2468: LDWRKR = M
! 2469: END IF
! 2470: ITAU = IR + LDWRKR*M
! 2471: IWORK = ITAU + M
! 2472: *
! 2473: * Compute A=L*Q
! 2474: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 2475: * (RWorkspace: 0)
! 2476: *
! 2477: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 2478: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2479: *
! 2480: * Copy L to WORK(IR), zeroing out above it
! 2481: *
! 2482: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ),
! 2483: $ LDWRKR )
! 2484: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 2485: $ WORK( IR+LDWRKR ), LDWRKR )
! 2486: *
! 2487: * Generate Q in A
! 2488: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 2489: * (RWorkspace: 0)
! 2490: *
! 2491: CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
! 2492: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2493: IE = 1
! 2494: ITAUQ = ITAU
! 2495: ITAUP = ITAUQ + M
! 2496: IWORK = ITAUP + M
! 2497: *
! 2498: * Bidiagonalize L in WORK(IR)
! 2499: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
! 2500: * (RWorkspace: need M)
! 2501: *
! 2502: CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S,
! 2503: $ RWORK( IE ), WORK( ITAUQ ),
! 2504: $ WORK( ITAUP ), WORK( IWORK ),
! 2505: $ LWORK-IWORK+1, IERR )
! 2506: *
! 2507: * Generate right vectors bidiagonalizing L in
! 2508: * WORK(IR)
! 2509: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
! 2510: * (RWorkspace: 0)
! 2511: *
! 2512: CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
! 2513: $ WORK( ITAUP ), WORK( IWORK ),
! 2514: $ LWORK-IWORK+1, IERR )
! 2515: IRWORK = IE + M
! 2516: *
! 2517: * Perform bidiagonal QR iteration, computing right
! 2518: * singular vectors of L in WORK(IR)
! 2519: * (CWorkspace: need M*M)
! 2520: * (RWorkspace: need BDSPAC)
! 2521: *
! 2522: CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
! 2523: $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
! 2524: $ RWORK( IRWORK ), INFO )
! 2525: *
! 2526: * Multiply right singular vectors of L in WORK(IR) by
! 2527: * Q in A, storing result in VT
! 2528: * (CWorkspace: need M*M)
! 2529: * (RWorkspace: 0)
! 2530: *
! 2531: CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ),
! 2532: $ LDWRKR, A, LDA, CZERO, VT, LDVT )
! 2533: *
! 2534: ELSE
! 2535: *
! 2536: * Insufficient workspace for a fast algorithm
! 2537: *
! 2538: ITAU = 1
! 2539: IWORK = ITAU + M
! 2540: *
! 2541: * Compute A=L*Q
! 2542: * (CWorkspace: need 2*M, prefer M+M*NB)
! 2543: * (RWorkspace: 0)
! 2544: *
! 2545: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 2546: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2547: *
! 2548: * Copy result to VT
! 2549: *
! 2550: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 2551: *
! 2552: * Generate Q in VT
! 2553: * (CWorkspace: need 2*M, prefer M+M*NB)
! 2554: * (RWorkspace: 0)
! 2555: *
! 2556: CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
! 2557: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2558: IE = 1
! 2559: ITAUQ = ITAU
! 2560: ITAUP = ITAUQ + M
! 2561: IWORK = ITAUP + M
! 2562: *
! 2563: * Zero out above L in A
! 2564: *
! 2565: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 2566: $ A( 1, 2 ), LDA )
! 2567: *
! 2568: * Bidiagonalize L in A
! 2569: * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
! 2570: * (RWorkspace: need M)
! 2571: *
! 2572: CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
! 2573: $ WORK( ITAUQ ), WORK( ITAUP ),
! 2574: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2575: *
! 2576: * Multiply right vectors bidiagonalizing L by Q in VT
! 2577: * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
! 2578: * (RWorkspace: 0)
! 2579: *
! 2580: CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
! 2581: $ WORK( ITAUP ), VT, LDVT,
! 2582: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2583: IRWORK = IE + M
! 2584: *
! 2585: * Perform bidiagonal QR iteration, computing right
! 2586: * singular vectors of A in VT
! 2587: * (CWorkspace: 0)
! 2588: * (RWorkspace: need BDSPAC)
! 2589: *
! 2590: CALL ZBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT,
! 2591: $ LDVT, CDUM, 1, CDUM, 1,
! 2592: $ RWORK( IRWORK ), INFO )
! 2593: *
! 2594: END IF
! 2595: *
! 2596: ELSE IF( WNTUO ) THEN
! 2597: *
! 2598: * Path 5t(N much larger than M, JOBU='O', JOBVT='S')
! 2599: * M right singular vectors to be computed in VT and
! 2600: * M left singular vectors to be overwritten on A
! 2601: *
! 2602: IF( LWORK.GE.2*M*M+3*M ) THEN
! 2603: *
! 2604: * Sufficient workspace for a fast algorithm
! 2605: *
! 2606: IU = 1
! 2607: IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
! 2608: *
! 2609: * WORK(IU) is LDA by M and WORK(IR) is LDA by M
! 2610: *
! 2611: LDWRKU = LDA
! 2612: IR = IU + LDWRKU*M
! 2613: LDWRKR = LDA
! 2614: ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
! 2615: *
! 2616: * WORK(IU) is LDA by M and WORK(IR) is M by M
! 2617: *
! 2618: LDWRKU = LDA
! 2619: IR = IU + LDWRKU*M
! 2620: LDWRKR = M
! 2621: ELSE
! 2622: *
! 2623: * WORK(IU) is M by M and WORK(IR) is M by M
! 2624: *
! 2625: LDWRKU = M
! 2626: IR = IU + LDWRKU*M
! 2627: LDWRKR = M
! 2628: END IF
! 2629: ITAU = IR + LDWRKR*M
! 2630: IWORK = ITAU + M
! 2631: *
! 2632: * Compute A=L*Q
! 2633: * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
! 2634: * (RWorkspace: 0)
! 2635: *
! 2636: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 2637: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2638: *
! 2639: * Copy L to WORK(IU), zeroing out below it
! 2640: *
! 2641: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
! 2642: $ LDWRKU )
! 2643: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 2644: $ WORK( IU+LDWRKU ), LDWRKU )
! 2645: *
! 2646: * Generate Q in A
! 2647: * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
! 2648: * (RWorkspace: 0)
! 2649: *
! 2650: CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
! 2651: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2652: IE = 1
! 2653: ITAUQ = ITAU
! 2654: ITAUP = ITAUQ + M
! 2655: IWORK = ITAUP + M
! 2656: *
! 2657: * Bidiagonalize L in WORK(IU), copying result to
! 2658: * WORK(IR)
! 2659: * (CWorkspace: need 2*M*M+3*M,
! 2660: * prefer 2*M*M+2*M+2*M*NB)
! 2661: * (RWorkspace: need M)
! 2662: *
! 2663: CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
! 2664: $ RWORK( IE ), WORK( ITAUQ ),
! 2665: $ WORK( ITAUP ), WORK( IWORK ),
! 2666: $ LWORK-IWORK+1, IERR )
! 2667: CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU,
! 2668: $ WORK( IR ), LDWRKR )
! 2669: *
! 2670: * Generate right bidiagonalizing vectors in WORK(IU)
! 2671: * (CWorkspace: need 2*M*M+3*M-1,
! 2672: * prefer 2*M*M+2*M+(M-1)*NB)
! 2673: * (RWorkspace: 0)
! 2674: *
! 2675: CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
! 2676: $ WORK( ITAUP ), WORK( IWORK ),
! 2677: $ LWORK-IWORK+1, IERR )
! 2678: *
! 2679: * Generate left bidiagonalizing vectors in WORK(IR)
! 2680: * (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
! 2681: * (RWorkspace: 0)
! 2682: *
! 2683: CALL ZUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
! 2684: $ WORK( ITAUQ ), WORK( IWORK ),
! 2685: $ LWORK-IWORK+1, IERR )
! 2686: IRWORK = IE + M
! 2687: *
! 2688: * Perform bidiagonal QR iteration, computing left
! 2689: * singular vectors of L in WORK(IR) and computing
! 2690: * right singular vectors of L in WORK(IU)
! 2691: * (CWorkspace: need 2*M*M)
! 2692: * (RWorkspace: need BDSPAC)
! 2693: *
! 2694: CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
! 2695: $ WORK( IU ), LDWRKU, WORK( IR ),
! 2696: $ LDWRKR, CDUM, 1, RWORK( IRWORK ),
! 2697: $ INFO )
! 2698: *
! 2699: * Multiply right singular vectors of L in WORK(IU) by
! 2700: * Q in A, storing result in VT
! 2701: * (CWorkspace: need M*M)
! 2702: * (RWorkspace: 0)
! 2703: *
! 2704: CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
! 2705: $ LDWRKU, A, LDA, CZERO, VT, LDVT )
! 2706: *
! 2707: * Copy left singular vectors of L to A
! 2708: * (CWorkspace: need M*M)
! 2709: * (RWorkspace: 0)
! 2710: *
! 2711: CALL ZLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
! 2712: $ LDA )
! 2713: *
! 2714: ELSE
! 2715: *
! 2716: * Insufficient workspace for a fast algorithm
! 2717: *
! 2718: ITAU = 1
! 2719: IWORK = ITAU + M
! 2720: *
! 2721: * Compute A=L*Q, copying result to VT
! 2722: * (CWorkspace: need 2*M, prefer M+M*NB)
! 2723: * (RWorkspace: 0)
! 2724: *
! 2725: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 2726: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2727: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 2728: *
! 2729: * Generate Q in VT
! 2730: * (CWorkspace: need 2*M, prefer M+M*NB)
! 2731: * (RWorkspace: 0)
! 2732: *
! 2733: CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
! 2734: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2735: IE = 1
! 2736: ITAUQ = ITAU
! 2737: ITAUP = ITAUQ + M
! 2738: IWORK = ITAUP + M
! 2739: *
! 2740: * Zero out above L in A
! 2741: *
! 2742: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 2743: $ A( 1, 2 ), LDA )
! 2744: *
! 2745: * Bidiagonalize L in A
! 2746: * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
! 2747: * (RWorkspace: need M)
! 2748: *
! 2749: CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
! 2750: $ WORK( ITAUQ ), WORK( ITAUP ),
! 2751: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2752: *
! 2753: * Multiply right vectors bidiagonalizing L by Q in VT
! 2754: * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
! 2755: * (RWorkspace: 0)
! 2756: *
! 2757: CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
! 2758: $ WORK( ITAUP ), VT, LDVT,
! 2759: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2760: *
! 2761: * Generate left bidiagonalizing vectors of L in A
! 2762: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 2763: * (RWorkspace: 0)
! 2764: *
! 2765: CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
! 2766: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2767: IRWORK = IE + M
! 2768: *
! 2769: * Perform bidiagonal QR iteration, computing left
! 2770: * singular vectors of A in A and computing right
! 2771: * singular vectors of A in VT
! 2772: * (CWorkspace: 0)
! 2773: * (RWorkspace: need BDSPAC)
! 2774: *
! 2775: CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
! 2776: $ LDVT, A, LDA, CDUM, 1,
! 2777: $ RWORK( IRWORK ), INFO )
! 2778: *
! 2779: END IF
! 2780: *
! 2781: ELSE IF( WNTUAS ) THEN
! 2782: *
! 2783: * Path 6t(N much larger than M, JOBU='S' or 'A',
! 2784: * JOBVT='S')
! 2785: * M right singular vectors to be computed in VT and
! 2786: * M left singular vectors to be computed in U
! 2787: *
! 2788: IF( LWORK.GE.M*M+3*M ) THEN
! 2789: *
! 2790: * Sufficient workspace for a fast algorithm
! 2791: *
! 2792: IU = 1
! 2793: IF( LWORK.GE.WRKBL+LDA*M ) THEN
! 2794: *
! 2795: * WORK(IU) is LDA by N
! 2796: *
! 2797: LDWRKU = LDA
! 2798: ELSE
! 2799: *
! 2800: * WORK(IU) is LDA by M
! 2801: *
! 2802: LDWRKU = M
! 2803: END IF
! 2804: ITAU = IU + LDWRKU*M
! 2805: IWORK = ITAU + M
! 2806: *
! 2807: * Compute A=L*Q
! 2808: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 2809: * (RWorkspace: 0)
! 2810: *
! 2811: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 2812: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2813: *
! 2814: * Copy L to WORK(IU), zeroing out above it
! 2815: *
! 2816: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
! 2817: $ LDWRKU )
! 2818: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 2819: $ WORK( IU+LDWRKU ), LDWRKU )
! 2820: *
! 2821: * Generate Q in A
! 2822: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 2823: * (RWorkspace: 0)
! 2824: *
! 2825: CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
! 2826: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2827: IE = 1
! 2828: ITAUQ = ITAU
! 2829: ITAUP = ITAUQ + M
! 2830: IWORK = ITAUP + M
! 2831: *
! 2832: * Bidiagonalize L in WORK(IU), copying result to U
! 2833: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
! 2834: * (RWorkspace: need M)
! 2835: *
! 2836: CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
! 2837: $ RWORK( IE ), WORK( ITAUQ ),
! 2838: $ WORK( ITAUP ), WORK( IWORK ),
! 2839: $ LWORK-IWORK+1, IERR )
! 2840: CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
! 2841: $ LDU )
! 2842: *
! 2843: * Generate right bidiagonalizing vectors in WORK(IU)
! 2844: * (CWorkspace: need M*M+3*M-1,
! 2845: * prefer M*M+2*M+(M-1)*NB)
! 2846: * (RWorkspace: 0)
! 2847: *
! 2848: CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
! 2849: $ WORK( ITAUP ), WORK( IWORK ),
! 2850: $ LWORK-IWORK+1, IERR )
! 2851: *
! 2852: * Generate left bidiagonalizing vectors in U
! 2853: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
! 2854: * (RWorkspace: 0)
! 2855: *
! 2856: CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
! 2857: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2858: IRWORK = IE + M
! 2859: *
! 2860: * Perform bidiagonal QR iteration, computing left
! 2861: * singular vectors of L in U and computing right
! 2862: * singular vectors of L in WORK(IU)
! 2863: * (CWorkspace: need M*M)
! 2864: * (RWorkspace: need BDSPAC)
! 2865: *
! 2866: CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
! 2867: $ WORK( IU ), LDWRKU, U, LDU, CDUM, 1,
! 2868: $ RWORK( IRWORK ), INFO )
! 2869: *
! 2870: * Multiply right singular vectors of L in WORK(IU) by
! 2871: * Q in A, storing result in VT
! 2872: * (CWorkspace: need M*M)
! 2873: * (RWorkspace: 0)
! 2874: *
! 2875: CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
! 2876: $ LDWRKU, A, LDA, CZERO, VT, LDVT )
! 2877: *
! 2878: ELSE
! 2879: *
! 2880: * Insufficient workspace for a fast algorithm
! 2881: *
! 2882: ITAU = 1
! 2883: IWORK = ITAU + M
! 2884: *
! 2885: * Compute A=L*Q, copying result to VT
! 2886: * (CWorkspace: need 2*M, prefer M+M*NB)
! 2887: * (RWorkspace: 0)
! 2888: *
! 2889: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 2890: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2891: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 2892: *
! 2893: * Generate Q in VT
! 2894: * (CWorkspace: need 2*M, prefer M+M*NB)
! 2895: * (RWorkspace: 0)
! 2896: *
! 2897: CALL ZUNGLQ( M, N, M, VT, LDVT, WORK( ITAU ),
! 2898: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2899: *
! 2900: * Copy L to U, zeroing out above it
! 2901: *
! 2902: CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
! 2903: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 2904: $ U( 1, 2 ), LDU )
! 2905: IE = 1
! 2906: ITAUQ = ITAU
! 2907: ITAUP = ITAUQ + M
! 2908: IWORK = ITAUP + M
! 2909: *
! 2910: * Bidiagonalize L in U
! 2911: * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
! 2912: * (RWorkspace: need M)
! 2913: *
! 2914: CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
! 2915: $ WORK( ITAUQ ), WORK( ITAUP ),
! 2916: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2917: *
! 2918: * Multiply right bidiagonalizing vectors in U by Q
! 2919: * in VT
! 2920: * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
! 2921: * (RWorkspace: 0)
! 2922: *
! 2923: CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
! 2924: $ WORK( ITAUP ), VT, LDVT,
! 2925: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2926: *
! 2927: * Generate left bidiagonalizing vectors in U
! 2928: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 2929: * (RWorkspace: 0)
! 2930: *
! 2931: CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
! 2932: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2933: IRWORK = IE + M
! 2934: *
! 2935: * Perform bidiagonal QR iteration, computing left
! 2936: * singular vectors of A in U and computing right
! 2937: * singular vectors of A in VT
! 2938: * (CWorkspace: 0)
! 2939: * (RWorkspace: need BDSPAC)
! 2940: *
! 2941: CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
! 2942: $ LDVT, U, LDU, CDUM, 1,
! 2943: $ RWORK( IRWORK ), INFO )
! 2944: *
! 2945: END IF
! 2946: *
! 2947: END IF
! 2948: *
! 2949: ELSE IF( WNTVA ) THEN
! 2950: *
! 2951: IF( WNTUN ) THEN
! 2952: *
! 2953: * Path 7t(N much larger than M, JOBU='N', JOBVT='A')
! 2954: * N right singular vectors to be computed in VT and
! 2955: * no left singular vectors to be computed
! 2956: *
! 2957: IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN
! 2958: *
! 2959: * Sufficient workspace for a fast algorithm
! 2960: *
! 2961: IR = 1
! 2962: IF( LWORK.GE.WRKBL+LDA*M ) THEN
! 2963: *
! 2964: * WORK(IR) is LDA by M
! 2965: *
! 2966: LDWRKR = LDA
! 2967: ELSE
! 2968: *
! 2969: * WORK(IR) is M by M
! 2970: *
! 2971: LDWRKR = M
! 2972: END IF
! 2973: ITAU = IR + LDWRKR*M
! 2974: IWORK = ITAU + M
! 2975: *
! 2976: * Compute A=L*Q, copying result to VT
! 2977: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 2978: * (RWorkspace: 0)
! 2979: *
! 2980: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 2981: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2982: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 2983: *
! 2984: * Copy L to WORK(IR), zeroing out above it
! 2985: *
! 2986: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IR ),
! 2987: $ LDWRKR )
! 2988: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 2989: $ WORK( IR+LDWRKR ), LDWRKR )
! 2990: *
! 2991: * Generate Q in VT
! 2992: * (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
! 2993: * (RWorkspace: 0)
! 2994: *
! 2995: CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
! 2996: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 2997: IE = 1
! 2998: ITAUQ = ITAU
! 2999: ITAUP = ITAUQ + M
! 3000: IWORK = ITAUP + M
! 3001: *
! 3002: * Bidiagonalize L in WORK(IR)
! 3003: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
! 3004: * (RWorkspace: need M)
! 3005: *
! 3006: CALL ZGEBRD( M, M, WORK( IR ), LDWRKR, S,
! 3007: $ RWORK( IE ), WORK( ITAUQ ),
! 3008: $ WORK( ITAUP ), WORK( IWORK ),
! 3009: $ LWORK-IWORK+1, IERR )
! 3010: *
! 3011: * Generate right bidiagonalizing vectors in WORK(IR)
! 3012: * (CWorkspace: need M*M+3*M-1,
! 3013: * prefer M*M+2*M+(M-1)*NB)
! 3014: * (RWorkspace: 0)
! 3015: *
! 3016: CALL ZUNGBR( 'P', M, M, M, WORK( IR ), LDWRKR,
! 3017: $ WORK( ITAUP ), WORK( IWORK ),
! 3018: $ LWORK-IWORK+1, IERR )
! 3019: IRWORK = IE + M
! 3020: *
! 3021: * Perform bidiagonal QR iteration, computing right
! 3022: * singular vectors of L in WORK(IR)
! 3023: * (CWorkspace: need M*M)
! 3024: * (RWorkspace: need BDSPAC)
! 3025: *
! 3026: CALL ZBDSQR( 'U', M, M, 0, 0, S, RWORK( IE ),
! 3027: $ WORK( IR ), LDWRKR, CDUM, 1, CDUM, 1,
! 3028: $ RWORK( IRWORK ), INFO )
! 3029: *
! 3030: * Multiply right singular vectors of L in WORK(IR) by
! 3031: * Q in VT, storing result in A
! 3032: * (CWorkspace: need M*M)
! 3033: * (RWorkspace: 0)
! 3034: *
! 3035: CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IR ),
! 3036: $ LDWRKR, VT, LDVT, CZERO, A, LDA )
! 3037: *
! 3038: * Copy right singular vectors of A from A to VT
! 3039: *
! 3040: CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
! 3041: *
! 3042: ELSE
! 3043: *
! 3044: * Insufficient workspace for a fast algorithm
! 3045: *
! 3046: ITAU = 1
! 3047: IWORK = ITAU + M
! 3048: *
! 3049: * Compute A=L*Q, copying result to VT
! 3050: * (CWorkspace: need 2*M, prefer M+M*NB)
! 3051: * (RWorkspace: 0)
! 3052: *
! 3053: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 3054: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3055: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 3056: *
! 3057: * Generate Q in VT
! 3058: * (CWorkspace: need M+N, prefer M+N*NB)
! 3059: * (RWorkspace: 0)
! 3060: *
! 3061: CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
! 3062: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3063: IE = 1
! 3064: ITAUQ = ITAU
! 3065: ITAUP = ITAUQ + M
! 3066: IWORK = ITAUP + M
! 3067: *
! 3068: * Zero out above L in A
! 3069: *
! 3070: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 3071: $ A( 1, 2 ), LDA )
! 3072: *
! 3073: * Bidiagonalize L in A
! 3074: * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
! 3075: * (RWorkspace: need M)
! 3076: *
! 3077: CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
! 3078: $ WORK( ITAUQ ), WORK( ITAUP ),
! 3079: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3080: *
! 3081: * Multiply right bidiagonalizing vectors in A by Q
! 3082: * in VT
! 3083: * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
! 3084: * (RWorkspace: 0)
! 3085: *
! 3086: CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
! 3087: $ WORK( ITAUP ), VT, LDVT,
! 3088: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3089: IRWORK = IE + M
! 3090: *
! 3091: * Perform bidiagonal QR iteration, computing right
! 3092: * singular vectors of A in VT
! 3093: * (CWorkspace: 0)
! 3094: * (RWorkspace: need BDSPAC)
! 3095: *
! 3096: CALL ZBDSQR( 'U', M, N, 0, 0, S, RWORK( IE ), VT,
! 3097: $ LDVT, CDUM, 1, CDUM, 1,
! 3098: $ RWORK( IRWORK ), INFO )
! 3099: *
! 3100: END IF
! 3101: *
! 3102: ELSE IF( WNTUO ) THEN
! 3103: *
! 3104: * Path 8t(N much larger than M, JOBU='O', JOBVT='A')
! 3105: * N right singular vectors to be computed in VT and
! 3106: * M left singular vectors to be overwritten on A
! 3107: *
! 3108: IF( LWORK.GE.2*M*M+MAX( N+M, 3*M ) ) THEN
! 3109: *
! 3110: * Sufficient workspace for a fast algorithm
! 3111: *
! 3112: IU = 1
! 3113: IF( LWORK.GE.WRKBL+2*LDA*M ) THEN
! 3114: *
! 3115: * WORK(IU) is LDA by M and WORK(IR) is LDA by M
! 3116: *
! 3117: LDWRKU = LDA
! 3118: IR = IU + LDWRKU*M
! 3119: LDWRKR = LDA
! 3120: ELSE IF( LWORK.GE.WRKBL+( LDA+M )*M ) THEN
! 3121: *
! 3122: * WORK(IU) is LDA by M and WORK(IR) is M by M
! 3123: *
! 3124: LDWRKU = LDA
! 3125: IR = IU + LDWRKU*M
! 3126: LDWRKR = M
! 3127: ELSE
! 3128: *
! 3129: * WORK(IU) is M by M and WORK(IR) is M by M
! 3130: *
! 3131: LDWRKU = M
! 3132: IR = IU + LDWRKU*M
! 3133: LDWRKR = M
! 3134: END IF
! 3135: ITAU = IR + LDWRKR*M
! 3136: IWORK = ITAU + M
! 3137: *
! 3138: * Compute A=L*Q, copying result to VT
! 3139: * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
! 3140: * (RWorkspace: 0)
! 3141: *
! 3142: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 3143: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3144: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 3145: *
! 3146: * Generate Q in VT
! 3147: * (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
! 3148: * (RWorkspace: 0)
! 3149: *
! 3150: CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
! 3151: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3152: *
! 3153: * Copy L to WORK(IU), zeroing out above it
! 3154: *
! 3155: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
! 3156: $ LDWRKU )
! 3157: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 3158: $ WORK( IU+LDWRKU ), LDWRKU )
! 3159: IE = 1
! 3160: ITAUQ = ITAU
! 3161: ITAUP = ITAUQ + M
! 3162: IWORK = ITAUP + M
! 3163: *
! 3164: * Bidiagonalize L in WORK(IU), copying result to
! 3165: * WORK(IR)
! 3166: * (CWorkspace: need 2*M*M+3*M,
! 3167: * prefer 2*M*M+2*M+2*M*NB)
! 3168: * (RWorkspace: need M)
! 3169: *
! 3170: CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
! 3171: $ RWORK( IE ), WORK( ITAUQ ),
! 3172: $ WORK( ITAUP ), WORK( IWORK ),
! 3173: $ LWORK-IWORK+1, IERR )
! 3174: CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU,
! 3175: $ WORK( IR ), LDWRKR )
! 3176: *
! 3177: * Generate right bidiagonalizing vectors in WORK(IU)
! 3178: * (CWorkspace: need 2*M*M+3*M-1,
! 3179: * prefer 2*M*M+2*M+(M-1)*NB)
! 3180: * (RWorkspace: 0)
! 3181: *
! 3182: CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
! 3183: $ WORK( ITAUP ), WORK( IWORK ),
! 3184: $ LWORK-IWORK+1, IERR )
! 3185: *
! 3186: * Generate left bidiagonalizing vectors in WORK(IR)
! 3187: * (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
! 3188: * (RWorkspace: 0)
! 3189: *
! 3190: CALL ZUNGBR( 'Q', M, M, M, WORK( IR ), LDWRKR,
! 3191: $ WORK( ITAUQ ), WORK( IWORK ),
! 3192: $ LWORK-IWORK+1, IERR )
! 3193: IRWORK = IE + M
! 3194: *
! 3195: * Perform bidiagonal QR iteration, computing left
! 3196: * singular vectors of L in WORK(IR) and computing
! 3197: * right singular vectors of L in WORK(IU)
! 3198: * (CWorkspace: need 2*M*M)
! 3199: * (RWorkspace: need BDSPAC)
! 3200: *
! 3201: CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
! 3202: $ WORK( IU ), LDWRKU, WORK( IR ),
! 3203: $ LDWRKR, CDUM, 1, RWORK( IRWORK ),
! 3204: $ INFO )
! 3205: *
! 3206: * Multiply right singular vectors of L in WORK(IU) by
! 3207: * Q in VT, storing result in A
! 3208: * (CWorkspace: need M*M)
! 3209: * (RWorkspace: 0)
! 3210: *
! 3211: CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
! 3212: $ LDWRKU, VT, LDVT, CZERO, A, LDA )
! 3213: *
! 3214: * Copy right singular vectors of A from A to VT
! 3215: *
! 3216: CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
! 3217: *
! 3218: * Copy left singular vectors of A from WORK(IR) to A
! 3219: *
! 3220: CALL ZLACPY( 'F', M, M, WORK( IR ), LDWRKR, A,
! 3221: $ LDA )
! 3222: *
! 3223: ELSE
! 3224: *
! 3225: * Insufficient workspace for a fast algorithm
! 3226: *
! 3227: ITAU = 1
! 3228: IWORK = ITAU + M
! 3229: *
! 3230: * Compute A=L*Q, copying result to VT
! 3231: * (CWorkspace: need 2*M, prefer M+M*NB)
! 3232: * (RWorkspace: 0)
! 3233: *
! 3234: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 3235: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3236: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 3237: *
! 3238: * Generate Q in VT
! 3239: * (CWorkspace: need M+N, prefer M+N*NB)
! 3240: * (RWorkspace: 0)
! 3241: *
! 3242: CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
! 3243: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3244: IE = 1
! 3245: ITAUQ = ITAU
! 3246: ITAUP = ITAUQ + M
! 3247: IWORK = ITAUP + M
! 3248: *
! 3249: * Zero out above L in A
! 3250: *
! 3251: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 3252: $ A( 1, 2 ), LDA )
! 3253: *
! 3254: * Bidiagonalize L in A
! 3255: * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
! 3256: * (RWorkspace: need M)
! 3257: *
! 3258: CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ),
! 3259: $ WORK( ITAUQ ), WORK( ITAUP ),
! 3260: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3261: *
! 3262: * Multiply right bidiagonalizing vectors in A by Q
! 3263: * in VT
! 3264: * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
! 3265: * (RWorkspace: 0)
! 3266: *
! 3267: CALL ZUNMBR( 'P', 'L', 'C', M, N, M, A, LDA,
! 3268: $ WORK( ITAUP ), VT, LDVT,
! 3269: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3270: *
! 3271: * Generate left bidiagonalizing vectors in A
! 3272: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 3273: * (RWorkspace: 0)
! 3274: *
! 3275: CALL ZUNGBR( 'Q', M, M, M, A, LDA, WORK( ITAUQ ),
! 3276: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3277: IRWORK = IE + M
! 3278: *
! 3279: * Perform bidiagonal QR iteration, computing left
! 3280: * singular vectors of A in A and computing right
! 3281: * singular vectors of A in VT
! 3282: * (CWorkspace: 0)
! 3283: * (RWorkspace: need BDSPAC)
! 3284: *
! 3285: CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
! 3286: $ LDVT, A, LDA, CDUM, 1,
! 3287: $ RWORK( IRWORK ), INFO )
! 3288: *
! 3289: END IF
! 3290: *
! 3291: ELSE IF( WNTUAS ) THEN
! 3292: *
! 3293: * Path 9t(N much larger than M, JOBU='S' or 'A',
! 3294: * JOBVT='A')
! 3295: * N right singular vectors to be computed in VT and
! 3296: * M left singular vectors to be computed in U
! 3297: *
! 3298: IF( LWORK.GE.M*M+MAX( N+M, 3*M ) ) THEN
! 3299: *
! 3300: * Sufficient workspace for a fast algorithm
! 3301: *
! 3302: IU = 1
! 3303: IF( LWORK.GE.WRKBL+LDA*M ) THEN
! 3304: *
! 3305: * WORK(IU) is LDA by M
! 3306: *
! 3307: LDWRKU = LDA
! 3308: ELSE
! 3309: *
! 3310: * WORK(IU) is M by M
! 3311: *
! 3312: LDWRKU = M
! 3313: END IF
! 3314: ITAU = IU + LDWRKU*M
! 3315: IWORK = ITAU + M
! 3316: *
! 3317: * Compute A=L*Q, copying result to VT
! 3318: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 3319: * (RWorkspace: 0)
! 3320: *
! 3321: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 3322: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3323: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 3324: *
! 3325: * Generate Q in VT
! 3326: * (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
! 3327: * (RWorkspace: 0)
! 3328: *
! 3329: CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
! 3330: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3331: *
! 3332: * Copy L to WORK(IU), zeroing out above it
! 3333: *
! 3334: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IU ),
! 3335: $ LDWRKU )
! 3336: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 3337: $ WORK( IU+LDWRKU ), LDWRKU )
! 3338: IE = 1
! 3339: ITAUQ = ITAU
! 3340: ITAUP = ITAUQ + M
! 3341: IWORK = ITAUP + M
! 3342: *
! 3343: * Bidiagonalize L in WORK(IU), copying result to U
! 3344: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
! 3345: * (RWorkspace: need M)
! 3346: *
! 3347: CALL ZGEBRD( M, M, WORK( IU ), LDWRKU, S,
! 3348: $ RWORK( IE ), WORK( ITAUQ ),
! 3349: $ WORK( ITAUP ), WORK( IWORK ),
! 3350: $ LWORK-IWORK+1, IERR )
! 3351: CALL ZLACPY( 'L', M, M, WORK( IU ), LDWRKU, U,
! 3352: $ LDU )
! 3353: *
! 3354: * Generate right bidiagonalizing vectors in WORK(IU)
! 3355: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
! 3356: * (RWorkspace: 0)
! 3357: *
! 3358: CALL ZUNGBR( 'P', M, M, M, WORK( IU ), LDWRKU,
! 3359: $ WORK( ITAUP ), WORK( IWORK ),
! 3360: $ LWORK-IWORK+1, IERR )
! 3361: *
! 3362: * Generate left bidiagonalizing vectors in U
! 3363: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
! 3364: * (RWorkspace: 0)
! 3365: *
! 3366: CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
! 3367: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3368: IRWORK = IE + M
! 3369: *
! 3370: * Perform bidiagonal QR iteration, computing left
! 3371: * singular vectors of L in U and computing right
! 3372: * singular vectors of L in WORK(IU)
! 3373: * (CWorkspace: need M*M)
! 3374: * (RWorkspace: need BDSPAC)
! 3375: *
! 3376: CALL ZBDSQR( 'U', M, M, M, 0, S, RWORK( IE ),
! 3377: $ WORK( IU ), LDWRKU, U, LDU, CDUM, 1,
! 3378: $ RWORK( IRWORK ), INFO )
! 3379: *
! 3380: * Multiply right singular vectors of L in WORK(IU) by
! 3381: * Q in VT, storing result in A
! 3382: * (CWorkspace: need M*M)
! 3383: * (RWorkspace: 0)
! 3384: *
! 3385: CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IU ),
! 3386: $ LDWRKU, VT, LDVT, CZERO, A, LDA )
! 3387: *
! 3388: * Copy right singular vectors of A from A to VT
! 3389: *
! 3390: CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
! 3391: *
! 3392: ELSE
! 3393: *
! 3394: * Insufficient workspace for a fast algorithm
! 3395: *
! 3396: ITAU = 1
! 3397: IWORK = ITAU + M
! 3398: *
! 3399: * Compute A=L*Q, copying result to VT
! 3400: * (CWorkspace: need 2*M, prefer M+M*NB)
! 3401: * (RWorkspace: 0)
! 3402: *
! 3403: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ),
! 3404: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3405: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 3406: *
! 3407: * Generate Q in VT
! 3408: * (CWorkspace: need M+N, prefer M+N*NB)
! 3409: * (RWorkspace: 0)
! 3410: *
! 3411: CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
! 3412: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3413: *
! 3414: * Copy L to U, zeroing out above it
! 3415: *
! 3416: CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
! 3417: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 3418: $ U( 1, 2 ), LDU )
! 3419: IE = 1
! 3420: ITAUQ = ITAU
! 3421: ITAUP = ITAUQ + M
! 3422: IWORK = ITAUP + M
! 3423: *
! 3424: * Bidiagonalize L in U
! 3425: * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
! 3426: * (RWorkspace: need M)
! 3427: *
! 3428: CALL ZGEBRD( M, M, U, LDU, S, RWORK( IE ),
! 3429: $ WORK( ITAUQ ), WORK( ITAUP ),
! 3430: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3431: *
! 3432: * Multiply right bidiagonalizing vectors in U by Q
! 3433: * in VT
! 3434: * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
! 3435: * (RWorkspace: 0)
! 3436: *
! 3437: CALL ZUNMBR( 'P', 'L', 'C', M, N, M, U, LDU,
! 3438: $ WORK( ITAUP ), VT, LDVT,
! 3439: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3440: *
! 3441: * Generate left bidiagonalizing vectors in U
! 3442: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 3443: * (RWorkspace: 0)
! 3444: *
! 3445: CALL ZUNGBR( 'Q', M, M, M, U, LDU, WORK( ITAUQ ),
! 3446: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3447: IRWORK = IE + M
! 3448: *
! 3449: * Perform bidiagonal QR iteration, computing left
! 3450: * singular vectors of A in U and computing right
! 3451: * singular vectors of A in VT
! 3452: * (CWorkspace: 0)
! 3453: * (RWorkspace: need BDSPAC)
! 3454: *
! 3455: CALL ZBDSQR( 'U', M, N, M, 0, S, RWORK( IE ), VT,
! 3456: $ LDVT, U, LDU, CDUM, 1,
! 3457: $ RWORK( IRWORK ), INFO )
! 3458: *
! 3459: END IF
! 3460: *
! 3461: END IF
! 3462: *
! 3463: END IF
! 3464: *
! 3465: ELSE
! 3466: *
! 3467: * N .LT. MNTHR
! 3468: *
! 3469: * Path 10t(N greater than M, but not much larger)
! 3470: * Reduce to bidiagonal form without LQ decomposition
! 3471: *
! 3472: IE = 1
! 3473: ITAUQ = 1
! 3474: ITAUP = ITAUQ + M
! 3475: IWORK = ITAUP + M
! 3476: *
! 3477: * Bidiagonalize A
! 3478: * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
! 3479: * (RWorkspace: M)
! 3480: *
! 3481: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 3482: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
! 3483: $ IERR )
! 3484: IF( WNTUAS ) THEN
! 3485: *
! 3486: * If left singular vectors desired in U, copy result to U
! 3487: * and generate left bidiagonalizing vectors in U
! 3488: * (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
! 3489: * (RWorkspace: 0)
! 3490: *
! 3491: CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
! 3492: CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
! 3493: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3494: END IF
! 3495: IF( WNTVAS ) THEN
! 3496: *
! 3497: * If right singular vectors desired in VT, copy result to
! 3498: * VT and generate right bidiagonalizing vectors in VT
! 3499: * (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB)
! 3500: * (RWorkspace: 0)
! 3501: *
! 3502: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 3503: IF( WNTVA )
! 3504: $ NRVT = N
! 3505: IF( WNTVS )
! 3506: $ NRVT = M
! 3507: CALL ZUNGBR( 'P', NRVT, N, M, VT, LDVT, WORK( ITAUP ),
! 3508: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3509: END IF
! 3510: IF( WNTUO ) THEN
! 3511: *
! 3512: * If left singular vectors desired in A, generate left
! 3513: * bidiagonalizing vectors in A
! 3514: * (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
! 3515: * (RWorkspace: 0)
! 3516: *
! 3517: CALL ZUNGBR( 'Q', M, M, N, A, LDA, WORK( ITAUQ ),
! 3518: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3519: END IF
! 3520: IF( WNTVO ) THEN
! 3521: *
! 3522: * If right singular vectors desired in A, generate right
! 3523: * bidiagonalizing vectors in A
! 3524: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 3525: * (RWorkspace: 0)
! 3526: *
! 3527: CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
! 3528: $ WORK( IWORK ), LWORK-IWORK+1, IERR )
! 3529: END IF
! 3530: IRWORK = IE + M
! 3531: IF( WNTUAS .OR. WNTUO )
! 3532: $ NRU = M
! 3533: IF( WNTUN )
! 3534: $ NRU = 0
! 3535: IF( WNTVAS .OR. WNTVO )
! 3536: $ NCVT = N
! 3537: IF( WNTVN )
! 3538: $ NCVT = 0
! 3539: IF( ( .NOT.WNTUO ) .AND. ( .NOT.WNTVO ) ) THEN
! 3540: *
! 3541: * Perform bidiagonal QR iteration, if desired, computing
! 3542: * left singular vectors in U and computing right singular
! 3543: * vectors in VT
! 3544: * (CWorkspace: 0)
! 3545: * (RWorkspace: need BDSPAC)
! 3546: *
! 3547: CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT,
! 3548: $ LDVT, U, LDU, CDUM, 1, RWORK( IRWORK ),
! 3549: $ INFO )
! 3550: ELSE IF( ( .NOT.WNTUO ) .AND. WNTVO ) THEN
! 3551: *
! 3552: * Perform bidiagonal QR iteration, if desired, computing
! 3553: * left singular vectors in U and computing right singular
! 3554: * vectors in A
! 3555: * (CWorkspace: 0)
! 3556: * (RWorkspace: need BDSPAC)
! 3557: *
! 3558: CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), A,
! 3559: $ LDA, U, LDU, CDUM, 1, RWORK( IRWORK ),
! 3560: $ INFO )
! 3561: ELSE
! 3562: *
! 3563: * Perform bidiagonal QR iteration, if desired, computing
! 3564: * left singular vectors in A and computing right singular
! 3565: * vectors in VT
! 3566: * (CWorkspace: 0)
! 3567: * (RWorkspace: need BDSPAC)
! 3568: *
! 3569: CALL ZBDSQR( 'L', M, NCVT, NRU, 0, S, RWORK( IE ), VT,
! 3570: $ LDVT, A, LDA, CDUM, 1, RWORK( IRWORK ),
! 3571: $ INFO )
! 3572: END IF
! 3573: *
! 3574: END IF
! 3575: *
! 3576: END IF
! 3577: *
! 3578: * Undo scaling if necessary
! 3579: *
! 3580: IF( ISCL.EQ.1 ) THEN
! 3581: IF( ANRM.GT.BIGNUM )
! 3582: $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
! 3583: $ IERR )
! 3584: IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
! 3585: $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1,
! 3586: $ RWORK( IE ), MINMN, IERR )
! 3587: IF( ANRM.LT.SMLNUM )
! 3588: $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
! 3589: $ IERR )
! 3590: IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
! 3591: $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1,
! 3592: $ RWORK( IE ), MINMN, IERR )
! 3593: END IF
! 3594: *
! 3595: * Return optimal workspace in WORK(1)
! 3596: *
! 3597: WORK( 1 ) = MAXWRK
! 3598: *
! 3599: RETURN
! 3600: *
! 3601: * End of ZGESVD
! 3602: *
! 3603: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgesvd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack