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