Annotation of rpl/lapack/lapack/dgelsd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
! 2: $ WORK, LWORK, IWORK, INFO )
! 3: *
! 4: * -- LAPACK driver routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
! 11: DOUBLE PRECISION RCOND
! 12: * ..
! 13: * .. Array Arguments ..
! 14: INTEGER IWORK( * )
! 15: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * DGELSD computes the minimum-norm solution to a real linear least
! 22: * squares problem:
! 23: * minimize 2-norm(| b - A*x |)
! 24: * using the singular value decomposition (SVD) of A. A is an M-by-N
! 25: * matrix which may be rank-deficient.
! 26: *
! 27: * Several right hand side vectors b and solution vectors x can be
! 28: * handled in a single call; they are stored as the columns of the
! 29: * M-by-NRHS right hand side matrix B and the N-by-NRHS solution
! 30: * matrix X.
! 31: *
! 32: * The problem is solved in three steps:
! 33: * (1) Reduce the coefficient matrix A to bidiagonal form with
! 34: * Householder transformations, reducing the original problem
! 35: * into a "bidiagonal least squares problem" (BLS)
! 36: * (2) Solve the BLS using a divide and conquer approach.
! 37: * (3) Apply back all the Householder tranformations to solve
! 38: * the original least squares problem.
! 39: *
! 40: * The effective rank of A is determined by treating as zero those
! 41: * singular values which are less than RCOND times the largest singular
! 42: * value.
! 43: *
! 44: * The divide and conquer algorithm makes very mild assumptions about
! 45: * floating point arithmetic. It will work on machines with a guard
! 46: * digit in add/subtract, or on those binary machines without guard
! 47: * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
! 48: * Cray-2. It could conceivably fail on hexadecimal or decimal machines
! 49: * without guard digits, but we know of none.
! 50: *
! 51: * Arguments
! 52: * =========
! 53: *
! 54: * M (input) INTEGER
! 55: * The number of rows of A. M >= 0.
! 56: *
! 57: * N (input) INTEGER
! 58: * The number of columns of A. N >= 0.
! 59: *
! 60: * NRHS (input) INTEGER
! 61: * The number of right hand sides, i.e., the number of columns
! 62: * of the matrices B and X. NRHS >= 0.
! 63: *
! 64: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
! 65: * On entry, the M-by-N matrix A.
! 66: * On exit, A has been destroyed.
! 67: *
! 68: * LDA (input) INTEGER
! 69: * The leading dimension of the array A. LDA >= max(1,M).
! 70: *
! 71: * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
! 72: * On entry, the M-by-NRHS right hand side matrix B.
! 73: * On exit, B is overwritten by the N-by-NRHS solution
! 74: * matrix X. If m >= n and RANK = n, the residual
! 75: * sum-of-squares for the solution in the i-th column is given
! 76: * by the sum of squares of elements n+1:m in that column.
! 77: *
! 78: * LDB (input) INTEGER
! 79: * The leading dimension of the array B. LDB >= max(1,max(M,N)).
! 80: *
! 81: * S (output) DOUBLE PRECISION array, dimension (min(M,N))
! 82: * The singular values of A in decreasing order.
! 83: * The condition number of A in the 2-norm = S(1)/S(min(m,n)).
! 84: *
! 85: * RCOND (input) DOUBLE PRECISION
! 86: * RCOND is used to determine the effective rank of A.
! 87: * Singular values S(i) <= RCOND*S(1) are treated as zero.
! 88: * If RCOND < 0, machine precision is used instead.
! 89: *
! 90: * RANK (output) INTEGER
! 91: * The effective rank of A, i.e., the number of singular values
! 92: * which are greater than RCOND*S(1).
! 93: *
! 94: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 95: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 96: *
! 97: * LWORK (input) INTEGER
! 98: * The dimension of the array WORK. LWORK must be at least 1.
! 99: * The exact minimum amount of workspace needed depends on M,
! 100: * N and NRHS. As long as LWORK is at least
! 101: * 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2,
! 102: * if M is greater than or equal to N or
! 103: * 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2,
! 104: * if M is less than N, the code will execute correctly.
! 105: * SMLSIZ is returned by ILAENV and is equal to the maximum
! 106: * size of the subproblems at the bottom of the computation
! 107: * tree (usually about 25), and
! 108: * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
! 109: * For good performance, LWORK should generally be larger.
! 110: *
! 111: * If LWORK = -1, then a workspace query is assumed; the routine
! 112: * only calculates the optimal size of the WORK array, returns
! 113: * this value as the first entry of the WORK array, and no error
! 114: * message related to LWORK is issued by XERBLA.
! 115: *
! 116: * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK))
! 117: * LIWORK >= 3 * MINMN * NLVL + 11 * MINMN,
! 118: * where MINMN = MIN( M,N ).
! 119: *
! 120: * INFO (output) INTEGER
! 121: * = 0: successful exit
! 122: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 123: * > 0: the algorithm for computing the SVD failed to converge;
! 124: * if INFO = i, i off-diagonal elements of an intermediate
! 125: * bidiagonal form did not converge to zero.
! 126: *
! 127: * Further Details
! 128: * ===============
! 129: *
! 130: * Based on contributions by
! 131: * Ming Gu and Ren-Cang Li, Computer Science Division, University of
! 132: * California at Berkeley, USA
! 133: * Osni Marques, LBNL/NERSC, USA
! 134: *
! 135: * =====================================================================
! 136: *
! 137: * .. Parameters ..
! 138: DOUBLE PRECISION ZERO, ONE, TWO
! 139: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
! 140: * ..
! 141: * .. Local Scalars ..
! 142: LOGICAL LQUERY
! 143: INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
! 144: $ LDWORK, MAXMN, MAXWRK, MINMN, MINWRK, MM,
! 145: $ MNTHR, NLVL, NWORK, SMLSIZ, WLALSD
! 146: DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
! 147: * ..
! 148: * .. External Subroutines ..
! 149: EXTERNAL DGEBRD, DGELQF, DGEQRF, DLABAD, DLACPY, DLALSD,
! 150: $ DLASCL, DLASET, DORMBR, DORMLQ, DORMQR, XERBLA
! 151: * ..
! 152: * .. External Functions ..
! 153: INTEGER ILAENV
! 154: DOUBLE PRECISION DLAMCH, DLANGE
! 155: EXTERNAL ILAENV, DLAMCH, DLANGE
! 156: * ..
! 157: * .. Intrinsic Functions ..
! 158: INTRINSIC DBLE, INT, LOG, MAX, MIN
! 159: * ..
! 160: * .. Executable Statements ..
! 161: *
! 162: * Test the input arguments.
! 163: *
! 164: INFO = 0
! 165: MINMN = MIN( M, N )
! 166: MAXMN = MAX( M, N )
! 167: MNTHR = ILAENV( 6, 'DGELSD', ' ', M, N, NRHS, -1 )
! 168: LQUERY = ( LWORK.EQ.-1 )
! 169: IF( M.LT.0 ) THEN
! 170: INFO = -1
! 171: ELSE IF( N.LT.0 ) THEN
! 172: INFO = -2
! 173: ELSE IF( NRHS.LT.0 ) THEN
! 174: INFO = -3
! 175: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 176: INFO = -5
! 177: ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
! 178: INFO = -7
! 179: END IF
! 180: *
! 181: SMLSIZ = ILAENV( 9, 'DGELSD', ' ', 0, 0, 0, 0 )
! 182: *
! 183: * Compute workspace.
! 184: * (Note: Comments in the code beginning "Workspace:" describe the
! 185: * minimal amount of workspace needed at that point in the code,
! 186: * as well as the preferred amount for good performance.
! 187: * NB refers to the optimal block size for the immediately
! 188: * following subroutine, as returned by ILAENV.)
! 189: *
! 190: MINWRK = 1
! 191: MINMN = MAX( 1, MINMN )
! 192: NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ+1 ) ) /
! 193: $ LOG( TWO ) ) + 1, 0 )
! 194: *
! 195: IF( INFO.EQ.0 ) THEN
! 196: MAXWRK = 0
! 197: MM = M
! 198: IF( M.GE.N .AND. M.GE.MNTHR ) THEN
! 199: *
! 200: * Path 1a - overdetermined, with many more rows than columns.
! 201: *
! 202: MM = N
! 203: MAXWRK = MAX( MAXWRK, N+N*ILAENV( 1, 'DGEQRF', ' ', M, N,
! 204: $ -1, -1 ) )
! 205: MAXWRK = MAX( MAXWRK, N+NRHS*
! 206: $ ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N, -1 ) )
! 207: END IF
! 208: IF( M.GE.N ) THEN
! 209: *
! 210: * Path 1 - overdetermined or exactly determined.
! 211: *
! 212: MAXWRK = MAX( MAXWRK, 3*N+( MM+N )*
! 213: $ ILAENV( 1, 'DGEBRD', ' ', MM, N, -1, -1 ) )
! 214: MAXWRK = MAX( MAXWRK, 3*N+NRHS*
! 215: $ ILAENV( 1, 'DORMBR', 'QLT', MM, NRHS, N, -1 ) )
! 216: MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
! 217: $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, N, -1 ) )
! 218: WLALSD = 9*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS+(SMLSIZ+1)**2
! 219: MAXWRK = MAX( MAXWRK, 3*N+WLALSD )
! 220: MINWRK = MAX( 3*N+MM, 3*N+NRHS, 3*N+WLALSD )
! 221: END IF
! 222: IF( N.GT.M ) THEN
! 223: WLALSD = 9*M+2*M*SMLSIZ+8*M*NLVL+M*NRHS+(SMLSIZ+1)**2
! 224: IF( N.GE.MNTHR ) THEN
! 225: *
! 226: * Path 2a - underdetermined, with many more columns
! 227: * than rows.
! 228: *
! 229: MAXWRK = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
! 230: MAXWRK = MAX( MAXWRK, M*M+4*M+2*M*
! 231: $ ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
! 232: MAXWRK = MAX( MAXWRK, M*M+4*M+NRHS*
! 233: $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, M, -1 ) )
! 234: MAXWRK = MAX( MAXWRK, M*M+4*M+( M-1 )*
! 235: $ ILAENV( 1, 'DORMBR', 'PLN', M, NRHS, M, -1 ) )
! 236: IF( NRHS.GT.1 ) THEN
! 237: MAXWRK = MAX( MAXWRK, M*M+M+M*NRHS )
! 238: ELSE
! 239: MAXWRK = MAX( MAXWRK, M*M+2*M )
! 240: END IF
! 241: MAXWRK = MAX( MAXWRK, M+NRHS*
! 242: $ ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M, -1 ) )
! 243: MAXWRK = MAX( MAXWRK, M*M+4*M+WLALSD )
! 244: ! XXX: Ensure the Path 2a case below is triggered. The workspace
! 245: ! calculation should use queries for all routines eventually.
! 246: MAXWRK = MAX( MAXWRK,
! 247: $ 4*M+M*M+MAX( M, 2*M-4, NRHS, N-3*M ) )
! 248: ELSE
! 249: *
! 250: * Path 2 - remaining underdetermined cases.
! 251: *
! 252: MAXWRK = 3*M + ( N+M )*ILAENV( 1, 'DGEBRD', ' ', M, N,
! 253: $ -1, -1 )
! 254: MAXWRK = MAX( MAXWRK, 3*M+NRHS*
! 255: $ ILAENV( 1, 'DORMBR', 'QLT', M, NRHS, N, -1 ) )
! 256: MAXWRK = MAX( MAXWRK, 3*M+M*
! 257: $ ILAENV( 1, 'DORMBR', 'PLN', N, NRHS, M, -1 ) )
! 258: MAXWRK = MAX( MAXWRK, 3*M+WLALSD )
! 259: END IF
! 260: MINWRK = MAX( 3*M+NRHS, 3*M+M, 3*M+WLALSD )
! 261: END IF
! 262: MINWRK = MIN( MINWRK, MAXWRK )
! 263: WORK( 1 ) = MAXWRK
! 264: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
! 265: INFO = -12
! 266: END IF
! 267: END IF
! 268: *
! 269: IF( INFO.NE.0 ) THEN
! 270: CALL XERBLA( 'DGELSD', -INFO )
! 271: RETURN
! 272: ELSE IF( LQUERY ) THEN
! 273: GO TO 10
! 274: END IF
! 275: *
! 276: * Quick return if possible.
! 277: *
! 278: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
! 279: RANK = 0
! 280: RETURN
! 281: END IF
! 282: *
! 283: * Get machine parameters.
! 284: *
! 285: EPS = DLAMCH( 'P' )
! 286: SFMIN = DLAMCH( 'S' )
! 287: SMLNUM = SFMIN / EPS
! 288: BIGNUM = ONE / SMLNUM
! 289: CALL DLABAD( SMLNUM, BIGNUM )
! 290: *
! 291: * Scale A if max entry outside range [SMLNUM,BIGNUM].
! 292: *
! 293: ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
! 294: IASCL = 0
! 295: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 296: *
! 297: * Scale matrix norm up to SMLNUM.
! 298: *
! 299: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
! 300: IASCL = 1
! 301: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 302: *
! 303: * Scale matrix norm down to BIGNUM.
! 304: *
! 305: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
! 306: IASCL = 2
! 307: ELSE IF( ANRM.EQ.ZERO ) THEN
! 308: *
! 309: * Matrix all zero. Return zero solution.
! 310: *
! 311: CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
! 312: CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, 1 )
! 313: RANK = 0
! 314: GO TO 10
! 315: END IF
! 316: *
! 317: * Scale B if max entry outside range [SMLNUM,BIGNUM].
! 318: *
! 319: BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
! 320: IBSCL = 0
! 321: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
! 322: *
! 323: * Scale matrix norm up to SMLNUM.
! 324: *
! 325: CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
! 326: IBSCL = 1
! 327: ELSE IF( BNRM.GT.BIGNUM ) THEN
! 328: *
! 329: * Scale matrix norm down to BIGNUM.
! 330: *
! 331: CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
! 332: IBSCL = 2
! 333: END IF
! 334: *
! 335: * If M < N make sure certain entries of B are zero.
! 336: *
! 337: IF( M.LT.N )
! 338: $ CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
! 339: *
! 340: * Overdetermined case.
! 341: *
! 342: IF( M.GE.N ) THEN
! 343: *
! 344: * Path 1 - overdetermined or exactly determined.
! 345: *
! 346: MM = M
! 347: IF( M.GE.MNTHR ) THEN
! 348: *
! 349: * Path 1a - overdetermined, with many more rows than columns.
! 350: *
! 351: MM = N
! 352: ITAU = 1
! 353: NWORK = ITAU + N
! 354: *
! 355: * Compute A=Q*R.
! 356: * (Workspace: need 2*N, prefer N+N*NB)
! 357: *
! 358: CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 359: $ LWORK-NWORK+1, INFO )
! 360: *
! 361: * Multiply B by transpose(Q).
! 362: * (Workspace: need N+NRHS, prefer N+NRHS*NB)
! 363: *
! 364: CALL DORMQR( 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAU ), B,
! 365: $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
! 366: *
! 367: * Zero out below R.
! 368: *
! 369: IF( N.GT.1 ) THEN
! 370: CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
! 371: END IF
! 372: END IF
! 373: *
! 374: IE = 1
! 375: ITAUQ = IE + N
! 376: ITAUP = ITAUQ + N
! 377: NWORK = ITAUP + N
! 378: *
! 379: * Bidiagonalize R in A.
! 380: * (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
! 381: *
! 382: CALL DGEBRD( MM, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
! 383: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 384: $ INFO )
! 385: *
! 386: * Multiply B by transpose of left bidiagonalizing vectors of R.
! 387: * (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
! 388: *
! 389: CALL DORMBR( 'Q', 'L', 'T', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
! 390: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
! 391: *
! 392: * Solve the bidiagonal least squares problem.
! 393: *
! 394: CALL DLALSD( 'U', SMLSIZ, N, NRHS, S, WORK( IE ), B, LDB,
! 395: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
! 396: IF( INFO.NE.0 ) THEN
! 397: GO TO 10
! 398: END IF
! 399: *
! 400: * Multiply B by right bidiagonalizing vectors of R.
! 401: *
! 402: CALL DORMBR( 'P', 'L', 'N', N, NRHS, N, A, LDA, WORK( ITAUP ),
! 403: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
! 404: *
! 405: ELSE IF( N.GE.MNTHR .AND. LWORK.GE.4*M+M*M+
! 406: $ MAX( M, 2*M-4, NRHS, N-3*M, WLALSD ) ) THEN
! 407: *
! 408: * Path 2a - underdetermined, with many more columns than rows
! 409: * and sufficient workspace for an efficient algorithm.
! 410: *
! 411: LDWORK = M
! 412: IF( LWORK.GE.MAX( 4*M+M*LDA+MAX( M, 2*M-4, NRHS, N-3*M ),
! 413: $ M*LDA+M+M*NRHS, 4*M+M*LDA+WLALSD ) )LDWORK = LDA
! 414: ITAU = 1
! 415: NWORK = M + 1
! 416: *
! 417: * Compute A=L*Q.
! 418: * (Workspace: need 2*M, prefer M+M*NB)
! 419: *
! 420: CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
! 421: $ LWORK-NWORK+1, INFO )
! 422: IL = NWORK
! 423: *
! 424: * Copy L to WORK(IL), zeroing out above its diagonal.
! 425: *
! 426: CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
! 427: CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, WORK( IL+LDWORK ),
! 428: $ LDWORK )
! 429: IE = IL + LDWORK*M
! 430: ITAUQ = IE + M
! 431: ITAUP = ITAUQ + M
! 432: NWORK = ITAUP + M
! 433: *
! 434: * Bidiagonalize L in WORK(IL).
! 435: * (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB)
! 436: *
! 437: CALL DGEBRD( M, M, WORK( IL ), LDWORK, S, WORK( IE ),
! 438: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
! 439: $ LWORK-NWORK+1, INFO )
! 440: *
! 441: * Multiply B by transpose of left bidiagonalizing vectors of L.
! 442: * (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB)
! 443: *
! 444: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, M, WORK( IL ), LDWORK,
! 445: $ WORK( ITAUQ ), B, LDB, WORK( NWORK ),
! 446: $ LWORK-NWORK+1, INFO )
! 447: *
! 448: * Solve the bidiagonal least squares problem.
! 449: *
! 450: CALL DLALSD( 'U', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB,
! 451: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
! 452: IF( INFO.NE.0 ) THEN
! 453: GO TO 10
! 454: END IF
! 455: *
! 456: * Multiply B by right bidiagonalizing vectors of L.
! 457: *
! 458: CALL DORMBR( 'P', 'L', 'N', M, NRHS, M, WORK( IL ), LDWORK,
! 459: $ WORK( ITAUP ), B, LDB, WORK( NWORK ),
! 460: $ LWORK-NWORK+1, INFO )
! 461: *
! 462: * Zero out below first M rows of B.
! 463: *
! 464: CALL DLASET( 'F', N-M, NRHS, ZERO, ZERO, B( M+1, 1 ), LDB )
! 465: NWORK = ITAU + M
! 466: *
! 467: * Multiply transpose(Q) by B.
! 468: * (Workspace: need M+NRHS, prefer M+NRHS*NB)
! 469: *
! 470: CALL DORMLQ( 'L', 'T', N, NRHS, M, A, LDA, WORK( ITAU ), B,
! 471: $ LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
! 472: *
! 473: ELSE
! 474: *
! 475: * Path 2 - remaining underdetermined cases.
! 476: *
! 477: IE = 1
! 478: ITAUQ = IE + M
! 479: ITAUP = ITAUQ + M
! 480: NWORK = ITAUP + M
! 481: *
! 482: * Bidiagonalize A.
! 483: * (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
! 484: *
! 485: CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
! 486: $ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
! 487: $ INFO )
! 488: *
! 489: * Multiply B by transpose of left bidiagonalizing vectors.
! 490: * (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB)
! 491: *
! 492: CALL DORMBR( 'Q', 'L', 'T', M, NRHS, N, A, LDA, WORK( ITAUQ ),
! 493: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
! 494: *
! 495: * Solve the bidiagonal least squares problem.
! 496: *
! 497: CALL DLALSD( 'L', SMLSIZ, M, NRHS, S, WORK( IE ), B, LDB,
! 498: $ RCOND, RANK, WORK( NWORK ), IWORK, INFO )
! 499: IF( INFO.NE.0 ) THEN
! 500: GO TO 10
! 501: END IF
! 502: *
! 503: * Multiply B by right bidiagonalizing vectors of A.
! 504: *
! 505: CALL DORMBR( 'P', 'L', 'N', N, NRHS, M, A, LDA, WORK( ITAUP ),
! 506: $ B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
! 507: *
! 508: END IF
! 509: *
! 510: * Undo scaling.
! 511: *
! 512: IF( IASCL.EQ.1 ) THEN
! 513: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
! 514: CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
! 515: $ INFO )
! 516: ELSE IF( IASCL.EQ.2 ) THEN
! 517: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
! 518: CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
! 519: $ INFO )
! 520: END IF
! 521: IF( IBSCL.EQ.1 ) THEN
! 522: CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
! 523: ELSE IF( IBSCL.EQ.2 ) THEN
! 524: CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
! 525: END IF
! 526: *
! 527: 10 CONTINUE
! 528: WORK( 1 ) = MAXWRK
! 529: RETURN
! 530: *
! 531: * End of DGELSD
! 532: *
! 533: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgelsd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack