Annotation of rpl/lapack/lapack/dgesdd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
! 2: $ LWORK, IWORK, INFO )
! 3: *
! 4: * -- LAPACK driver routine (version 3.2.1) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * March 2009
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER JOBZ
! 11: INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
! 12: * ..
! 13: * .. Array Arguments ..
! 14: INTEGER IWORK( * )
! 15: DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
! 16: $ VT( LDVT, * ), WORK( * )
! 17: * ..
! 18: *
! 19: * Purpose
! 20: * =======
! 21: *
! 22: * DGESDD computes the singular value decomposition (SVD) of a real
! 23: * M-by-N matrix A, optionally computing the left and right singular
! 24: * vectors. If singular vectors are desired, it uses a
! 25: * divide-and-conquer algorithm.
! 26: *
! 27: * The SVD is written
! 28: *
! 29: * A = U * SIGMA * transpose(V)
! 30: *
! 31: * where SIGMA is an M-by-N matrix which is zero except for its
! 32: * min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
! 33: * V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
! 34: * are the singular values of A; they are real and non-negative, and
! 35: * are returned in descending order. The first min(m,n) columns of
! 36: * U and V are the left and right singular vectors of A.
! 37: *
! 38: * Note that the routine returns VT = V**T, not V.
! 39: *
! 40: * The divide and conquer algorithm makes very mild assumptions about
! 41: * floating point arithmetic. It will work on machines with a guard
! 42: * digit in add/subtract, or on those binary machines without guard
! 43: * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
! 44: * Cray-2. It could conceivably fail on hexadecimal or decimal machines
! 45: * without guard digits, but we know of none.
! 46: *
! 47: * Arguments
! 48: * =========
! 49: *
! 50: * JOBZ (input) CHARACTER*1
! 51: * Specifies options for computing all or part of the matrix U:
! 52: * = 'A': all M columns of U and all N rows of V**T are
! 53: * returned in the arrays U and VT;
! 54: * = 'S': the first min(M,N) columns of U and the first
! 55: * min(M,N) rows of V**T are returned in the arrays U
! 56: * and VT;
! 57: * = 'O': If M >= N, the first N columns of U are overwritten
! 58: * on the array A and all rows of V**T are returned in
! 59: * the array VT;
! 60: * otherwise, all columns of U are returned in the
! 61: * array U and the first M rows of V**T are overwritten
! 62: * in the array A;
! 63: * = 'N': no columns of U or rows of V**T are computed.
! 64: *
! 65: * M (input) INTEGER
! 66: * The number of rows of the input matrix A. M >= 0.
! 67: *
! 68: * N (input) INTEGER
! 69: * The number of columns of the input matrix A. N >= 0.
! 70: *
! 71: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
! 72: * On entry, the M-by-N matrix A.
! 73: * On exit,
! 74: * if JOBZ = 'O', A is overwritten with the first N columns
! 75: * of U (the left singular vectors, stored
! 76: * columnwise) if M >= N;
! 77: * A is overwritten with the first M rows
! 78: * of V**T (the right singular vectors, stored
! 79: * rowwise) otherwise.
! 80: * if JOBZ .ne. 'O', the contents of A are destroyed.
! 81: *
! 82: * LDA (input) INTEGER
! 83: * The leading dimension of the array A. LDA >= max(1,M).
! 84: *
! 85: * S (output) DOUBLE PRECISION array, dimension (min(M,N))
! 86: * The singular values of A, sorted so that S(i) >= S(i+1).
! 87: *
! 88: * U (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
! 89: * UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
! 90: * UCOL = min(M,N) if JOBZ = 'S'.
! 91: * If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
! 92: * orthogonal matrix U;
! 93: * if JOBZ = 'S', U contains the first min(M,N) columns of U
! 94: * (the left singular vectors, stored columnwise);
! 95: * if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
! 96: *
! 97: * LDU (input) INTEGER
! 98: * The leading dimension of the array U. LDU >= 1; if
! 99: * JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
! 100: *
! 101: * VT (output) DOUBLE PRECISION array, dimension (LDVT,N)
! 102: * If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
! 103: * N-by-N orthogonal matrix V**T;
! 104: * if JOBZ = 'S', VT contains the first min(M,N) rows of
! 105: * V**T (the right singular vectors, stored rowwise);
! 106: * if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
! 107: *
! 108: * LDVT (input) INTEGER
! 109: * The leading dimension of the array VT. LDVT >= 1; if
! 110: * JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
! 111: * if JOBZ = 'S', LDVT >= min(M,N).
! 112: *
! 113: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 114: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
! 115: *
! 116: * LWORK (input) INTEGER
! 117: * The dimension of the array WORK. LWORK >= 1.
! 118: * If JOBZ = 'N',
! 119: * LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)).
! 120: * If JOBZ = 'O',
! 121: * LWORK >= 3*min(M,N) +
! 122: * max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
! 123: * If JOBZ = 'S' or 'A'
! 124: * LWORK >= 3*min(M,N) +
! 125: * max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).
! 126: * For good performance, LWORK should generally be larger.
! 127: * If LWORK = -1 but other input arguments are legal, WORK(1)
! 128: * returns the optimal LWORK.
! 129: *
! 130: * IWORK (workspace) INTEGER array, dimension (8*min(M,N))
! 131: *
! 132: * INFO (output) INTEGER
! 133: * = 0: successful exit.
! 134: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 135: * > 0: DBDSDC did not converge, updating process failed.
! 136: *
! 137: * Further Details
! 138: * ===============
! 139: *
! 140: * Based on contributions by
! 141: * Ming Gu and Huan Ren, Computer Science Division, University of
! 142: * California at Berkeley, USA
! 143: *
! 144: * =====================================================================
! 145: *
! 146: * .. Parameters ..
! 147: DOUBLE PRECISION ZERO, ONE
! 148: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 149: * ..
! 150: * .. Local Scalars ..
! 151: LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
! 152: INTEGER BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
! 153: $ IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
! 154: $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
! 155: $ MNTHR, NWORK, WRKBL
! 156: DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
! 157: * ..
! 158: * .. Local Arrays ..
! 159: INTEGER IDUM( 1 )
! 160: DOUBLE PRECISION DUM( 1 )
! 161: * ..
! 162: * .. External Subroutines ..
! 163: EXTERNAL DBDSDC, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY,
! 164: $ DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR,
! 165: $ XERBLA
! 166: * ..
! 167: * .. External Functions ..
! 168: LOGICAL LSAME
! 169: INTEGER ILAENV
! 170: DOUBLE PRECISION DLAMCH, DLANGE
! 171: EXTERNAL DLAMCH, DLANGE, ILAENV, LSAME
! 172: * ..
! 173: * .. Intrinsic Functions ..
! 174: INTRINSIC INT, MAX, MIN, SQRT
! 175: * ..
! 176: * .. Executable Statements ..
! 177: *
! 178: * Test the input arguments
! 179: *
! 180: INFO = 0
! 181: MINMN = MIN( M, N )
! 182: WNTQA = LSAME( JOBZ, 'A' )
! 183: WNTQS = LSAME( JOBZ, 'S' )
! 184: WNTQAS = WNTQA .OR. WNTQS
! 185: WNTQO = LSAME( JOBZ, 'O' )
! 186: WNTQN = LSAME( JOBZ, 'N' )
! 187: LQUERY = ( LWORK.EQ.-1 )
! 188: *
! 189: IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
! 190: INFO = -1
! 191: ELSE IF( M.LT.0 ) THEN
! 192: INFO = -2
! 193: ELSE IF( N.LT.0 ) THEN
! 194: INFO = -3
! 195: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 196: INFO = -5
! 197: ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
! 198: $ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
! 199: INFO = -8
! 200: ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
! 201: $ ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
! 202: $ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
! 203: INFO = -10
! 204: END IF
! 205: *
! 206: * Compute workspace
! 207: * (Note: Comments in the code beginning "Workspace:" describe the
! 208: * minimal amount of workspace needed at that point in the code,
! 209: * as well as the preferred amount for good performance.
! 210: * NB refers to the optimal block size for the immediately
! 211: * following subroutine, as returned by ILAENV.)
! 212: *
! 213: IF( INFO.EQ.0 ) THEN
! 214: MINWRK = 1
! 215: MAXWRK = 1
! 216: IF( M.GE.N .AND. MINMN.GT.0 ) THEN
! 217: *
! 218: * Compute space needed for DBDSDC
! 219: *
! 220: MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
! 221: IF( WNTQN ) THEN
! 222: BDSPAC = 7*N
! 223: ELSE
! 224: BDSPAC = 3*N*N + 4*N
! 225: END IF
! 226: IF( M.GE.MNTHR ) THEN
! 227: IF( WNTQN ) THEN
! 228: *
! 229: * Path 1 (M much larger than N, JOBZ='N')
! 230: *
! 231: WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1,
! 232: $ -1 )
! 233: WRKBL = MAX( WRKBL, 3*N+2*N*
! 234: $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
! 235: MAXWRK = MAX( WRKBL, BDSPAC+N )
! 236: MINWRK = BDSPAC + N
! 237: ELSE IF( WNTQO ) THEN
! 238: *
! 239: * Path 2 (M much larger than N, JOBZ='O')
! 240: *
! 241: WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
! 242: WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
! 243: $ N, N, -1 ) )
! 244: WRKBL = MAX( WRKBL, 3*N+2*N*
! 245: $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
! 246: WRKBL = MAX( WRKBL, 3*N+N*
! 247: $ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
! 248: WRKBL = MAX( WRKBL, 3*N+N*
! 249: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
! 250: WRKBL = MAX( WRKBL, BDSPAC+3*N )
! 251: MAXWRK = WRKBL + 2*N*N
! 252: MINWRK = BDSPAC + 2*N*N + 3*N
! 253: ELSE IF( WNTQS ) THEN
! 254: *
! 255: * Path 3 (M much larger than N, JOBZ='S')
! 256: *
! 257: WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
! 258: WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
! 259: $ N, N, -1 ) )
! 260: WRKBL = MAX( WRKBL, 3*N+2*N*
! 261: $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
! 262: WRKBL = MAX( WRKBL, 3*N+N*
! 263: $ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
! 264: WRKBL = MAX( WRKBL, 3*N+N*
! 265: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
! 266: WRKBL = MAX( WRKBL, BDSPAC+3*N )
! 267: MAXWRK = WRKBL + N*N
! 268: MINWRK = BDSPAC + N*N + 3*N
! 269: ELSE IF( WNTQA ) THEN
! 270: *
! 271: * Path 4 (M much larger than N, JOBZ='A')
! 272: *
! 273: WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
! 274: WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,
! 275: $ M, N, -1 ) )
! 276: WRKBL = MAX( WRKBL, 3*N+2*N*
! 277: $ ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
! 278: WRKBL = MAX( WRKBL, 3*N+N*
! 279: $ ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
! 280: WRKBL = MAX( WRKBL, 3*N+N*
! 281: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
! 282: WRKBL = MAX( WRKBL, BDSPAC+3*N )
! 283: MAXWRK = WRKBL + N*N
! 284: MINWRK = BDSPAC + N*N + 3*N
! 285: END IF
! 286: ELSE
! 287: *
! 288: * Path 5 (M at least N, but not much larger)
! 289: *
! 290: WRKBL = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
! 291: $ -1 )
! 292: IF( WNTQN ) THEN
! 293: MAXWRK = MAX( WRKBL, BDSPAC+3*N )
! 294: MINWRK = 3*N + MAX( M, BDSPAC )
! 295: ELSE IF( WNTQO ) THEN
! 296: WRKBL = MAX( WRKBL, 3*N+N*
! 297: $ ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
! 298: WRKBL = MAX( WRKBL, 3*N+N*
! 299: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
! 300: WRKBL = MAX( WRKBL, BDSPAC+3*N )
! 301: MAXWRK = WRKBL + M*N
! 302: MINWRK = 3*N + MAX( M, N*N+BDSPAC )
! 303: ELSE IF( WNTQS ) THEN
! 304: WRKBL = MAX( WRKBL, 3*N+N*
! 305: $ ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
! 306: WRKBL = MAX( WRKBL, 3*N+N*
! 307: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
! 308: MAXWRK = MAX( WRKBL, BDSPAC+3*N )
! 309: MINWRK = 3*N + MAX( M, BDSPAC )
! 310: ELSE IF( WNTQA ) THEN
! 311: WRKBL = MAX( WRKBL, 3*N+M*
! 312: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
! 313: WRKBL = MAX( WRKBL, 3*N+N*
! 314: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
! 315: MAXWRK = MAX( MAXWRK, BDSPAC+3*N )
! 316: MINWRK = 3*N + MAX( M, BDSPAC )
! 317: END IF
! 318: END IF
! 319: ELSE IF( MINMN.GT.0 ) THEN
! 320: *
! 321: * Compute space needed for DBDSDC
! 322: *
! 323: MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
! 324: IF( WNTQN ) THEN
! 325: BDSPAC = 7*M
! 326: ELSE
! 327: BDSPAC = 3*M*M + 4*M
! 328: END IF
! 329: IF( N.GE.MNTHR ) THEN
! 330: IF( WNTQN ) THEN
! 331: *
! 332: * Path 1t (N much larger than M, JOBZ='N')
! 333: *
! 334: WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1,
! 335: $ -1 )
! 336: WRKBL = MAX( WRKBL, 3*M+2*M*
! 337: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
! 338: MAXWRK = MAX( WRKBL, BDSPAC+M )
! 339: MINWRK = BDSPAC + M
! 340: ELSE IF( WNTQO ) THEN
! 341: *
! 342: * Path 2t (N much larger than M, JOBZ='O')
! 343: *
! 344: WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
! 345: WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
! 346: $ N, M, -1 ) )
! 347: WRKBL = MAX( WRKBL, 3*M+2*M*
! 348: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
! 349: WRKBL = MAX( WRKBL, 3*M+M*
! 350: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
! 351: WRKBL = MAX( WRKBL, 3*M+M*
! 352: $ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
! 353: WRKBL = MAX( WRKBL, BDSPAC+3*M )
! 354: MAXWRK = WRKBL + 2*M*M
! 355: MINWRK = BDSPAC + 2*M*M + 3*M
! 356: ELSE IF( WNTQS ) THEN
! 357: *
! 358: * Path 3t (N much larger than M, JOBZ='S')
! 359: *
! 360: WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
! 361: WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
! 362: $ N, M, -1 ) )
! 363: WRKBL = MAX( WRKBL, 3*M+2*M*
! 364: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
! 365: WRKBL = MAX( WRKBL, 3*M+M*
! 366: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
! 367: WRKBL = MAX( WRKBL, 3*M+M*
! 368: $ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
! 369: WRKBL = MAX( WRKBL, BDSPAC+3*M )
! 370: MAXWRK = WRKBL + M*M
! 371: MINWRK = BDSPAC + M*M + 3*M
! 372: ELSE IF( WNTQA ) THEN
! 373: *
! 374: * Path 4t (N much larger than M, JOBZ='A')
! 375: *
! 376: WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
! 377: WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,
! 378: $ N, M, -1 ) )
! 379: WRKBL = MAX( WRKBL, 3*M+2*M*
! 380: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
! 381: WRKBL = MAX( WRKBL, 3*M+M*
! 382: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
! 383: WRKBL = MAX( WRKBL, 3*M+M*
! 384: $ ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
! 385: WRKBL = MAX( WRKBL, BDSPAC+3*M )
! 386: MAXWRK = WRKBL + M*M
! 387: MINWRK = BDSPAC + M*M + 3*M
! 388: END IF
! 389: ELSE
! 390: *
! 391: * Path 5t (N greater than M, but not much larger)
! 392: *
! 393: WRKBL = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
! 394: $ -1 )
! 395: IF( WNTQN ) THEN
! 396: MAXWRK = MAX( WRKBL, BDSPAC+3*M )
! 397: MINWRK = 3*M + MAX( N, BDSPAC )
! 398: ELSE IF( WNTQO ) THEN
! 399: WRKBL = MAX( WRKBL, 3*M+M*
! 400: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
! 401: WRKBL = MAX( WRKBL, 3*M+M*
! 402: $ ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
! 403: WRKBL = MAX( WRKBL, BDSPAC+3*M )
! 404: MAXWRK = WRKBL + M*N
! 405: MINWRK = 3*M + MAX( N, M*M+BDSPAC )
! 406: ELSE IF( WNTQS ) THEN
! 407: WRKBL = MAX( WRKBL, 3*M+M*
! 408: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
! 409: WRKBL = MAX( WRKBL, 3*M+M*
! 410: $ ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
! 411: MAXWRK = MAX( WRKBL, BDSPAC+3*M )
! 412: MINWRK = 3*M + MAX( N, BDSPAC )
! 413: ELSE IF( WNTQA ) THEN
! 414: WRKBL = MAX( WRKBL, 3*M+M*
! 415: $ ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
! 416: WRKBL = MAX( WRKBL, 3*M+M*
! 417: $ ILAENV( 1, 'DORMBR', 'PRT', N, N, M, -1 ) )
! 418: MAXWRK = MAX( WRKBL, BDSPAC+3*M )
! 419: MINWRK = 3*M + MAX( N, BDSPAC )
! 420: END IF
! 421: END IF
! 422: END IF
! 423: MAXWRK = MAX( MAXWRK, MINWRK )
! 424: WORK( 1 ) = MAXWRK
! 425: *
! 426: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
! 427: INFO = -12
! 428: END IF
! 429: END IF
! 430: *
! 431: IF( INFO.NE.0 ) THEN
! 432: CALL XERBLA( 'DGESDD', -INFO )
! 433: RETURN
! 434: ELSE IF( LQUERY ) THEN
! 435: RETURN
! 436: END IF
! 437: *
! 438: * Quick return if possible
! 439: *
! 440: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
! 441: RETURN
! 442: END IF
! 443: *
! 444: * Get machine constants
! 445: *
! 446: EPS = DLAMCH( 'P' )
! 447: SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
! 448: BIGNUM = ONE / SMLNUM
! 449: *
! 450: * Scale A if max element outside range [SMLNUM,BIGNUM]
! 451: *
! 452: ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
! 453: ISCL = 0
! 454: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 455: ISCL = 1
! 456: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
! 457: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 458: ISCL = 1
! 459: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
! 460: END IF
! 461: *
! 462: IF( M.GE.N ) THEN
! 463: *
! 464: * A has at least as many rows as columns. If A has sufficiently
! 465: * more rows than columns, first reduce using the QR
! 466: * decomposition (if sufficient workspace available)
! 467: *
! 468: IF( M.GE.MNTHR ) THEN
! 469: *
! 470: IF( WNTQN ) THEN
! 471: *
! 472: * Path 1 (M much larger than N, JOBZ='N')
! 473: * No singular vectors to be computed
! 474: *
! 475: ITAU = 1
! 476: NWORK = ITAU + N
! 477: *
! 478: * Compute A=Q*R
! 479: * (Workspace: need 2*N, prefer N+N*NB)
! 480: *
! 481: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 482: $ LWORK-NWORK+1, IERR )
! 483: *
! 484: * Zero out below R
! 485: *
! 486: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
! 487: IE = 1
! 488: ITAUQ = IE + N
! 489: ITAUP = ITAUQ + N
! 490: NWORK = ITAUP + N
! 491: *
! 492: * Bidiagonalize R in A
! 493: * (Workspace: need 4*N, prefer 3*N+2*N*NB)
! 494: *
! 495: CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
! 496: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 497: $ IERR )
! 498: NWORK = IE + N
! 499: *
! 500: * Perform bidiagonal SVD, computing singular values only
! 501: * (Workspace: need N+BDSPAC)
! 502: *
! 503: CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
! 504: $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
! 505: *
! 506: ELSE IF( WNTQO ) THEN
! 507: *
! 508: * Path 2 (M much larger than N, JOBZ = 'O')
! 509: * N left singular vectors to be overwritten on A and
! 510: * N right singular vectors to be computed in VT
! 511: *
! 512: IR = 1
! 513: *
! 514: * WORK(IR) is LDWRKR by N
! 515: *
! 516: IF( LWORK.GE.LDA*N+N*N+3*N+BDSPAC ) THEN
! 517: LDWRKR = LDA
! 518: ELSE
! 519: LDWRKR = ( LWORK-N*N-3*N-BDSPAC ) / N
! 520: END IF
! 521: ITAU = IR + LDWRKR*N
! 522: NWORK = ITAU + N
! 523: *
! 524: * Compute A=Q*R
! 525: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
! 526: *
! 527: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 528: $ LWORK-NWORK+1, IERR )
! 529: *
! 530: * Copy R to WORK(IR), zeroing out below it
! 531: *
! 532: CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
! 533: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
! 534: $ LDWRKR )
! 535: *
! 536: * Generate Q in A
! 537: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
! 538: *
! 539: CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
! 540: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 541: IE = ITAU
! 542: ITAUQ = IE + N
! 543: ITAUP = ITAUQ + N
! 544: NWORK = ITAUP + N
! 545: *
! 546: * Bidiagonalize R in VT, copying result to WORK(IR)
! 547: * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
! 548: *
! 549: CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
! 550: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
! 551: $ LWORK-NWORK+1, IERR )
! 552: *
! 553: * WORK(IU) is N by N
! 554: *
! 555: IU = NWORK
! 556: NWORK = IU + N*N
! 557: *
! 558: * Perform bidiagonal SVD, computing left singular vectors
! 559: * of bidiagonal matrix in WORK(IU) and computing right
! 560: * singular vectors of bidiagonal matrix in VT
! 561: * (Workspace: need N+N*N+BDSPAC)
! 562: *
! 563: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
! 564: $ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
! 565: $ INFO )
! 566: *
! 567: * Overwrite WORK(IU) by left singular vectors of R
! 568: * and VT by right singular vectors of R
! 569: * (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
! 570: *
! 571: CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
! 572: $ WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
! 573: $ LWORK-NWORK+1, IERR )
! 574: CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
! 575: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 576: $ LWORK-NWORK+1, IERR )
! 577: *
! 578: * Multiply Q in A by left singular vectors of R in
! 579: * WORK(IU), storing result in WORK(IR) and copying to A
! 580: * (Workspace: need 2*N*N, prefer N*N+M*N)
! 581: *
! 582: DO 10 I = 1, M, LDWRKR
! 583: CHUNK = MIN( M-I+1, LDWRKR )
! 584: CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
! 585: $ LDA, WORK( IU ), N, ZERO, WORK( IR ),
! 586: $ LDWRKR )
! 587: CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
! 588: $ A( I, 1 ), LDA )
! 589: 10 CONTINUE
! 590: *
! 591: ELSE IF( WNTQS ) THEN
! 592: *
! 593: * Path 3 (M much larger than N, JOBZ='S')
! 594: * N left singular vectors to be computed in U and
! 595: * N right singular vectors to be computed in VT
! 596: *
! 597: IR = 1
! 598: *
! 599: * WORK(IR) is N by N
! 600: *
! 601: LDWRKR = N
! 602: ITAU = IR + LDWRKR*N
! 603: NWORK = ITAU + N
! 604: *
! 605: * Compute A=Q*R
! 606: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
! 607: *
! 608: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 609: $ LWORK-NWORK+1, IERR )
! 610: *
! 611: * Copy R to WORK(IR), zeroing out below it
! 612: *
! 613: CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
! 614: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
! 615: $ LDWRKR )
! 616: *
! 617: * Generate Q in A
! 618: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
! 619: *
! 620: CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
! 621: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 622: IE = ITAU
! 623: ITAUQ = IE + N
! 624: ITAUP = ITAUQ + N
! 625: NWORK = ITAUP + N
! 626: *
! 627: * Bidiagonalize R in WORK(IR)
! 628: * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
! 629: *
! 630: CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
! 631: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
! 632: $ LWORK-NWORK+1, IERR )
! 633: *
! 634: * Perform bidiagonal SVD, computing left singular vectors
! 635: * of bidiagoal matrix in U and computing right singular
! 636: * vectors of bidiagonal matrix in VT
! 637: * (Workspace: need N+BDSPAC)
! 638: *
! 639: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
! 640: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
! 641: $ INFO )
! 642: *
! 643: * Overwrite U by left singular vectors of R and VT
! 644: * by right singular vectors of R
! 645: * (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
! 646: *
! 647: CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
! 648: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 649: $ LWORK-NWORK+1, IERR )
! 650: *
! 651: CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
! 652: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 653: $ LWORK-NWORK+1, IERR )
! 654: *
! 655: * Multiply Q in A by left singular vectors of R in
! 656: * WORK(IR), storing result in U
! 657: * (Workspace: need N*N)
! 658: *
! 659: CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
! 660: CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),
! 661: $ LDWRKR, ZERO, U, LDU )
! 662: *
! 663: ELSE IF( WNTQA ) THEN
! 664: *
! 665: * Path 4 (M much larger than N, JOBZ='A')
! 666: * M left singular vectors to be computed in U and
! 667: * N right singular vectors to be computed in VT
! 668: *
! 669: IU = 1
! 670: *
! 671: * WORK(IU) is N by N
! 672: *
! 673: LDWRKU = N
! 674: ITAU = IU + LDWRKU*N
! 675: NWORK = ITAU + N
! 676: *
! 677: * Compute A=Q*R, copying result to U
! 678: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
! 679: *
! 680: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 681: $ LWORK-NWORK+1, IERR )
! 682: CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
! 683: *
! 684: * Generate Q in U
! 685: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
! 686: CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
! 687: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 688: *
! 689: * Produce R in A, zeroing out other entries
! 690: *
! 691: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
! 692: IE = ITAU
! 693: ITAUQ = IE + N
! 694: ITAUP = ITAUQ + N
! 695: NWORK = ITAUP + N
! 696: *
! 697: * Bidiagonalize R in A
! 698: * (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
! 699: *
! 700: CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
! 701: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 702: $ IERR )
! 703: *
! 704: * Perform bidiagonal SVD, computing left singular vectors
! 705: * of bidiagonal matrix in WORK(IU) and computing right
! 706: * singular vectors of bidiagonal matrix in VT
! 707: * (Workspace: need N+N*N+BDSPAC)
! 708: *
! 709: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
! 710: $ VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
! 711: $ INFO )
! 712: *
! 713: * Overwrite WORK(IU) by left singular vectors of R and VT
! 714: * by right singular vectors of R
! 715: * (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
! 716: *
! 717: CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
! 718: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
! 719: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 720: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
! 721: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 722: $ LWORK-NWORK+1, IERR )
! 723: *
! 724: * Multiply Q in U by left singular vectors of R in
! 725: * WORK(IU), storing result in A
! 726: * (Workspace: need N*N)
! 727: *
! 728: CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),
! 729: $ LDWRKU, ZERO, A, LDA )
! 730: *
! 731: * Copy left singular vectors of A from A to U
! 732: *
! 733: CALL DLACPY( 'F', M, N, A, LDA, U, LDU )
! 734: *
! 735: END IF
! 736: *
! 737: ELSE
! 738: *
! 739: * M .LT. MNTHR
! 740: *
! 741: * Path 5 (M at least N, but not much larger)
! 742: * Reduce to bidiagonal form without QR decomposition
! 743: *
! 744: IE = 1
! 745: ITAUQ = IE + N
! 746: ITAUP = ITAUQ + N
! 747: NWORK = ITAUP + N
! 748: *
! 749: * Bidiagonalize A
! 750: * (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
! 751: *
! 752: CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
! 753: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 754: $ IERR )
! 755: IF( WNTQN ) THEN
! 756: *
! 757: * Perform bidiagonal SVD, only computing singular values
! 758: * (Workspace: need N+BDSPAC)
! 759: *
! 760: CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
! 761: $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
! 762: ELSE IF( WNTQO ) THEN
! 763: IU = NWORK
! 764: IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
! 765: *
! 766: * WORK( IU ) is M by N
! 767: *
! 768: LDWRKU = M
! 769: NWORK = IU + LDWRKU*N
! 770: CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),
! 771: $ LDWRKU )
! 772: ELSE
! 773: *
! 774: * WORK( IU ) is N by N
! 775: *
! 776: LDWRKU = N
! 777: NWORK = IU + LDWRKU*N
! 778: *
! 779: * WORK(IR) is LDWRKR by N
! 780: *
! 781: IR = NWORK
! 782: LDWRKR = ( LWORK-N*N-3*N ) / N
! 783: END IF
! 784: NWORK = IU + LDWRKU*N
! 785: *
! 786: * Perform bidiagonal SVD, computing left singular vectors
! 787: * of bidiagonal matrix in WORK(IU) and computing right
! 788: * singular vectors of bidiagonal matrix in VT
! 789: * (Workspace: need N+N*N+BDSPAC)
! 790: *
! 791: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
! 792: $ LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
! 793: $ IWORK, INFO )
! 794: *
! 795: * Overwrite VT by right singular vectors of A
! 796: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
! 797: *
! 798: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
! 799: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 800: $ LWORK-NWORK+1, IERR )
! 801: *
! 802: IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
! 803: *
! 804: * Overwrite WORK(IU) by left singular vectors of A
! 805: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
! 806: *
! 807: CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
! 808: $ WORK( ITAUQ ), WORK( IU ), LDWRKU,
! 809: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 810: *
! 811: * Copy left singular vectors of A from WORK(IU) to A
! 812: *
! 813: CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
! 814: ELSE
! 815: *
! 816: * Generate Q in A
! 817: * (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
! 818: *
! 819: CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
! 820: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 821: *
! 822: * Multiply Q in A by left singular vectors of
! 823: * bidiagonal matrix in WORK(IU), storing result in
! 824: * WORK(IR) and copying to A
! 825: * (Workspace: need 2*N*N, prefer N*N+M*N)
! 826: *
! 827: DO 20 I = 1, M, LDWRKR
! 828: CHUNK = MIN( M-I+1, LDWRKR )
! 829: CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
! 830: $ LDA, WORK( IU ), LDWRKU, ZERO,
! 831: $ WORK( IR ), LDWRKR )
! 832: CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
! 833: $ A( I, 1 ), LDA )
! 834: 20 CONTINUE
! 835: END IF
! 836: *
! 837: ELSE IF( WNTQS ) THEN
! 838: *
! 839: * Perform bidiagonal SVD, computing left singular vectors
! 840: * of bidiagonal matrix in U and computing right singular
! 841: * vectors of bidiagonal matrix in VT
! 842: * (Workspace: need N+BDSPAC)
! 843: *
! 844: CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU )
! 845: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
! 846: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
! 847: $ INFO )
! 848: *
! 849: * Overwrite U by left singular vectors of A and VT
! 850: * by right singular vectors of A
! 851: * (Workspace: need 3*N, prefer 2*N+N*NB)
! 852: *
! 853: CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
! 854: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 855: $ LWORK-NWORK+1, IERR )
! 856: CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
! 857: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 858: $ LWORK-NWORK+1, IERR )
! 859: ELSE IF( WNTQA ) THEN
! 860: *
! 861: * Perform bidiagonal SVD, computing left singular vectors
! 862: * of bidiagonal matrix in U and computing right singular
! 863: * vectors of bidiagonal matrix in VT
! 864: * (Workspace: need N+BDSPAC)
! 865: *
! 866: CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU )
! 867: CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
! 868: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
! 869: $ INFO )
! 870: *
! 871: * Set the right corner of U to identity matrix
! 872: *
! 873: IF( M.GT.N ) THEN
! 874: CALL DLASET( 'F', M-N, M-N, ZERO, ONE, U( N+1, N+1 ),
! 875: $ LDU )
! 876: END IF
! 877: *
! 878: * Overwrite U by left singular vectors of A and VT
! 879: * by right singular vectors of A
! 880: * (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB)
! 881: *
! 882: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
! 883: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 884: $ LWORK-NWORK+1, IERR )
! 885: CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
! 886: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 887: $ LWORK-NWORK+1, IERR )
! 888: END IF
! 889: *
! 890: END IF
! 891: *
! 892: ELSE
! 893: *
! 894: * A has more columns than rows. If A has sufficiently more
! 895: * columns than rows, first reduce using the LQ decomposition (if
! 896: * sufficient workspace available)
! 897: *
! 898: IF( N.GE.MNTHR ) THEN
! 899: *
! 900: IF( WNTQN ) THEN
! 901: *
! 902: * Path 1t (N much larger than M, JOBZ='N')
! 903: * No singular vectors to be computed
! 904: *
! 905: ITAU = 1
! 906: NWORK = ITAU + M
! 907: *
! 908: * Compute A=L*Q
! 909: * (Workspace: need 2*M, prefer M+M*NB)
! 910: *
! 911: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 912: $ LWORK-NWORK+1, IERR )
! 913: *
! 914: * Zero out above L
! 915: *
! 916: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
! 917: IE = 1
! 918: ITAUQ = IE + M
! 919: ITAUP = ITAUQ + M
! 920: NWORK = ITAUP + M
! 921: *
! 922: * Bidiagonalize L in A
! 923: * (Workspace: need 4*M, prefer 3*M+2*M*NB)
! 924: *
! 925: CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
! 926: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 927: $ IERR )
! 928: NWORK = IE + M
! 929: *
! 930: * Perform bidiagonal SVD, computing singular values only
! 931: * (Workspace: need M+BDSPAC)
! 932: *
! 933: CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
! 934: $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
! 935: *
! 936: ELSE IF( WNTQO ) THEN
! 937: *
! 938: * Path 2t (N much larger than M, JOBZ='O')
! 939: * M right singular vectors to be overwritten on A and
! 940: * M left singular vectors to be computed in U
! 941: *
! 942: IVT = 1
! 943: *
! 944: * IVT is M by M
! 945: *
! 946: IL = IVT + M*M
! 947: IF( LWORK.GE.M*N+M*M+3*M+BDSPAC ) THEN
! 948: *
! 949: * WORK(IL) is M by N
! 950: *
! 951: LDWRKL = M
! 952: CHUNK = N
! 953: ELSE
! 954: LDWRKL = M
! 955: CHUNK = ( LWORK-M*M ) / M
! 956: END IF
! 957: ITAU = IL + LDWRKL*M
! 958: NWORK = ITAU + M
! 959: *
! 960: * Compute A=L*Q
! 961: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
! 962: *
! 963: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 964: $ LWORK-NWORK+1, IERR )
! 965: *
! 966: * Copy L to WORK(IL), zeroing about above it
! 967: *
! 968: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
! 969: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
! 970: $ WORK( IL+LDWRKL ), LDWRKL )
! 971: *
! 972: * Generate Q in A
! 973: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
! 974: *
! 975: CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
! 976: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 977: IE = ITAU
! 978: ITAUQ = IE + M
! 979: ITAUP = ITAUQ + M
! 980: NWORK = ITAUP + M
! 981: *
! 982: * Bidiagonalize L in WORK(IL)
! 983: * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
! 984: *
! 985: CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
! 986: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
! 987: $ LWORK-NWORK+1, IERR )
! 988: *
! 989: * Perform bidiagonal SVD, computing left singular vectors
! 990: * of bidiagonal matrix in U, and computing right singular
! 991: * vectors of bidiagonal matrix in WORK(IVT)
! 992: * (Workspace: need M+M*M+BDSPAC)
! 993: *
! 994: CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
! 995: $ WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
! 996: $ IWORK, INFO )
! 997: *
! 998: * Overwrite U by left singular vectors of L and WORK(IVT)
! 999: * by right singular vectors of L
! 1000: * (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
! 1001: *
! 1002: CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
! 1003: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1004: $ LWORK-NWORK+1, IERR )
! 1005: CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
! 1006: $ WORK( ITAUP ), WORK( IVT ), M,
! 1007: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1008: *
! 1009: * Multiply right singular vectors of L in WORK(IVT) by Q
! 1010: * in A, storing result in WORK(IL) and copying to A
! 1011: * (Workspace: need 2*M*M, prefer M*M+M*N)
! 1012: *
! 1013: DO 30 I = 1, N, CHUNK
! 1014: BLK = MIN( N-I+1, CHUNK )
! 1015: CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,
! 1016: $ A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
! 1017: CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
! 1018: $ A( 1, I ), LDA )
! 1019: 30 CONTINUE
! 1020: *
! 1021: ELSE IF( WNTQS ) THEN
! 1022: *
! 1023: * Path 3t (N much larger than M, JOBZ='S')
! 1024: * M right singular vectors to be computed in VT and
! 1025: * M left singular vectors to be computed in U
! 1026: *
! 1027: IL = 1
! 1028: *
! 1029: * WORK(IL) is M by M
! 1030: *
! 1031: LDWRKL = M
! 1032: ITAU = IL + LDWRKL*M
! 1033: NWORK = ITAU + M
! 1034: *
! 1035: * Compute A=L*Q
! 1036: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
! 1037: *
! 1038: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 1039: $ LWORK-NWORK+1, IERR )
! 1040: *
! 1041: * Copy L to WORK(IL), zeroing out above it
! 1042: *
! 1043: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
! 1044: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
! 1045: $ WORK( IL+LDWRKL ), LDWRKL )
! 1046: *
! 1047: * Generate Q in A
! 1048: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
! 1049: *
! 1050: CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
! 1051: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1052: IE = ITAU
! 1053: ITAUQ = IE + M
! 1054: ITAUP = ITAUQ + M
! 1055: NWORK = ITAUP + M
! 1056: *
! 1057: * Bidiagonalize L in WORK(IU), copying result to U
! 1058: * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
! 1059: *
! 1060: CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
! 1061: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
! 1062: $ LWORK-NWORK+1, IERR )
! 1063: *
! 1064: * Perform bidiagonal SVD, computing left singular vectors
! 1065: * of bidiagonal matrix in U and computing right singular
! 1066: * vectors of bidiagonal matrix in VT
! 1067: * (Workspace: need M+BDSPAC)
! 1068: *
! 1069: CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
! 1070: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
! 1071: $ INFO )
! 1072: *
! 1073: * Overwrite U by left singular vectors of L and VT
! 1074: * by right singular vectors of L
! 1075: * (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
! 1076: *
! 1077: CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
! 1078: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1079: $ LWORK-NWORK+1, IERR )
! 1080: CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
! 1081: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 1082: $ LWORK-NWORK+1, IERR )
! 1083: *
! 1084: * Multiply right singular vectors of L in WORK(IL) by
! 1085: * Q in A, storing result in VT
! 1086: * (Workspace: need M*M)
! 1087: *
! 1088: CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
! 1089: CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,
! 1090: $ A, LDA, ZERO, VT, LDVT )
! 1091: *
! 1092: ELSE IF( WNTQA ) THEN
! 1093: *
! 1094: * Path 4t (N much larger than M, JOBZ='A')
! 1095: * N right singular vectors to be computed in VT and
! 1096: * M left singular vectors to be computed in U
! 1097: *
! 1098: IVT = 1
! 1099: *
! 1100: * WORK(IVT) is M by M
! 1101: *
! 1102: LDWKVT = M
! 1103: ITAU = IVT + LDWKVT*M
! 1104: NWORK = ITAU + M
! 1105: *
! 1106: * Compute A=L*Q, copying result to VT
! 1107: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
! 1108: *
! 1109: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 1110: $ LWORK-NWORK+1, IERR )
! 1111: CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
! 1112: *
! 1113: * Generate Q in VT
! 1114: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
! 1115: *
! 1116: CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
! 1117: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1118: *
! 1119: * Produce L in A, zeroing out other entries
! 1120: *
! 1121: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
! 1122: IE = ITAU
! 1123: ITAUQ = IE + M
! 1124: ITAUP = ITAUQ + M
! 1125: NWORK = ITAUP + M
! 1126: *
! 1127: * Bidiagonalize L in A
! 1128: * (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
! 1129: *
! 1130: CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
! 1131: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 1132: $ IERR )
! 1133: *
! 1134: * Perform bidiagonal SVD, computing left singular vectors
! 1135: * of bidiagonal matrix in U and computing right singular
! 1136: * vectors of bidiagonal matrix in WORK(IVT)
! 1137: * (Workspace: need M+M*M+BDSPAC)
! 1138: *
! 1139: CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
! 1140: $ WORK( IVT ), LDWKVT, DUM, IDUM,
! 1141: $ WORK( NWORK ), IWORK, INFO )
! 1142: *
! 1143: * Overwrite U by left singular vectors of L and WORK(IVT)
! 1144: * by right singular vectors of L
! 1145: * (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
! 1146: *
! 1147: CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
! 1148: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1149: $ LWORK-NWORK+1, IERR )
! 1150: CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
! 1151: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
! 1152: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1153: *
! 1154: * Multiply right singular vectors of L in WORK(IVT) by
! 1155: * Q in VT, storing result in A
! 1156: * (Workspace: need M*M)
! 1157: *
! 1158: CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,
! 1159: $ VT, LDVT, ZERO, A, LDA )
! 1160: *
! 1161: * Copy right singular vectors of A from A to VT
! 1162: *
! 1163: CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT )
! 1164: *
! 1165: END IF
! 1166: *
! 1167: ELSE
! 1168: *
! 1169: * N .LT. MNTHR
! 1170: *
! 1171: * Path 5t (N greater than M, but not much larger)
! 1172: * Reduce to bidiagonal form without LQ decomposition
! 1173: *
! 1174: IE = 1
! 1175: ITAUQ = IE + M
! 1176: ITAUP = ITAUQ + M
! 1177: NWORK = ITAUP + M
! 1178: *
! 1179: * Bidiagonalize A
! 1180: * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
! 1181: *
! 1182: CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
! 1183: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 1184: $ IERR )
! 1185: IF( WNTQN ) THEN
! 1186: *
! 1187: * Perform bidiagonal SVD, only computing singular values
! 1188: * (Workspace: need M+BDSPAC)
! 1189: *
! 1190: CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
! 1191: $ DUM, IDUM, WORK( NWORK ), IWORK, INFO )
! 1192: ELSE IF( WNTQO ) THEN
! 1193: LDWKVT = M
! 1194: IVT = NWORK
! 1195: IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
! 1196: *
! 1197: * WORK( IVT ) is M by N
! 1198: *
! 1199: CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),
! 1200: $ LDWKVT )
! 1201: NWORK = IVT + LDWKVT*N
! 1202: ELSE
! 1203: *
! 1204: * WORK( IVT ) is M by M
! 1205: *
! 1206: NWORK = IVT + LDWKVT*M
! 1207: IL = NWORK
! 1208: *
! 1209: * WORK(IL) is M by CHUNK
! 1210: *
! 1211: CHUNK = ( LWORK-M*M-3*M ) / M
! 1212: END IF
! 1213: *
! 1214: * Perform bidiagonal SVD, computing left singular vectors
! 1215: * of bidiagonal matrix in U and computing right singular
! 1216: * vectors of bidiagonal matrix in WORK(IVT)
! 1217: * (Workspace: need M*M+BDSPAC)
! 1218: *
! 1219: CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,
! 1220: $ WORK( IVT ), LDWKVT, DUM, IDUM,
! 1221: $ WORK( NWORK ), IWORK, INFO )
! 1222: *
! 1223: * Overwrite U by left singular vectors of A
! 1224: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
! 1225: *
! 1226: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
! 1227: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1228: $ LWORK-NWORK+1, IERR )
! 1229: *
! 1230: IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
! 1231: *
! 1232: * Overwrite WORK(IVT) by left singular vectors of A
! 1233: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
! 1234: *
! 1235: CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
! 1236: $ WORK( ITAUP ), WORK( IVT ), LDWKVT,
! 1237: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1238: *
! 1239: * Copy right singular vectors of A from WORK(IVT) to A
! 1240: *
! 1241: CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
! 1242: ELSE
! 1243: *
! 1244: * Generate P**T in A
! 1245: * (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
! 1246: *
! 1247: CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
! 1248: $ WORK( NWORK ), LWORK-NWORK+1, IERR )
! 1249: *
! 1250: * Multiply Q in A by right singular vectors of
! 1251: * bidiagonal matrix in WORK(IVT), storing result in
! 1252: * WORK(IL) and copying to A
! 1253: * (Workspace: need 2*M*M, prefer M*M+M*N)
! 1254: *
! 1255: DO 40 I = 1, N, CHUNK
! 1256: BLK = MIN( N-I+1, CHUNK )
! 1257: CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),
! 1258: $ LDWKVT, A( 1, I ), LDA, ZERO,
! 1259: $ WORK( IL ), M )
! 1260: CALL DLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ),
! 1261: $ LDA )
! 1262: 40 CONTINUE
! 1263: END IF
! 1264: ELSE IF( WNTQS ) THEN
! 1265: *
! 1266: * Perform bidiagonal SVD, computing left singular vectors
! 1267: * of bidiagonal matrix in U and computing right singular
! 1268: * vectors of bidiagonal matrix in VT
! 1269: * (Workspace: need M+BDSPAC)
! 1270: *
! 1271: CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
! 1272: CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
! 1273: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
! 1274: $ INFO )
! 1275: *
! 1276: * Overwrite U by left singular vectors of A and VT
! 1277: * by right singular vectors of A
! 1278: * (Workspace: need 3*M, prefer 2*M+M*NB)
! 1279: *
! 1280: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
! 1281: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1282: $ LWORK-NWORK+1, IERR )
! 1283: CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
! 1284: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 1285: $ LWORK-NWORK+1, IERR )
! 1286: ELSE IF( WNTQA ) THEN
! 1287: *
! 1288: * Perform bidiagonal SVD, computing left singular vectors
! 1289: * of bidiagonal matrix in U and computing right singular
! 1290: * vectors of bidiagonal matrix in VT
! 1291: * (Workspace: need M+BDSPAC)
! 1292: *
! 1293: CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
! 1294: CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
! 1295: $ LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
! 1296: $ INFO )
! 1297: *
! 1298: * Set the right corner of VT to identity matrix
! 1299: *
! 1300: IF( N.GT.M ) THEN
! 1301: CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT( M+1, M+1 ),
! 1302: $ LDVT )
! 1303: END IF
! 1304: *
! 1305: * Overwrite U by left singular vectors of A and VT
! 1306: * by right singular vectors of A
! 1307: * (Workspace: need 2*M+N, prefer 2*M+N*NB)
! 1308: *
! 1309: CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
! 1310: $ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
! 1311: $ LWORK-NWORK+1, IERR )
! 1312: CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
! 1313: $ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
! 1314: $ LWORK-NWORK+1, IERR )
! 1315: END IF
! 1316: *
! 1317: END IF
! 1318: *
! 1319: END IF
! 1320: *
! 1321: * Undo scaling if necessary
! 1322: *
! 1323: IF( ISCL.EQ.1 ) THEN
! 1324: IF( ANRM.GT.BIGNUM )
! 1325: $ CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
! 1326: $ IERR )
! 1327: IF( ANRM.LT.SMLNUM )
! 1328: $ CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
! 1329: $ IERR )
! 1330: END IF
! 1331: *
! 1332: * Return optimal workspace in WORK(1)
! 1333: *
! 1334: WORK( 1 ) = MAXWRK
! 1335: *
! 1336: RETURN
! 1337: *
! 1338: * End of DGESDD
! 1339: *
! 1340: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgesdd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack