Annotation of rpl/lapack/lapack/zgesdd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
! 2: $ LWORK, RWORK, IWORK, 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: * 8-15-00: Improve consistency of WS calculations (eca)
! 9: *
! 10: * .. Scalar Arguments ..
! 11: CHARACTER JOBZ
! 12: INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
! 13: * ..
! 14: * .. Array Arguments ..
! 15: INTEGER IWORK( * )
! 16: DOUBLE PRECISION RWORK( * ), S( * )
! 17: COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
! 18: $ WORK( * )
! 19: * ..
! 20: *
! 21: * Purpose
! 22: * =======
! 23: *
! 24: * ZGESDD computes the singular value decomposition (SVD) of a complex
! 25: * M-by-N matrix A, optionally computing the left and/or right singular
! 26: * vectors, by using divide-and-conquer method. The SVD is written
! 27: *
! 28: * A = U * SIGMA * conjugate-transpose(V)
! 29: *
! 30: * where SIGMA is an M-by-N matrix which is zero except for its
! 31: * min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
! 32: * V is an N-by-N unitary matrix. The diagonal elements of SIGMA
! 33: * are the singular values of A; they are real and non-negative, and
! 34: * are returned in descending order. The first min(m,n) columns of
! 35: * U and V are the left and right singular vectors of A.
! 36: *
! 37: * Note that the routine returns VT = V**H, not V.
! 38: *
! 39: * The divide and conquer algorithm makes very mild assumptions about
! 40: * floating point arithmetic. It will work on machines with a guard
! 41: * digit in add/subtract, or on those binary machines without guard
! 42: * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
! 43: * Cray-2. It could conceivably fail on hexadecimal or decimal machines
! 44: * without guard digits, but we know of none.
! 45: *
! 46: * Arguments
! 47: * =========
! 48: *
! 49: * JOBZ (input) CHARACTER*1
! 50: * Specifies options for computing all or part of the matrix U:
! 51: * = 'A': all M columns of U and all N rows of V**H are
! 52: * returned in the arrays U and VT;
! 53: * = 'S': the first min(M,N) columns of U and the first
! 54: * min(M,N) rows of V**H are returned in the arrays U
! 55: * and VT;
! 56: * = 'O': If M >= N, the first N columns of U are overwritten
! 57: * in the array A and all rows of V**H are returned in
! 58: * the array VT;
! 59: * otherwise, all columns of U are returned in the
! 60: * array U and the first M rows of V**H are overwritten
! 61: * in the array A;
! 62: * = 'N': no columns of U or rows of V**H are computed.
! 63: *
! 64: * M (input) INTEGER
! 65: * The number of rows of the input matrix A. M >= 0.
! 66: *
! 67: * N (input) INTEGER
! 68: * The number of columns of the input matrix A. N >= 0.
! 69: *
! 70: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
! 71: * On entry, the M-by-N matrix A.
! 72: * On exit,
! 73: * if JOBZ = 'O', A is overwritten with the first N columns
! 74: * of U (the left singular vectors, stored
! 75: * columnwise) if M >= N;
! 76: * A is overwritten with the first M rows
! 77: * of V**H (the right singular vectors, stored
! 78: * rowwise) otherwise.
! 79: * if JOBZ .ne. 'O', the contents of A 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: * UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
! 89: * UCOL = min(M,N) if JOBZ = 'S'.
! 90: * If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
! 91: * unitary matrix U;
! 92: * if JOBZ = 'S', U contains the first min(M,N) columns of U
! 93: * (the left singular vectors, stored columnwise);
! 94: * if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
! 95: *
! 96: * LDU (input) INTEGER
! 97: * The leading dimension of the array U. LDU >= 1; if
! 98: * JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
! 99: *
! 100: * VT (output) COMPLEX*16 array, dimension (LDVT,N)
! 101: * If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
! 102: * N-by-N unitary matrix V**H;
! 103: * if JOBZ = 'S', VT contains the first min(M,N) rows of
! 104: * V**H (the right singular vectors, stored rowwise);
! 105: * if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
! 106: *
! 107: * LDVT (input) INTEGER
! 108: * The leading dimension of the array VT. LDVT >= 1; if
! 109: * JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
! 110: * if JOBZ = 'S', LDVT >= min(M,N).
! 111: *
! 112: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
! 113: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 114: *
! 115: * LWORK (input) INTEGER
! 116: * The dimension of the array WORK. LWORK >= 1.
! 117: * if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N).
! 118: * if JOBZ = 'O',
! 119: * LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N).
! 120: * if JOBZ = 'S' or 'A',
! 121: * LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N).
! 122: * For good performance, LWORK should generally be larger.
! 123: *
! 124: * If LWORK = -1, a workspace query is assumed. The optimal
! 125: * size for the WORK array is calculated and stored in WORK(1),
! 126: * and no other work except argument checking is performed.
! 127: *
! 128: * RWORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
! 129: * If JOBZ = 'N', LRWORK >= 5*min(M,N).
! 130: * Otherwise, LRWORK >= 5*min(M,N)*min(M,N) + 7*min(M,N)
! 131: *
! 132: * IWORK (workspace) INTEGER array, dimension (8*min(M,N))
! 133: *
! 134: * INFO (output) INTEGER
! 135: * = 0: successful exit.
! 136: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 137: * > 0: The updating process of DBDSDC did not converge.
! 138: *
! 139: * Further Details
! 140: * ===============
! 141: *
! 142: * Based on contributions by
! 143: * Ming Gu and Huan Ren, Computer Science Division, University of
! 144: * California at Berkeley, USA
! 145: *
! 146: * =====================================================================
! 147: *
! 148: * .. Parameters ..
! 149: INTEGER LQUERV
! 150: PARAMETER ( LQUERV = -1 )
! 151: COMPLEX*16 CZERO, CONE
! 152: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
! 153: $ CONE = ( 1.0D+0, 0.0D+0 ) )
! 154: DOUBLE PRECISION ZERO, ONE
! 155: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 156: * ..
! 157: * .. Local Scalars ..
! 158: LOGICAL WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
! 159: INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
! 160: $ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
! 161: $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
! 162: $ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL
! 163: DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
! 164: * ..
! 165: * .. Local Arrays ..
! 166: INTEGER IDUM( 1 )
! 167: DOUBLE PRECISION DUM( 1 )
! 168: * ..
! 169: * .. External Subroutines ..
! 170: EXTERNAL DBDSDC, DLASCL, XERBLA, ZGEBRD, ZGELQF, ZGEMM,
! 171: $ ZGEQRF, ZLACP2, ZLACPY, ZLACRM, ZLARCM, ZLASCL,
! 172: $ ZLASET, ZUNGBR, ZUNGLQ, ZUNGQR, ZUNMBR
! 173: * ..
! 174: * .. External Functions ..
! 175: LOGICAL LSAME
! 176: INTEGER ILAENV
! 177: DOUBLE PRECISION DLAMCH, ZLANGE
! 178: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
! 179: * ..
! 180: * .. Intrinsic Functions ..
! 181: INTRINSIC INT, MAX, MIN, SQRT
! 182: * ..
! 183: * .. Executable Statements ..
! 184: *
! 185: * Test the input arguments
! 186: *
! 187: INFO = 0
! 188: MINMN = MIN( M, N )
! 189: MNTHR1 = INT( MINMN*17.0D0 / 9.0D0 )
! 190: MNTHR2 = INT( MINMN*5.0D0 / 3.0D0 )
! 191: WNTQA = LSAME( JOBZ, 'A' )
! 192: WNTQS = LSAME( JOBZ, 'S' )
! 193: WNTQAS = WNTQA .OR. WNTQS
! 194: WNTQO = LSAME( JOBZ, 'O' )
! 195: WNTQN = LSAME( JOBZ, 'N' )
! 196: MINWRK = 1
! 197: MAXWRK = 1
! 198: *
! 199: IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
! 200: INFO = -1
! 201: ELSE IF( M.LT.0 ) THEN
! 202: INFO = -2
! 203: ELSE IF( N.LT.0 ) THEN
! 204: INFO = -3
! 205: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 206: INFO = -5
! 207: ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
! 208: $ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
! 209: INFO = -8
! 210: ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
! 211: $ ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
! 212: $ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
! 213: INFO = -10
! 214: END IF
! 215: *
! 216: * Compute workspace
! 217: * (Note: Comments in the code beginning "Workspace:" describe the
! 218: * minimal amount of workspace needed at that point in the code,
! 219: * as well as the preferred amount for good performance.
! 220: * CWorkspace refers to complex workspace, and RWorkspace to
! 221: * real workspace. NB refers to the optimal block size for the
! 222: * immediately following subroutine, as returned by ILAENV.)
! 223: *
! 224: IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN
! 225: IF( M.GE.N ) THEN
! 226: *
! 227: * There is no complex work space needed for bidiagonal SVD
! 228: * The real work space needed for bidiagonal SVD is BDSPAC
! 229: * for computing singular values and singular vectors; BDSPAN
! 230: * for computing singular values only.
! 231: * BDSPAC = 5*N*N + 7*N
! 232: * BDSPAN = MAX(7*N+4, 3*N+2+SMLSIZ*(SMLSIZ+8))
! 233: *
! 234: IF( M.GE.MNTHR1 ) THEN
! 235: IF( WNTQN ) THEN
! 236: *
! 237: * Path 1 (M much larger than N, JOBZ='N')
! 238: *
! 239: MAXWRK = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1,
! 240: $ -1 )
! 241: MAXWRK = MAX( MAXWRK, 2*N+2*N*
! 242: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 243: MINWRK = 3*N
! 244: ELSE IF( WNTQO ) THEN
! 245: *
! 246: * Path 2 (M much larger than N, JOBZ='O')
! 247: *
! 248: WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 249: WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
! 250: $ N, N, -1 ) )
! 251: WRKBL = MAX( WRKBL, 2*N+2*N*
! 252: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 253: WRKBL = MAX( WRKBL, 2*N+N*
! 254: $ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )
! 255: WRKBL = MAX( WRKBL, 2*N+N*
! 256: $ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )
! 257: MAXWRK = M*N + N*N + WRKBL
! 258: MINWRK = 2*N*N + 3*N
! 259: ELSE IF( WNTQS ) THEN
! 260: *
! 261: * Path 3 (M much larger than N, JOBZ='S')
! 262: *
! 263: WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 264: WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'ZUNGQR', ' ', M,
! 265: $ N, N, -1 ) )
! 266: WRKBL = MAX( WRKBL, 2*N+2*N*
! 267: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 268: WRKBL = MAX( WRKBL, 2*N+N*
! 269: $ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )
! 270: WRKBL = MAX( WRKBL, 2*N+N*
! 271: $ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )
! 272: MAXWRK = N*N + WRKBL
! 273: MINWRK = N*N + 3*N
! 274: ELSE IF( WNTQA ) THEN
! 275: *
! 276: * Path 4 (M much larger than N, JOBZ='A')
! 277: *
! 278: WRKBL = N + N*ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 279: WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'ZUNGQR', ' ', M,
! 280: $ M, N, -1 ) )
! 281: WRKBL = MAX( WRKBL, 2*N+2*N*
! 282: $ ILAENV( 1, 'ZGEBRD', ' ', N, N, -1, -1 ) )
! 283: WRKBL = MAX( WRKBL, 2*N+N*
! 284: $ ILAENV( 1, 'ZUNMBR', 'QLN', N, N, N, -1 ) )
! 285: WRKBL = MAX( WRKBL, 2*N+N*
! 286: $ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )
! 287: MAXWRK = N*N + WRKBL
! 288: MINWRK = N*N + 2*N + M
! 289: END IF
! 290: ELSE IF( M.GE.MNTHR2 ) THEN
! 291: *
! 292: * Path 5 (M much larger than N, but not as much as MNTHR1)
! 293: *
! 294: MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,
! 295: $ -1, -1 )
! 296: MINWRK = 2*N + M
! 297: IF( WNTQO ) THEN
! 298: MAXWRK = MAX( MAXWRK, 2*N+N*
! 299: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
! 300: MAXWRK = MAX( MAXWRK, 2*N+N*
! 301: $ ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) )
! 302: MAXWRK = MAXWRK + M*N
! 303: MINWRK = MINWRK + N*N
! 304: ELSE IF( WNTQS ) THEN
! 305: MAXWRK = MAX( MAXWRK, 2*N+N*
! 306: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
! 307: MAXWRK = MAX( MAXWRK, 2*N+N*
! 308: $ ILAENV( 1, 'ZUNGBR', 'Q', M, N, N, -1 ) )
! 309: ELSE IF( WNTQA ) THEN
! 310: MAXWRK = MAX( MAXWRK, 2*N+N*
! 311: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, N, -1 ) )
! 312: MAXWRK = MAX( MAXWRK, 2*N+M*
! 313: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )
! 314: END IF
! 315: ELSE
! 316: *
! 317: * Path 6 (M at least N, but not much larger)
! 318: *
! 319: MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,
! 320: $ -1, -1 )
! 321: MINWRK = 2*N + M
! 322: IF( WNTQO ) THEN
! 323: MAXWRK = MAX( MAXWRK, 2*N+N*
! 324: $ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )
! 325: MAXWRK = MAX( MAXWRK, 2*N+N*
! 326: $ ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 ) )
! 327: MAXWRK = MAXWRK + M*N
! 328: MINWRK = MINWRK + N*N
! 329: ELSE IF( WNTQS ) THEN
! 330: MAXWRK = MAX( MAXWRK, 2*N+N*
! 331: $ ILAENV( 1, 'ZUNMBR', 'PRC', N, N, N, -1 ) )
! 332: MAXWRK = MAX( MAXWRK, 2*N+N*
! 333: $ ILAENV( 1, 'ZUNMBR', 'QLN', M, N, N, -1 ) )
! 334: ELSE IF( WNTQA ) THEN
! 335: MAXWRK = MAX( MAXWRK, 2*N+N*
! 336: $ ILAENV( 1, 'ZUNGBR', 'PRC', N, N, N, -1 ) )
! 337: MAXWRK = MAX( MAXWRK, 2*N+M*
! 338: $ ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) )
! 339: END IF
! 340: END IF
! 341: ELSE
! 342: *
! 343: * There is no complex work space needed for bidiagonal SVD
! 344: * The real work space needed for bidiagonal SVD is BDSPAC
! 345: * for computing singular values and singular vectors; BDSPAN
! 346: * for computing singular values only.
! 347: * BDSPAC = 5*M*M + 7*M
! 348: * BDSPAN = MAX(7*M+4, 3*M+2+SMLSIZ*(SMLSIZ+8))
! 349: *
! 350: IF( N.GE.MNTHR1 ) THEN
! 351: IF( WNTQN ) THEN
! 352: *
! 353: * Path 1t (N much larger than M, JOBZ='N')
! 354: *
! 355: MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
! 356: $ -1 )
! 357: MAXWRK = MAX( MAXWRK, 2*M+2*M*
! 358: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 359: MINWRK = 3*M
! 360: ELSE IF( WNTQO ) THEN
! 361: *
! 362: * Path 2t (N much larger than M, JOBZ='O')
! 363: *
! 364: WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 365: WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
! 366: $ N, M, -1 ) )
! 367: WRKBL = MAX( WRKBL, 2*M+2*M*
! 368: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 369: WRKBL = MAX( WRKBL, 2*M+M*
! 370: $ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )
! 371: WRKBL = MAX( WRKBL, 2*M+M*
! 372: $ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )
! 373: MAXWRK = M*N + M*M + WRKBL
! 374: MINWRK = 2*M*M + 3*M
! 375: ELSE IF( WNTQS ) THEN
! 376: *
! 377: * Path 3t (N much larger than M, JOBZ='S')
! 378: *
! 379: WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 380: WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'ZUNGLQ', ' ', M,
! 381: $ N, M, -1 ) )
! 382: WRKBL = MAX( WRKBL, 2*M+2*M*
! 383: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 384: WRKBL = MAX( WRKBL, 2*M+M*
! 385: $ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )
! 386: WRKBL = MAX( WRKBL, 2*M+M*
! 387: $ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )
! 388: MAXWRK = M*M + WRKBL
! 389: MINWRK = M*M + 3*M
! 390: ELSE IF( WNTQA ) THEN
! 391: *
! 392: * Path 4t (N much larger than M, JOBZ='A')
! 393: *
! 394: WRKBL = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 395: WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'ZUNGLQ', ' ', N,
! 396: $ N, M, -1 ) )
! 397: WRKBL = MAX( WRKBL, 2*M+2*M*
! 398: $ ILAENV( 1, 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 399: WRKBL = MAX( WRKBL, 2*M+M*
! 400: $ ILAENV( 1, 'ZUNMBR', 'PRC', M, M, M, -1 ) )
! 401: WRKBL = MAX( WRKBL, 2*M+M*
! 402: $ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, M, -1 ) )
! 403: MAXWRK = M*M + WRKBL
! 404: MINWRK = M*M + 2*M + N
! 405: END IF
! 406: ELSE IF( N.GE.MNTHR2 ) THEN
! 407: *
! 408: * Path 5t (N much larger than M, but not as much as MNTHR1)
! 409: *
! 410: MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,
! 411: $ -1, -1 )
! 412: MINWRK = 2*M + N
! 413: IF( WNTQO ) THEN
! 414: MAXWRK = MAX( MAXWRK, 2*M+M*
! 415: $ ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) )
! 416: MAXWRK = MAX( MAXWRK, 2*M+M*
! 417: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )
! 418: MAXWRK = MAXWRK + M*N
! 419: MINWRK = MINWRK + M*M
! 420: ELSE IF( WNTQS ) THEN
! 421: MAXWRK = MAX( MAXWRK, 2*M+M*
! 422: $ ILAENV( 1, 'ZUNGBR', 'P', M, N, M, -1 ) )
! 423: MAXWRK = MAX( MAXWRK, 2*M+M*
! 424: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )
! 425: ELSE IF( WNTQA ) THEN
! 426: MAXWRK = MAX( MAXWRK, 2*M+N*
! 427: $ ILAENV( 1, 'ZUNGBR', 'P', N, N, M, -1 ) )
! 428: MAXWRK = MAX( MAXWRK, 2*M+M*
! 429: $ ILAENV( 1, 'ZUNGBR', 'Q', M, M, N, -1 ) )
! 430: END IF
! 431: ELSE
! 432: *
! 433: * Path 6t (N greater than M, but not much larger)
! 434: *
! 435: MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'ZGEBRD', ' ', M, N,
! 436: $ -1, -1 )
! 437: MINWRK = 2*M + N
! 438: IF( WNTQO ) THEN
! 439: MAXWRK = MAX( MAXWRK, 2*M+M*
! 440: $ ILAENV( 1, 'ZUNMBR', 'PRC', M, N, M, -1 ) )
! 441: MAXWRK = MAX( MAXWRK, 2*M+M*
! 442: $ ILAENV( 1, 'ZUNMBR', 'QLN', M, M, N, -1 ) )
! 443: MAXWRK = MAXWRK + M*N
! 444: MINWRK = MINWRK + M*M
! 445: ELSE IF( WNTQS ) THEN
! 446: MAXWRK = MAX( MAXWRK, 2*M+M*
! 447: $ ILAENV( 1, 'ZUNGBR', 'PRC', M, N, M, -1 ) )
! 448: MAXWRK = MAX( MAXWRK, 2*M+M*
! 449: $ ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) )
! 450: ELSE IF( WNTQA ) THEN
! 451: MAXWRK = MAX( MAXWRK, 2*M+N*
! 452: $ ILAENV( 1, 'ZUNGBR', 'PRC', N, N, M, -1 ) )
! 453: MAXWRK = MAX( MAXWRK, 2*M+M*
! 454: $ ILAENV( 1, 'ZUNGBR', 'QLN', M, M, N, -1 ) )
! 455: END IF
! 456: END IF
! 457: END IF
! 458: MAXWRK = MAX( MAXWRK, MINWRK )
! 459: END IF
! 460: IF( INFO.EQ.0 ) THEN
! 461: WORK( 1 ) = MAXWRK
! 462: IF( LWORK.LT.MINWRK .AND. LWORK.NE.LQUERV )
! 463: $ INFO = -13
! 464: END IF
! 465: *
! 466: * Quick returns
! 467: *
! 468: IF( INFO.NE.0 ) THEN
! 469: CALL XERBLA( 'ZGESDD', -INFO )
! 470: RETURN
! 471: END IF
! 472: IF( LWORK.EQ.LQUERV )
! 473: $ RETURN
! 474: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
! 475: RETURN
! 476: END IF
! 477: *
! 478: * Get machine constants
! 479: *
! 480: EPS = DLAMCH( 'P' )
! 481: SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
! 482: BIGNUM = ONE / SMLNUM
! 483: *
! 484: * Scale A if max element outside range [SMLNUM,BIGNUM]
! 485: *
! 486: ANRM = ZLANGE( 'M', M, N, A, LDA, DUM )
! 487: ISCL = 0
! 488: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 489: ISCL = 1
! 490: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
! 491: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 492: ISCL = 1
! 493: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
! 494: END IF
! 495: *
! 496: IF( M.GE.N ) THEN
! 497: *
! 498: * A has at least as many rows as columns. If A has sufficiently
! 499: * more rows than columns, first reduce using the QR
! 500: * decomposition (if sufficient workspace available)
! 501: *
! 502: IF( M.GE.MNTHR1 ) THEN
! 503: *
! 504: IF( WNTQN ) THEN
! 505: *
! 506: * Path 1 (M much larger than N, JOBZ='N')
! 507: * No singular vectors to be computed
! 508: *
! 509: ITAU = 1
! 510: NWORK = ITAU + N
! 511: *
! 512: * Compute A=Q*R
! 513: * (CWorkspace: need 2*N, prefer N+N*NB)
! 514: * (RWorkspace: need 0)
! 515: *
! 516: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 517: $ LWORK-NWORK+1, IERR )
! 518: *
! 519: * Zero out below R
! 520: *
! 521: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
! 522: $ LDA )
! 523: IE = 1
! 524: ITAUQ = 1
! 525: ITAUP = ITAUQ + N
! 526: NWORK = ITAUP + N
! 527: *
! 528: * Bidiagonalize R in A
! 529: * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
! 530: * (RWorkspace: need N)
! 531: *
! 532: CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 533: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 534: $ IERR )
! 535: NRWORK = IE + N
! 536: *
! 537: * Perform bidiagonal SVD, compute singular values only
! 538: * (CWorkspace: 0)
! 539: * (RWorkspace: need BDSPAN)
! 540: *
! 541: CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
! 542: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
! 543: *
! 544: ELSE IF( WNTQO ) THEN
! 545: *
! 546: * Path 2 (M much larger than N, JOBZ='O')
! 547: * N left singular vectors to be overwritten on A and
! 548: * N right singular vectors to be computed in VT
! 549: *
! 550: IU = 1
! 551: *
! 552: * WORK(IU) is N by N
! 553: *
! 554: LDWRKU = N
! 555: IR = IU + LDWRKU*N
! 556: IF( LWORK.GE.M*N+N*N+3*N ) THEN
! 557: *
! 558: * WORK(IR) is M by N
! 559: *
! 560: LDWRKR = M
! 561: ELSE
! 562: LDWRKR = ( LWORK-N*N-3*N ) / N
! 563: END IF
! 564: ITAU = IR + LDWRKR*N
! 565: NWORK = ITAU + N
! 566: *
! 567: * Compute A=Q*R
! 568: * (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB)
! 569: * (RWorkspace: 0)
! 570: *
! 571: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 572: $ LWORK-NWORK+1, IERR )
! 573: *
! 574: * Copy R to WORK( IR ), zeroing out below it
! 575: *
! 576: CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
! 577: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
! 578: $ LDWRKR )
! 579: *
! 580: * Generate Q in A
! 581: * (CWorkspace: need 2*N, prefer N+N*NB)
! 582: * (RWorkspace: 0)
! 583: *
! 584: CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
! 585: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 586: IE = 1
! 587: ITAUQ = ITAU
! 588: ITAUP = ITAUQ + N
! 589: NWORK = ITAUP + N
! 590: *
! 591: * Bidiagonalize R in WORK(IR)
! 592: * (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB)
! 593: * (RWorkspace: need N)
! 594: *
! 595: CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
! 596: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
! 597: $ LWORK-NWORK+1, IERR )
! 598: *
! 599: * Perform bidiagonal SVD, computing left singular vectors
! 600: * of R in WORK(IRU) and computing right singular vectors
! 601: * of R in WORK(IRVT)
! 602: * (CWorkspace: need 0)
! 603: * (RWorkspace: need BDSPAC)
! 604: *
! 605: IRU = IE + N
! 606: IRVT = IRU + N*N
! 607: NRWORK = IRVT + N*N
! 608: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
! 609: $ N, RWORK( IRVT ), N, DUM, IDUM,
! 610: $ RWORK( NRWORK ), IWORK, INFO )
! 611: *
! 612: * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
! 613: * Overwrite WORK(IU) by the left singular vectors of R
! 614: * (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB)
! 615: * (RWorkspace: 0)
! 616: *
! 617: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
! 618: $ LDWRKU )
! 619: CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
! 620: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
! 621: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 622: *
! 623: * Copy real matrix RWORK(IRVT) to complex matrix VT
! 624: * Overwrite VT by the right singular vectors of R
! 625: * (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
! 626: * (RWorkspace: 0)
! 627: *
! 628: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
! 629: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
! 630: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 631: $ LWORK-NWORK+1, IERR )
! 632: *
! 633: * Multiply Q in A by left singular vectors of R in
! 634: * WORK(IU), storing result in WORK(IR) and copying to A
! 635: * (CWorkspace: need 2*N*N, prefer N*N+M*N)
! 636: * (RWorkspace: 0)
! 637: *
! 638: DO 10 I = 1, M, LDWRKR
! 639: CHUNK = MIN( M-I+1, LDWRKR )
! 640: CALL ZGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
! 641: $ LDA, WORK( IU ), LDWRKU, CZERO,
! 642: $ WORK( IR ), LDWRKR )
! 643: CALL ZLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
! 644: $ A( I, 1 ), LDA )
! 645: 10 CONTINUE
! 646: *
! 647: ELSE IF( WNTQS ) THEN
! 648: *
! 649: * Path 3 (M much larger than N, JOBZ='S')
! 650: * N left singular vectors to be computed in U and
! 651: * N right singular vectors to be computed in VT
! 652: *
! 653: IR = 1
! 654: *
! 655: * WORK(IR) is N by N
! 656: *
! 657: LDWRKR = N
! 658: ITAU = IR + LDWRKR*N
! 659: NWORK = ITAU + N
! 660: *
! 661: * Compute A=Q*R
! 662: * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
! 663: * (RWorkspace: 0)
! 664: *
! 665: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 666: $ LWORK-NWORK+1, IERR )
! 667: *
! 668: * Copy R to WORK(IR), zeroing out below it
! 669: *
! 670: CALL ZLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
! 671: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
! 672: $ LDWRKR )
! 673: *
! 674: * Generate Q in A
! 675: * (CWorkspace: need 2*N, prefer N+N*NB)
! 676: * (RWorkspace: 0)
! 677: *
! 678: CALL ZUNGQR( M, N, N, A, LDA, WORK( ITAU ),
! 679: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 680: IE = 1
! 681: ITAUQ = ITAU
! 682: ITAUP = ITAUQ + N
! 683: NWORK = ITAUP + N
! 684: *
! 685: * Bidiagonalize R in WORK(IR)
! 686: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
! 687: * (RWorkspace: need N)
! 688: *
! 689: CALL ZGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
! 690: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
! 691: $ LWORK-NWORK+1, IERR )
! 692: *
! 693: * Perform bidiagonal SVD, computing left singular vectors
! 694: * of bidiagonal matrix in RWORK(IRU) and computing right
! 695: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 696: * (CWorkspace: need 0)
! 697: * (RWorkspace: need BDSPAC)
! 698: *
! 699: IRU = IE + N
! 700: IRVT = IRU + N*N
! 701: NRWORK = IRVT + N*N
! 702: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
! 703: $ N, RWORK( IRVT ), N, DUM, IDUM,
! 704: $ RWORK( NRWORK ), IWORK, INFO )
! 705: *
! 706: * Copy real matrix RWORK(IRU) to complex matrix U
! 707: * Overwrite U by left singular vectors of R
! 708: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
! 709: * (RWorkspace: 0)
! 710: *
! 711: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
! 712: CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
! 713: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 714: $ LWORK-NWORK+1, IERR )
! 715: *
! 716: * Copy real matrix RWORK(IRVT) to complex matrix VT
! 717: * Overwrite VT by right singular vectors of R
! 718: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
! 719: * (RWorkspace: 0)
! 720: *
! 721: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
! 722: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
! 723: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 724: $ LWORK-NWORK+1, IERR )
! 725: *
! 726: * Multiply Q in A by left singular vectors of R in
! 727: * WORK(IR), storing result in U
! 728: * (CWorkspace: need N*N)
! 729: * (RWorkspace: 0)
! 730: *
! 731: CALL ZLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
! 732: CALL ZGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),
! 733: $ LDWRKR, CZERO, U, LDU )
! 734: *
! 735: ELSE IF( WNTQA ) THEN
! 736: *
! 737: * Path 4 (M much larger than N, JOBZ='A')
! 738: * M left singular vectors to be computed in U and
! 739: * N right singular vectors to be computed in VT
! 740: *
! 741: IU = 1
! 742: *
! 743: * WORK(IU) is N by N
! 744: *
! 745: LDWRKU = N
! 746: ITAU = IU + LDWRKU*N
! 747: NWORK = ITAU + N
! 748: *
! 749: * Compute A=Q*R, copying result to U
! 750: * (CWorkspace: need 2*N, prefer N+N*NB)
! 751: * (RWorkspace: 0)
! 752: *
! 753: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 754: $ LWORK-NWORK+1, IERR )
! 755: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 756: *
! 757: * Generate Q in U
! 758: * (CWorkspace: need N+M, prefer N+M*NB)
! 759: * (RWorkspace: 0)
! 760: *
! 761: CALL ZUNGQR( M, M, N, U, LDU, WORK( ITAU ),
! 762: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 763: *
! 764: * Produce R in A, zeroing out below it
! 765: *
! 766: CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
! 767: $ LDA )
! 768: IE = 1
! 769: ITAUQ = ITAU
! 770: ITAUP = ITAUQ + N
! 771: NWORK = ITAUP + N
! 772: *
! 773: * Bidiagonalize R in A
! 774: * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
! 775: * (RWorkspace: need N)
! 776: *
! 777: CALL ZGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 778: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 779: $ IERR )
! 780: IRU = IE + N
! 781: IRVT = IRU + N*N
! 782: NRWORK = IRVT + N*N
! 783: *
! 784: * Perform bidiagonal SVD, computing left singular vectors
! 785: * of bidiagonal matrix in RWORK(IRU) and computing right
! 786: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 787: * (CWorkspace: need 0)
! 788: * (RWorkspace: need BDSPAC)
! 789: *
! 790: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
! 791: $ N, RWORK( IRVT ), N, DUM, IDUM,
! 792: $ RWORK( NRWORK ), IWORK, INFO )
! 793: *
! 794: * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
! 795: * Overwrite WORK(IU) by left singular vectors of R
! 796: * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
! 797: * (RWorkspace: 0)
! 798: *
! 799: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
! 800: $ LDWRKU )
! 801: CALL ZUNMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
! 802: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
! 803: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 804: *
! 805: * Copy real matrix RWORK(IRVT) to complex matrix VT
! 806: * Overwrite VT by right singular vectors of R
! 807: * (CWorkspace: need 3*N, prefer 2*N+N*NB)
! 808: * (RWorkspace: 0)
! 809: *
! 810: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
! 811: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
! 812: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 813: $ LWORK-NWORK+1, IERR )
! 814: *
! 815: * Multiply Q in U by left singular vectors of R in
! 816: * WORK(IU), storing result in A
! 817: * (CWorkspace: need N*N)
! 818: * (RWorkspace: 0)
! 819: *
! 820: CALL ZGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),
! 821: $ LDWRKU, CZERO, A, LDA )
! 822: *
! 823: * Copy left singular vectors of A from A to U
! 824: *
! 825: CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
! 826: *
! 827: END IF
! 828: *
! 829: ELSE IF( M.GE.MNTHR2 ) THEN
! 830: *
! 831: * MNTHR2 <= M < MNTHR1
! 832: *
! 833: * Path 5 (M much larger than N, but not as much as MNTHR1)
! 834: * Reduce to bidiagonal form without QR decomposition, use
! 835: * ZUNGBR and matrix multiplication to compute singular vectors
! 836: *
! 837: IE = 1
! 838: NRWORK = IE + N
! 839: ITAUQ = 1
! 840: ITAUP = ITAUQ + N
! 841: NWORK = ITAUP + N
! 842: *
! 843: * Bidiagonalize A
! 844: * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
! 845: * (RWorkspace: need N)
! 846: *
! 847: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 848: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 849: $ IERR )
! 850: IF( WNTQN ) THEN
! 851: *
! 852: * Compute singular values only
! 853: * (Cworkspace: 0)
! 854: * (Rworkspace: need BDSPAN)
! 855: *
! 856: CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
! 857: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
! 858: ELSE IF( WNTQO ) THEN
! 859: IU = NWORK
! 860: IRU = NRWORK
! 861: IRVT = IRU + N*N
! 862: NRWORK = IRVT + N*N
! 863: *
! 864: * Copy A to VT, generate P**H
! 865: * (Cworkspace: need 2*N, prefer N+N*NB)
! 866: * (Rworkspace: 0)
! 867: *
! 868: CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
! 869: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
! 870: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 871: *
! 872: * Generate Q in A
! 873: * (CWorkspace: need 2*N, prefer N+N*NB)
! 874: * (RWorkspace: 0)
! 875: *
! 876: CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
! 877: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 878: *
! 879: IF( LWORK.GE.M*N+3*N ) THEN
! 880: *
! 881: * WORK( IU ) is M by N
! 882: *
! 883: LDWRKU = M
! 884: ELSE
! 885: *
! 886: * WORK(IU) is LDWRKU by N
! 887: *
! 888: LDWRKU = ( LWORK-3*N ) / N
! 889: END IF
! 890: NWORK = IU + LDWRKU*N
! 891: *
! 892: * Perform bidiagonal SVD, computing left singular vectors
! 893: * of bidiagonal matrix in RWORK(IRU) and computing right
! 894: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 895: * (CWorkspace: need 0)
! 896: * (RWorkspace: need BDSPAC)
! 897: *
! 898: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
! 899: $ N, RWORK( IRVT ), N, DUM, IDUM,
! 900: $ RWORK( NRWORK ), IWORK, INFO )
! 901: *
! 902: * Multiply real matrix RWORK(IRVT) by P**H in VT,
! 903: * storing the result in WORK(IU), copying to VT
! 904: * (Cworkspace: need 0)
! 905: * (Rworkspace: need 3*N*N)
! 906: *
! 907: CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,
! 908: $ WORK( IU ), LDWRKU, RWORK( NRWORK ) )
! 909: CALL ZLACPY( 'F', N, N, WORK( IU ), LDWRKU, VT, LDVT )
! 910: *
! 911: * Multiply Q in A by real matrix RWORK(IRU), storing the
! 912: * result in WORK(IU), copying to A
! 913: * (CWorkspace: need N*N, prefer M*N)
! 914: * (Rworkspace: need 3*N*N, prefer N*N+2*M*N)
! 915: *
! 916: NRWORK = IRVT
! 917: DO 20 I = 1, M, LDWRKU
! 918: CHUNK = MIN( M-I+1, LDWRKU )
! 919: CALL ZLACRM( CHUNK, N, A( I, 1 ), LDA, RWORK( IRU ),
! 920: $ N, WORK( IU ), LDWRKU, RWORK( NRWORK ) )
! 921: CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
! 922: $ A( I, 1 ), LDA )
! 923: 20 CONTINUE
! 924: *
! 925: ELSE IF( WNTQS ) THEN
! 926: *
! 927: * Copy A to VT, generate P**H
! 928: * (Cworkspace: need 2*N, prefer N+N*NB)
! 929: * (Rworkspace: 0)
! 930: *
! 931: CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
! 932: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
! 933: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 934: *
! 935: * Copy A to U, generate Q
! 936: * (Cworkspace: need 2*N, prefer N+N*NB)
! 937: * (Rworkspace: 0)
! 938: *
! 939: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 940: CALL ZUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),
! 941: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 942: *
! 943: * Perform bidiagonal SVD, computing left singular vectors
! 944: * of bidiagonal matrix in RWORK(IRU) and computing right
! 945: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 946: * (CWorkspace: need 0)
! 947: * (RWorkspace: need BDSPAC)
! 948: *
! 949: IRU = NRWORK
! 950: IRVT = IRU + N*N
! 951: NRWORK = IRVT + N*N
! 952: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
! 953: $ N, RWORK( IRVT ), N, DUM, IDUM,
! 954: $ RWORK( NRWORK ), IWORK, INFO )
! 955: *
! 956: * Multiply real matrix RWORK(IRVT) by P**H in VT,
! 957: * storing the result in A, copying to VT
! 958: * (Cworkspace: need 0)
! 959: * (Rworkspace: need 3*N*N)
! 960: *
! 961: CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
! 962: $ RWORK( NRWORK ) )
! 963: CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT )
! 964: *
! 965: * Multiply Q in U by real matrix RWORK(IRU), storing the
! 966: * result in A, copying to U
! 967: * (CWorkspace: need 0)
! 968: * (Rworkspace: need N*N+2*M*N)
! 969: *
! 970: NRWORK = IRVT
! 971: CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
! 972: $ RWORK( NRWORK ) )
! 973: CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
! 974: ELSE
! 975: *
! 976: * Copy A to VT, generate P**H
! 977: * (Cworkspace: need 2*N, prefer N+N*NB)
! 978: * (Rworkspace: 0)
! 979: *
! 980: CALL ZLACPY( 'U', N, N, A, LDA, VT, LDVT )
! 981: CALL ZUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
! 982: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 983: *
! 984: * Copy A to U, generate Q
! 985: * (Cworkspace: need 2*N, prefer N+N*NB)
! 986: * (Rworkspace: 0)
! 987: *
! 988: CALL ZLACPY( 'L', M, N, A, LDA, U, LDU )
! 989: CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
! 990: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 991: *
! 992: * Perform bidiagonal SVD, computing left singular vectors
! 993: * of bidiagonal matrix in RWORK(IRU) and computing right
! 994: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 995: * (CWorkspace: need 0)
! 996: * (RWorkspace: need BDSPAC)
! 997: *
! 998: IRU = NRWORK
! 999: IRVT = IRU + N*N
! 1000: NRWORK = IRVT + N*N
! 1001: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
! 1002: $ N, RWORK( IRVT ), N, DUM, IDUM,
! 1003: $ RWORK( NRWORK ), IWORK, INFO )
! 1004: *
! 1005: * Multiply real matrix RWORK(IRVT) by P**H in VT,
! 1006: * storing the result in A, copying to VT
! 1007: * (Cworkspace: need 0)
! 1008: * (Rworkspace: need 3*N*N)
! 1009: *
! 1010: CALL ZLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
! 1011: $ RWORK( NRWORK ) )
! 1012: CALL ZLACPY( 'F', N, N, A, LDA, VT, LDVT )
! 1013: *
! 1014: * Multiply Q in U by real matrix RWORK(IRU), storing the
! 1015: * result in A, copying to U
! 1016: * (CWorkspace: 0)
! 1017: * (Rworkspace: need 3*N*N)
! 1018: *
! 1019: NRWORK = IRVT
! 1020: CALL ZLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
! 1021: $ RWORK( NRWORK ) )
! 1022: CALL ZLACPY( 'F', M, N, A, LDA, U, LDU )
! 1023: END IF
! 1024: *
! 1025: ELSE
! 1026: *
! 1027: * M .LT. MNTHR2
! 1028: *
! 1029: * Path 6 (M at least N, but not much larger)
! 1030: * Reduce to bidiagonal form without QR decomposition
! 1031: * Use ZUNMBR to compute singular vectors
! 1032: *
! 1033: IE = 1
! 1034: NRWORK = IE + N
! 1035: ITAUQ = 1
! 1036: ITAUP = ITAUQ + N
! 1037: NWORK = ITAUP + N
! 1038: *
! 1039: * Bidiagonalize A
! 1040: * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
! 1041: * (RWorkspace: need N)
! 1042: *
! 1043: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 1044: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 1045: $ IERR )
! 1046: IF( WNTQN ) THEN
! 1047: *
! 1048: * Compute singular values only
! 1049: * (Cworkspace: 0)
! 1050: * (Rworkspace: need BDSPAN)
! 1051: *
! 1052: CALL DBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
! 1053: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
! 1054: ELSE IF( WNTQO ) THEN
! 1055: IU = NWORK
! 1056: IRU = NRWORK
! 1057: IRVT = IRU + N*N
! 1058: NRWORK = IRVT + N*N
! 1059: IF( LWORK.GE.M*N+3*N ) THEN
! 1060: *
! 1061: * WORK( IU ) is M by N
! 1062: *
! 1063: LDWRKU = M
! 1064: ELSE
! 1065: *
! 1066: * WORK( IU ) is LDWRKU by N
! 1067: *
! 1068: LDWRKU = ( LWORK-3*N ) / N
! 1069: END IF
! 1070: NWORK = IU + LDWRKU*N
! 1071: *
! 1072: * Perform bidiagonal SVD, computing left singular vectors
! 1073: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1074: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1075: * (CWorkspace: need 0)
! 1076: * (RWorkspace: need BDSPAC)
! 1077: *
! 1078: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
! 1079: $ N, RWORK( IRVT ), N, DUM, IDUM,
! 1080: $ RWORK( NRWORK ), IWORK, INFO )
! 1081: *
! 1082: * Copy real matrix RWORK(IRVT) to complex matrix VT
! 1083: * Overwrite VT by right singular vectors of A
! 1084: * (Cworkspace: need 2*N, prefer N+N*NB)
! 1085: * (Rworkspace: need 0)
! 1086: *
! 1087: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
! 1088: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
! 1089: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 1090: $ LWORK-NWORK+1, IERR )
! 1091: *
! 1092: IF( LWORK.GE.M*N+3*N ) THEN
! 1093: *
! 1094: * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
! 1095: * Overwrite WORK(IU) by left singular vectors of A, copying
! 1096: * to A
! 1097: * (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB)
! 1098: * (Rworkspace: need 0)
! 1099: *
! 1100: CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),
! 1101: $ LDWRKU )
! 1102: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
! 1103: $ LDWRKU )
! 1104: CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
! 1105: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
! 1106: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1107: CALL ZLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
! 1108: ELSE
! 1109: *
! 1110: * Generate Q in A
! 1111: * (Cworkspace: need 2*N, prefer N+N*NB)
! 1112: * (Rworkspace: need 0)
! 1113: *
! 1114: CALL ZUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
! 1115: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1116: *
! 1117: * Multiply Q in A by real matrix RWORK(IRU), storing the
! 1118: * result in WORK(IU), copying to A
! 1119: * (CWorkspace: need N*N, prefer M*N)
! 1120: * (Rworkspace: need 3*N*N, prefer N*N+2*M*N)
! 1121: *
! 1122: NRWORK = IRVT
! 1123: DO 30 I = 1, M, LDWRKU
! 1124: CHUNK = MIN( M-I+1, LDWRKU )
! 1125: CALL ZLACRM( CHUNK, N, A( I, 1 ), LDA,
! 1126: $ RWORK( IRU ), N, WORK( IU ), LDWRKU,
! 1127: $ RWORK( NRWORK ) )
! 1128: CALL ZLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
! 1129: $ A( I, 1 ), LDA )
! 1130: 30 CONTINUE
! 1131: END IF
! 1132: *
! 1133: ELSE IF( WNTQS ) THEN
! 1134: *
! 1135: * Perform bidiagonal SVD, computing left singular vectors
! 1136: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1137: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1138: * (CWorkspace: need 0)
! 1139: * (RWorkspace: need BDSPAC)
! 1140: *
! 1141: IRU = NRWORK
! 1142: IRVT = IRU + N*N
! 1143: NRWORK = IRVT + N*N
! 1144: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
! 1145: $ N, RWORK( IRVT ), N, DUM, IDUM,
! 1146: $ RWORK( NRWORK ), IWORK, INFO )
! 1147: *
! 1148: * Copy real matrix RWORK(IRU) to complex matrix U
! 1149: * Overwrite U by left singular vectors of A
! 1150: * (CWorkspace: need 3*N, prefer 2*N+N*NB)
! 1151: * (RWorkspace: 0)
! 1152: *
! 1153: CALL ZLASET( 'F', M, N, CZERO, CZERO, U, LDU )
! 1154: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
! 1155: CALL ZUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
! 1156: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1157: $ LWORK-NWORK+1, IERR )
! 1158: *
! 1159: * Copy real matrix RWORK(IRVT) to complex matrix VT
! 1160: * Overwrite VT by right singular vectors of A
! 1161: * (CWorkspace: need 3*N, prefer 2*N+N*NB)
! 1162: * (RWorkspace: 0)
! 1163: *
! 1164: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
! 1165: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
! 1166: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 1167: $ LWORK-NWORK+1, IERR )
! 1168: ELSE
! 1169: *
! 1170: * Perform bidiagonal SVD, computing left singular vectors
! 1171: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1172: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1173: * (CWorkspace: need 0)
! 1174: * (RWorkspace: need BDSPAC)
! 1175: *
! 1176: IRU = NRWORK
! 1177: IRVT = IRU + N*N
! 1178: NRWORK = IRVT + N*N
! 1179: CALL DBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
! 1180: $ N, RWORK( IRVT ), N, DUM, IDUM,
! 1181: $ RWORK( NRWORK ), IWORK, INFO )
! 1182: *
! 1183: * Set the right corner of U to identity matrix
! 1184: *
! 1185: CALL ZLASET( 'F', M, M, CZERO, CZERO, U, LDU )
! 1186: IF( M.GT.N ) THEN
! 1187: CALL ZLASET( 'F', M-N, M-N, CZERO, CONE,
! 1188: $ U( N+1, N+1 ), LDU )
! 1189: END IF
! 1190: *
! 1191: * Copy real matrix RWORK(IRU) to complex matrix U
! 1192: * Overwrite U by left singular vectors of A
! 1193: * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
! 1194: * (RWorkspace: 0)
! 1195: *
! 1196: CALL ZLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
! 1197: CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
! 1198: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1199: $ LWORK-NWORK+1, IERR )
! 1200: *
! 1201: * Copy real matrix RWORK(IRVT) to complex matrix VT
! 1202: * Overwrite VT by right singular vectors of A
! 1203: * (CWorkspace: need 3*N, prefer 2*N+N*NB)
! 1204: * (RWorkspace: 0)
! 1205: *
! 1206: CALL ZLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
! 1207: CALL ZUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
! 1208: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 1209: $ LWORK-NWORK+1, IERR )
! 1210: END IF
! 1211: *
! 1212: END IF
! 1213: *
! 1214: ELSE
! 1215: *
! 1216: * A has more columns than rows. If A has sufficiently more
! 1217: * columns than rows, first reduce using the LQ decomposition (if
! 1218: * sufficient workspace available)
! 1219: *
! 1220: IF( N.GE.MNTHR1 ) THEN
! 1221: *
! 1222: IF( WNTQN ) THEN
! 1223: *
! 1224: * Path 1t (N much larger than M, JOBZ='N')
! 1225: * No singular vectors to be computed
! 1226: *
! 1227: ITAU = 1
! 1228: NWORK = ITAU + M
! 1229: *
! 1230: * Compute A=L*Q
! 1231: * (CWorkspace: need 2*M, prefer M+M*NB)
! 1232: * (RWorkspace: 0)
! 1233: *
! 1234: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 1235: $ LWORK-NWORK+1, IERR )
! 1236: *
! 1237: * Zero out above L
! 1238: *
! 1239: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
! 1240: $ LDA )
! 1241: IE = 1
! 1242: ITAUQ = 1
! 1243: ITAUP = ITAUQ + M
! 1244: NWORK = ITAUP + M
! 1245: *
! 1246: * Bidiagonalize L in A
! 1247: * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
! 1248: * (RWorkspace: need M)
! 1249: *
! 1250: CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 1251: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 1252: $ IERR )
! 1253: NRWORK = IE + M
! 1254: *
! 1255: * Perform bidiagonal SVD, compute singular values only
! 1256: * (CWorkspace: 0)
! 1257: * (RWorkspace: need BDSPAN)
! 1258: *
! 1259: CALL DBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
! 1260: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
! 1261: *
! 1262: ELSE IF( WNTQO ) THEN
! 1263: *
! 1264: * Path 2t (N much larger than M, JOBZ='O')
! 1265: * M right singular vectors to be overwritten on A and
! 1266: * M left singular vectors to be computed in U
! 1267: *
! 1268: IVT = 1
! 1269: LDWKVT = M
! 1270: *
! 1271: * WORK(IVT) is M by M
! 1272: *
! 1273: IL = IVT + LDWKVT*M
! 1274: IF( LWORK.GE.M*N+M*M+3*M ) THEN
! 1275: *
! 1276: * WORK(IL) M by N
! 1277: *
! 1278: LDWRKL = M
! 1279: CHUNK = N
! 1280: ELSE
! 1281: *
! 1282: * WORK(IL) is M by CHUNK
! 1283: *
! 1284: LDWRKL = M
! 1285: CHUNK = ( LWORK-M*M-3*M ) / M
! 1286: END IF
! 1287: ITAU = IL + LDWRKL*CHUNK
! 1288: NWORK = ITAU + M
! 1289: *
! 1290: * Compute A=L*Q
! 1291: * (CWorkspace: need 2*M, prefer M+M*NB)
! 1292: * (RWorkspace: 0)
! 1293: *
! 1294: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 1295: $ LWORK-NWORK+1, IERR )
! 1296: *
! 1297: * Copy L to WORK(IL), zeroing about above it
! 1298: *
! 1299: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
! 1300: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 1301: $ WORK( IL+LDWRKL ), LDWRKL )
! 1302: *
! 1303: * Generate Q in A
! 1304: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 1305: * (RWorkspace: 0)
! 1306: *
! 1307: CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
! 1308: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1309: IE = 1
! 1310: ITAUQ = ITAU
! 1311: ITAUP = ITAUQ + M
! 1312: NWORK = ITAUP + M
! 1313: *
! 1314: * Bidiagonalize L in WORK(IL)
! 1315: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
! 1316: * (RWorkspace: need M)
! 1317: *
! 1318: CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
! 1319: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
! 1320: $ LWORK-NWORK+1, IERR )
! 1321: *
! 1322: * Perform bidiagonal SVD, computing left singular vectors
! 1323: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1324: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1325: * (CWorkspace: need 0)
! 1326: * (RWorkspace: need BDSPAC)
! 1327: *
! 1328: IRU = IE + M
! 1329: IRVT = IRU + M*M
! 1330: NRWORK = IRVT + M*M
! 1331: CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
! 1332: $ M, RWORK( IRVT ), M, DUM, IDUM,
! 1333: $ RWORK( NRWORK ), IWORK, INFO )
! 1334: *
! 1335: * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
! 1336: * Overwrite WORK(IU) by the left singular vectors of L
! 1337: * (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
! 1338: * (RWorkspace: 0)
! 1339: *
! 1340: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
! 1341: CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
! 1342: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1343: $ LWORK-NWORK+1, IERR )
! 1344: *
! 1345: * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
! 1346: * Overwrite WORK(IVT) by the right singular vectors of L
! 1347: * (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
! 1348: * (RWorkspace: 0)
! 1349: *
! 1350: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
! 1351: $ LDWKVT )
! 1352: CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
! 1353: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
! 1354: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1355: *
! 1356: * Multiply right singular vectors of L in WORK(IL) by Q
! 1357: * in A, storing result in WORK(IL) and copying to A
! 1358: * (CWorkspace: need 2*M*M, prefer M*M+M*N))
! 1359: * (RWorkspace: 0)
! 1360: *
! 1361: DO 40 I = 1, N, CHUNK
! 1362: BLK = MIN( N-I+1, CHUNK )
! 1363: CALL ZGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IVT ), M,
! 1364: $ A( 1, I ), LDA, CZERO, WORK( IL ),
! 1365: $ LDWRKL )
! 1366: CALL ZLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
! 1367: $ A( 1, I ), LDA )
! 1368: 40 CONTINUE
! 1369: *
! 1370: ELSE IF( WNTQS ) THEN
! 1371: *
! 1372: * Path 3t (N much larger than M, JOBZ='S')
! 1373: * M right singular vectors to be computed in VT and
! 1374: * M left singular vectors to be computed in U
! 1375: *
! 1376: IL = 1
! 1377: *
! 1378: * WORK(IL) is M by M
! 1379: *
! 1380: LDWRKL = M
! 1381: ITAU = IL + LDWRKL*M
! 1382: NWORK = ITAU + M
! 1383: *
! 1384: * Compute A=L*Q
! 1385: * (CWorkspace: need 2*M, prefer M+M*NB)
! 1386: * (RWorkspace: 0)
! 1387: *
! 1388: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 1389: $ LWORK-NWORK+1, IERR )
! 1390: *
! 1391: * Copy L to WORK(IL), zeroing out above it
! 1392: *
! 1393: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
! 1394: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO,
! 1395: $ WORK( IL+LDWRKL ), LDWRKL )
! 1396: *
! 1397: * Generate Q in A
! 1398: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
! 1399: * (RWorkspace: 0)
! 1400: *
! 1401: CALL ZUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
! 1402: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1403: IE = 1
! 1404: ITAUQ = ITAU
! 1405: ITAUP = ITAUQ + M
! 1406: NWORK = ITAUP + M
! 1407: *
! 1408: * Bidiagonalize L in WORK(IL)
! 1409: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
! 1410: * (RWorkspace: need M)
! 1411: *
! 1412: CALL ZGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
! 1413: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
! 1414: $ LWORK-NWORK+1, IERR )
! 1415: *
! 1416: * Perform bidiagonal SVD, computing left singular vectors
! 1417: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1418: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1419: * (CWorkspace: need 0)
! 1420: * (RWorkspace: need BDSPAC)
! 1421: *
! 1422: IRU = IE + M
! 1423: IRVT = IRU + M*M
! 1424: NRWORK = IRVT + M*M
! 1425: CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
! 1426: $ M, RWORK( IRVT ), M, DUM, IDUM,
! 1427: $ RWORK( NRWORK ), IWORK, INFO )
! 1428: *
! 1429: * Copy real matrix RWORK(IRU) to complex matrix U
! 1430: * Overwrite U by left singular vectors of L
! 1431: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
! 1432: * (RWorkspace: 0)
! 1433: *
! 1434: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
! 1435: CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
! 1436: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1437: $ LWORK-NWORK+1, IERR )
! 1438: *
! 1439: * Copy real matrix RWORK(IRVT) to complex matrix VT
! 1440: * Overwrite VT by left singular vectors of L
! 1441: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
! 1442: * (RWorkspace: 0)
! 1443: *
! 1444: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
! 1445: CALL ZUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
! 1446: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 1447: $ LWORK-NWORK+1, IERR )
! 1448: *
! 1449: * Copy VT to WORK(IL), multiply right singular vectors of L
! 1450: * in WORK(IL) by Q in A, storing result in VT
! 1451: * (CWorkspace: need M*M)
! 1452: * (RWorkspace: 0)
! 1453: *
! 1454: CALL ZLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
! 1455: CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,
! 1456: $ A, LDA, CZERO, VT, LDVT )
! 1457: *
! 1458: ELSE IF( WNTQA ) THEN
! 1459: *
! 1460: * Path 9t (N much larger than M, JOBZ='A')
! 1461: * N right singular vectors to be computed in VT and
! 1462: * M left singular vectors to be computed in U
! 1463: *
! 1464: IVT = 1
! 1465: *
! 1466: * WORK(IVT) is M by M
! 1467: *
! 1468: LDWKVT = M
! 1469: ITAU = IVT + LDWKVT*M
! 1470: NWORK = ITAU + M
! 1471: *
! 1472: * Compute A=L*Q, copying result to VT
! 1473: * (CWorkspace: need 2*M, prefer M+M*NB)
! 1474: * (RWorkspace: 0)
! 1475: *
! 1476: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 1477: $ LWORK-NWORK+1, IERR )
! 1478: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 1479: *
! 1480: * Generate Q in VT
! 1481: * (CWorkspace: need M+N, prefer M+N*NB)
! 1482: * (RWorkspace: 0)
! 1483: *
! 1484: CALL ZUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
! 1485: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1486: *
! 1487: * Produce L in A, zeroing out above it
! 1488: *
! 1489: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
! 1490: $ LDA )
! 1491: IE = 1
! 1492: ITAUQ = ITAU
! 1493: ITAUP = ITAUQ + M
! 1494: NWORK = ITAUP + M
! 1495: *
! 1496: * Bidiagonalize L in A
! 1497: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
! 1498: * (RWorkspace: need M)
! 1499: *
! 1500: CALL ZGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 1501: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 1502: $ IERR )
! 1503: *
! 1504: * Perform bidiagonal SVD, computing left singular vectors
! 1505: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1506: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1507: * (CWorkspace: need 0)
! 1508: * (RWorkspace: need BDSPAC)
! 1509: *
! 1510: IRU = IE + M
! 1511: IRVT = IRU + M*M
! 1512: NRWORK = IRVT + M*M
! 1513: CALL DBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
! 1514: $ M, RWORK( IRVT ), M, DUM, IDUM,
! 1515: $ RWORK( NRWORK ), IWORK, INFO )
! 1516: *
! 1517: * Copy real matrix RWORK(IRU) to complex matrix U
! 1518: * Overwrite U by left singular vectors of L
! 1519: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 1520: * (RWorkspace: 0)
! 1521: *
! 1522: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
! 1523: CALL ZUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
! 1524: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1525: $ LWORK-NWORK+1, IERR )
! 1526: *
! 1527: * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
! 1528: * Overwrite WORK(IVT) by right singular vectors of L
! 1529: * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
! 1530: * (RWorkspace: 0)
! 1531: *
! 1532: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
! 1533: $ LDWKVT )
! 1534: CALL ZUNMBR( 'P', 'R', 'C', M, M, M, A, LDA,
! 1535: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
! 1536: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1537: *
! 1538: * Multiply right singular vectors of L in WORK(IVT) by
! 1539: * Q in VT, storing result in A
! 1540: * (CWorkspace: need M*M)
! 1541: * (RWorkspace: 0)
! 1542: *
! 1543: CALL ZGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT,
! 1544: $ VT, LDVT, CZERO, A, LDA )
! 1545: *
! 1546: * Copy right singular vectors of A from A to VT
! 1547: *
! 1548: CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
! 1549: *
! 1550: END IF
! 1551: *
! 1552: ELSE IF( N.GE.MNTHR2 ) THEN
! 1553: *
! 1554: * MNTHR2 <= N < MNTHR1
! 1555: *
! 1556: * Path 5t (N much larger than M, but not as much as MNTHR1)
! 1557: * Reduce to bidiagonal form without QR decomposition, use
! 1558: * ZUNGBR and matrix multiplication to compute singular vectors
! 1559: *
! 1560: *
! 1561: IE = 1
! 1562: NRWORK = IE + M
! 1563: ITAUQ = 1
! 1564: ITAUP = ITAUQ + M
! 1565: NWORK = ITAUP + M
! 1566: *
! 1567: * Bidiagonalize A
! 1568: * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
! 1569: * (RWorkspace: M)
! 1570: *
! 1571: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 1572: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 1573: $ IERR )
! 1574: *
! 1575: IF( WNTQN ) THEN
! 1576: *
! 1577: * Compute singular values only
! 1578: * (Cworkspace: 0)
! 1579: * (Rworkspace: need BDSPAN)
! 1580: *
! 1581: CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
! 1582: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
! 1583: ELSE IF( WNTQO ) THEN
! 1584: IRVT = NRWORK
! 1585: IRU = IRVT + M*M
! 1586: NRWORK = IRU + M*M
! 1587: IVT = NWORK
! 1588: *
! 1589: * Copy A to U, generate Q
! 1590: * (Cworkspace: need 2*M, prefer M+M*NB)
! 1591: * (Rworkspace: 0)
! 1592: *
! 1593: CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
! 1594: CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
! 1595: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1596: *
! 1597: * Generate P**H in A
! 1598: * (Cworkspace: need 2*M, prefer M+M*NB)
! 1599: * (Rworkspace: 0)
! 1600: *
! 1601: CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
! 1602: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1603: *
! 1604: LDWKVT = M
! 1605: IF( LWORK.GE.M*N+3*M ) THEN
! 1606: *
! 1607: * WORK( IVT ) is M by N
! 1608: *
! 1609: NWORK = IVT + LDWKVT*N
! 1610: CHUNK = N
! 1611: ELSE
! 1612: *
! 1613: * WORK( IVT ) is M by CHUNK
! 1614: *
! 1615: CHUNK = ( LWORK-3*M ) / M
! 1616: NWORK = IVT + LDWKVT*CHUNK
! 1617: END IF
! 1618: *
! 1619: * Perform bidiagonal SVD, computing left singular vectors
! 1620: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1621: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1622: * (CWorkspace: need 0)
! 1623: * (RWorkspace: need BDSPAC)
! 1624: *
! 1625: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
! 1626: $ M, RWORK( IRVT ), M, DUM, IDUM,
! 1627: $ RWORK( NRWORK ), IWORK, INFO )
! 1628: *
! 1629: * Multiply Q in U by real matrix RWORK(IRVT)
! 1630: * storing the result in WORK(IVT), copying to U
! 1631: * (Cworkspace: need 0)
! 1632: * (Rworkspace: need 2*M*M)
! 1633: *
! 1634: CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ),
! 1635: $ LDWKVT, RWORK( NRWORK ) )
! 1636: CALL ZLACPY( 'F', M, M, WORK( IVT ), LDWKVT, U, LDU )
! 1637: *
! 1638: * Multiply RWORK(IRVT) by P**H in A, storing the
! 1639: * result in WORK(IVT), copying to A
! 1640: * (CWorkspace: need M*M, prefer M*N)
! 1641: * (Rworkspace: need 2*M*M, prefer 2*M*N)
! 1642: *
! 1643: NRWORK = IRU
! 1644: DO 50 I = 1, N, CHUNK
! 1645: BLK = MIN( N-I+1, CHUNK )
! 1646: CALL ZLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), LDA,
! 1647: $ WORK( IVT ), LDWKVT, RWORK( NRWORK ) )
! 1648: CALL ZLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
! 1649: $ A( 1, I ), LDA )
! 1650: 50 CONTINUE
! 1651: ELSE IF( WNTQS ) THEN
! 1652: *
! 1653: * Copy A to U, generate Q
! 1654: * (Cworkspace: need 2*M, prefer M+M*NB)
! 1655: * (Rworkspace: 0)
! 1656: *
! 1657: CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
! 1658: CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
! 1659: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1660: *
! 1661: * Copy A to VT, generate P**H
! 1662: * (Cworkspace: need 2*M, prefer M+M*NB)
! 1663: * (Rworkspace: 0)
! 1664: *
! 1665: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 1666: CALL ZUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ),
! 1667: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1668: *
! 1669: * Perform bidiagonal SVD, computing left singular vectors
! 1670: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1671: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1672: * (CWorkspace: need 0)
! 1673: * (RWorkspace: need BDSPAC)
! 1674: *
! 1675: IRVT = NRWORK
! 1676: IRU = IRVT + M*M
! 1677: NRWORK = IRU + M*M
! 1678: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
! 1679: $ M, RWORK( IRVT ), M, DUM, IDUM,
! 1680: $ RWORK( NRWORK ), IWORK, INFO )
! 1681: *
! 1682: * Multiply Q in U by real matrix RWORK(IRU), storing the
! 1683: * result in A, copying to U
! 1684: * (CWorkspace: need 0)
! 1685: * (Rworkspace: need 3*M*M)
! 1686: *
! 1687: CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
! 1688: $ RWORK( NRWORK ) )
! 1689: CALL ZLACPY( 'F', M, M, A, LDA, U, LDU )
! 1690: *
! 1691: * Multiply real matrix RWORK(IRVT) by P**H in VT,
! 1692: * storing the result in A, copying to VT
! 1693: * (Cworkspace: need 0)
! 1694: * (Rworkspace: need M*M+2*M*N)
! 1695: *
! 1696: NRWORK = IRU
! 1697: CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
! 1698: $ RWORK( NRWORK ) )
! 1699: CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
! 1700: ELSE
! 1701: *
! 1702: * Copy A to U, generate Q
! 1703: * (Cworkspace: need 2*M, prefer M+M*NB)
! 1704: * (Rworkspace: 0)
! 1705: *
! 1706: CALL ZLACPY( 'L', M, M, A, LDA, U, LDU )
! 1707: CALL ZUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
! 1708: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1709: *
! 1710: * Copy A to VT, generate P**H
! 1711: * (Cworkspace: need 2*M, prefer M+M*NB)
! 1712: * (Rworkspace: 0)
! 1713: *
! 1714: CALL ZLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 1715: CALL ZUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ),
! 1716: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1717: *
! 1718: * Perform bidiagonal SVD, computing left singular vectors
! 1719: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1720: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1721: * (CWorkspace: need 0)
! 1722: * (RWorkspace: need BDSPAC)
! 1723: *
! 1724: IRVT = NRWORK
! 1725: IRU = IRVT + M*M
! 1726: NRWORK = IRU + M*M
! 1727: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
! 1728: $ M, RWORK( IRVT ), M, DUM, IDUM,
! 1729: $ RWORK( NRWORK ), IWORK, INFO )
! 1730: *
! 1731: * Multiply Q in U by real matrix RWORK(IRU), storing the
! 1732: * result in A, copying to U
! 1733: * (CWorkspace: need 0)
! 1734: * (Rworkspace: need 3*M*M)
! 1735: *
! 1736: CALL ZLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
! 1737: $ RWORK( NRWORK ) )
! 1738: CALL ZLACPY( 'F', M, M, A, LDA, U, LDU )
! 1739: *
! 1740: * Multiply real matrix RWORK(IRVT) by P**H in VT,
! 1741: * storing the result in A, copying to VT
! 1742: * (Cworkspace: need 0)
! 1743: * (Rworkspace: need M*M+2*M*N)
! 1744: *
! 1745: CALL ZLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
! 1746: $ RWORK( NRWORK ) )
! 1747: CALL ZLACPY( 'F', M, N, A, LDA, VT, LDVT )
! 1748: END IF
! 1749: *
! 1750: ELSE
! 1751: *
! 1752: * N .LT. MNTHR2
! 1753: *
! 1754: * Path 6t (N greater than M, but not much larger)
! 1755: * Reduce to bidiagonal form without LQ decomposition
! 1756: * Use ZUNMBR to compute singular vectors
! 1757: *
! 1758: IE = 1
! 1759: NRWORK = IE + M
! 1760: ITAUQ = 1
! 1761: ITAUP = ITAUQ + M
! 1762: NWORK = ITAUP + M
! 1763: *
! 1764: * Bidiagonalize A
! 1765: * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
! 1766: * (RWorkspace: M)
! 1767: *
! 1768: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 1769: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 1770: $ IERR )
! 1771: IF( WNTQN ) THEN
! 1772: *
! 1773: * Compute singular values only
! 1774: * (Cworkspace: 0)
! 1775: * (Rworkspace: need BDSPAN)
! 1776: *
! 1777: CALL DBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
! 1778: $ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
! 1779: ELSE IF( WNTQO ) THEN
! 1780: LDWKVT = M
! 1781: IVT = NWORK
! 1782: IF( LWORK.GE.M*N+3*M ) THEN
! 1783: *
! 1784: * WORK( IVT ) is M by N
! 1785: *
! 1786: CALL ZLASET( 'F', M, N, CZERO, CZERO, WORK( IVT ),
! 1787: $ LDWKVT )
! 1788: NWORK = IVT + LDWKVT*N
! 1789: ELSE
! 1790: *
! 1791: * WORK( IVT ) is M by CHUNK
! 1792: *
! 1793: CHUNK = ( LWORK-3*M ) / M
! 1794: NWORK = IVT + LDWKVT*CHUNK
! 1795: END IF
! 1796: *
! 1797: * Perform bidiagonal SVD, computing left singular vectors
! 1798: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1799: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1800: * (CWorkspace: need 0)
! 1801: * (RWorkspace: need BDSPAC)
! 1802: *
! 1803: IRVT = NRWORK
! 1804: IRU = IRVT + M*M
! 1805: NRWORK = IRU + M*M
! 1806: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
! 1807: $ M, RWORK( IRVT ), M, DUM, IDUM,
! 1808: $ RWORK( NRWORK ), IWORK, INFO )
! 1809: *
! 1810: * Copy real matrix RWORK(IRU) to complex matrix U
! 1811: * Overwrite U by left singular vectors of A
! 1812: * (Cworkspace: need 2*M, prefer M+M*NB)
! 1813: * (Rworkspace: need 0)
! 1814: *
! 1815: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
! 1816: CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
! 1817: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1818: $ LWORK-NWORK+1, IERR )
! 1819: *
! 1820: IF( LWORK.GE.M*N+3*M ) THEN
! 1821: *
! 1822: * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
! 1823: * Overwrite WORK(IVT) by right singular vectors of A,
! 1824: * copying to A
! 1825: * (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB)
! 1826: * (Rworkspace: need 0)
! 1827: *
! 1828: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
! 1829: $ LDWKVT )
! 1830: CALL ZUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
! 1831: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
! 1832: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1833: CALL ZLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
! 1834: ELSE
! 1835: *
! 1836: * Generate P**H in A
! 1837: * (Cworkspace: need 2*M, prefer M+M*NB)
! 1838: * (Rworkspace: need 0)
! 1839: *
! 1840: CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
! 1841: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1842: *
! 1843: * Multiply Q in A by real matrix RWORK(IRU), storing the
! 1844: * result in WORK(IU), copying to A
! 1845: * (CWorkspace: need M*M, prefer M*N)
! 1846: * (Rworkspace: need 3*M*M, prefer M*M+2*M*N)
! 1847: *
! 1848: NRWORK = IRU
! 1849: DO 60 I = 1, N, CHUNK
! 1850: BLK = MIN( N-I+1, CHUNK )
! 1851: CALL ZLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ),
! 1852: $ LDA, WORK( IVT ), LDWKVT,
! 1853: $ RWORK( NRWORK ) )
! 1854: CALL ZLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
! 1855: $ A( 1, I ), LDA )
! 1856: 60 CONTINUE
! 1857: END IF
! 1858: ELSE IF( WNTQS ) THEN
! 1859: *
! 1860: * Perform bidiagonal SVD, computing left singular vectors
! 1861: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1862: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1863: * (CWorkspace: need 0)
! 1864: * (RWorkspace: need BDSPAC)
! 1865: *
! 1866: IRVT = NRWORK
! 1867: IRU = IRVT + M*M
! 1868: NRWORK = IRU + M*M
! 1869: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
! 1870: $ M, RWORK( IRVT ), M, DUM, IDUM,
! 1871: $ RWORK( NRWORK ), IWORK, INFO )
! 1872: *
! 1873: * Copy real matrix RWORK(IRU) to complex matrix U
! 1874: * Overwrite U by left singular vectors of A
! 1875: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 1876: * (RWorkspace: M*M)
! 1877: *
! 1878: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
! 1879: CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
! 1880: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1881: $ LWORK-NWORK+1, IERR )
! 1882: *
! 1883: * Copy real matrix RWORK(IRVT) to complex matrix VT
! 1884: * Overwrite VT by right singular vectors of A
! 1885: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 1886: * (RWorkspace: M*M)
! 1887: *
! 1888: CALL ZLASET( 'F', M, N, CZERO, CZERO, VT, LDVT )
! 1889: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
! 1890: CALL ZUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
! 1891: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 1892: $ LWORK-NWORK+1, IERR )
! 1893: ELSE
! 1894: *
! 1895: * Perform bidiagonal SVD, computing left singular vectors
! 1896: * of bidiagonal matrix in RWORK(IRU) and computing right
! 1897: * singular vectors of bidiagonal matrix in RWORK(IRVT)
! 1898: * (CWorkspace: need 0)
! 1899: * (RWorkspace: need BDSPAC)
! 1900: *
! 1901: IRVT = NRWORK
! 1902: IRU = IRVT + M*M
! 1903: NRWORK = IRU + M*M
! 1904: *
! 1905: CALL DBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
! 1906: $ M, RWORK( IRVT ), M, DUM, IDUM,
! 1907: $ RWORK( NRWORK ), IWORK, INFO )
! 1908: *
! 1909: * Copy real matrix RWORK(IRU) to complex matrix U
! 1910: * Overwrite U by left singular vectors of A
! 1911: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 1912: * (RWorkspace: M*M)
! 1913: *
! 1914: CALL ZLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
! 1915: CALL ZUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
! 1916: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1917: $ LWORK-NWORK+1, IERR )
! 1918: *
! 1919: * Set all of VT to identity matrix
! 1920: *
! 1921: CALL ZLASET( 'F', N, N, CZERO, CONE, VT, LDVT )
! 1922: *
! 1923: * Copy real matrix RWORK(IRVT) to complex matrix VT
! 1924: * Overwrite VT by right singular vectors of A
! 1925: * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
! 1926: * (RWorkspace: M*M)
! 1927: *
! 1928: CALL ZLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
! 1929: CALL ZUNMBR( 'P', 'R', 'C', N, N, M, A, LDA,
! 1930: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 1931: $ LWORK-NWORK+1, IERR )
! 1932: END IF
! 1933: *
! 1934: END IF
! 1935: *
! 1936: END IF
! 1937: *
! 1938: * Undo scaling if necessary
! 1939: *
! 1940: IF( ISCL.EQ.1 ) THEN
! 1941: IF( ANRM.GT.BIGNUM )
! 1942: $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
! 1943: $ IERR )
! 1944: IF( INFO.NE.0 .AND. ANRM.GT.BIGNUM )
! 1945: $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN-1, 1,
! 1946: $ RWORK( IE ), MINMN, IERR )
! 1947: IF( ANRM.LT.SMLNUM )
! 1948: $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
! 1949: $ IERR )
! 1950: IF( INFO.NE.0 .AND. ANRM.LT.SMLNUM )
! 1951: $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN-1, 1,
! 1952: $ RWORK( IE ), MINMN, IERR )
! 1953: END IF
! 1954: *
! 1955: * Return optimal workspace in WORK(1)
! 1956: *
! 1957: WORK( 1 ) = MAXWRK
! 1958: *
! 1959: RETURN
! 1960: *
! 1961: * End of ZGESDD
! 1962: *
! 1963: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgesdd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack