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