Annotation of rpl/lapack/lapack/zgelss.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
! 2: $ WORK, LWORK, RWORK, 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: DOUBLE PRECISION RWORK( * ), S( * )
! 15: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * ZGELSS computes the minimum norm solution to a complex linear
! 22: * least squares problem:
! 23: *
! 24: * Minimize 2-norm(| b - A*x |).
! 25: *
! 26: * using the singular value decomposition (SVD) of A. A is an M-by-N
! 27: * matrix which may be rank-deficient.
! 28: *
! 29: * Several right hand side vectors b and solution vectors x can be
! 30: * handled in a single call; they are stored as the columns of the
! 31: * M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix
! 32: * X.
! 33: *
! 34: * The effective rank of A is determined by treating as zero those
! 35: * singular values which are less than RCOND times the largest singular
! 36: * value.
! 37: *
! 38: * Arguments
! 39: * =========
! 40: *
! 41: * M (input) INTEGER
! 42: * The number of rows of the matrix A. M >= 0.
! 43: *
! 44: * N (input) INTEGER
! 45: * The number of columns of the matrix A. N >= 0.
! 46: *
! 47: * NRHS (input) INTEGER
! 48: * The number of right hand sides, i.e., the number of columns
! 49: * of the matrices B and X. NRHS >= 0.
! 50: *
! 51: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
! 52: * On entry, the M-by-N matrix A.
! 53: * On exit, the first min(m,n) rows of A are overwritten with
! 54: * its right singular vectors, stored rowwise.
! 55: *
! 56: * LDA (input) INTEGER
! 57: * The leading dimension of the array A. LDA >= max(1,M).
! 58: *
! 59: * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
! 60: * On entry, the M-by-NRHS right hand side matrix B.
! 61: * On exit, B is overwritten by the N-by-NRHS solution matrix X.
! 62: * If m >= n and RANK = n, the residual sum-of-squares for
! 63: * the solution in the i-th column is given by the sum of
! 64: * squares of the modulus of elements n+1:m in that column.
! 65: *
! 66: * LDB (input) INTEGER
! 67: * The leading dimension of the array B. LDB >= max(1,M,N).
! 68: *
! 69: * S (output) DOUBLE PRECISION array, dimension (min(M,N))
! 70: * The singular values of A in decreasing order.
! 71: * The condition number of A in the 2-norm = S(1)/S(min(m,n)).
! 72: *
! 73: * RCOND (input) DOUBLE PRECISION
! 74: * RCOND is used to determine the effective rank of A.
! 75: * Singular values S(i) <= RCOND*S(1) are treated as zero.
! 76: * If RCOND < 0, machine precision is used instead.
! 77: *
! 78: * RANK (output) INTEGER
! 79: * The effective rank of A, i.e., the number of singular values
! 80: * which are greater than RCOND*S(1).
! 81: *
! 82: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
! 83: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 84: *
! 85: * LWORK (input) INTEGER
! 86: * The dimension of the array WORK. LWORK >= 1, and also:
! 87: * LWORK >= 2*min(M,N) + max(M,N,NRHS)
! 88: * For good performance, LWORK should generally be larger.
! 89: *
! 90: * If LWORK = -1, then a workspace query is assumed; the routine
! 91: * only calculates the optimal size of the WORK array, returns
! 92: * this value as the first entry of the WORK array, and no error
! 93: * message related to LWORK is issued by XERBLA.
! 94: *
! 95: * RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,N))
! 96: *
! 97: * INFO (output) INTEGER
! 98: * = 0: successful exit
! 99: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 100: * > 0: the algorithm for computing the SVD failed to converge;
! 101: * if INFO = i, i off-diagonal elements of an intermediate
! 102: * bidiagonal form did not converge to zero.
! 103: *
! 104: * =====================================================================
! 105: *
! 106: * .. Parameters ..
! 107: DOUBLE PRECISION ZERO, ONE
! 108: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 109: COMPLEX*16 CZERO, CONE
! 110: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
! 111: $ CONE = ( 1.0D+0, 0.0D+0 ) )
! 112: * ..
! 113: * .. Local Scalars ..
! 114: LOGICAL LQUERY
! 115: INTEGER BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK,
! 116: $ ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
! 117: $ MAXWRK, MINMN, MINWRK, MM, MNTHR
! 118: DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
! 119: * ..
! 120: * .. Local Arrays ..
! 121: COMPLEX*16 VDUM( 1 )
! 122: * ..
! 123: * .. External Subroutines ..
! 124: EXTERNAL DLABAD, DLASCL, DLASET, XERBLA, ZBDSQR, ZCOPY,
! 125: $ ZDRSCL, ZGEBRD, ZGELQF, ZGEMM, ZGEMV, ZGEQRF,
! 126: $ ZLACPY, ZLASCL, ZLASET, ZUNGBR, ZUNMBR, ZUNMLQ,
! 127: $ ZUNMQR
! 128: * ..
! 129: * .. External Functions ..
! 130: INTEGER ILAENV
! 131: DOUBLE PRECISION DLAMCH, ZLANGE
! 132: EXTERNAL ILAENV, DLAMCH, ZLANGE
! 133: * ..
! 134: * .. Intrinsic Functions ..
! 135: INTRINSIC MAX, MIN
! 136: * ..
! 137: * .. Executable Statements ..
! 138: *
! 139: * Test the input arguments
! 140: *
! 141: INFO = 0
! 142: MINMN = MIN( M, N )
! 143: MAXMN = MAX( M, N )
! 144: LQUERY = ( LWORK.EQ.-1 )
! 145: IF( M.LT.0 ) THEN
! 146: INFO = -1
! 147: ELSE IF( N.LT.0 ) THEN
! 148: INFO = -2
! 149: ELSE IF( NRHS.LT.0 ) THEN
! 150: INFO = -3
! 151: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 152: INFO = -5
! 153: ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
! 154: INFO = -7
! 155: END IF
! 156: *
! 157: * Compute workspace
! 158: * (Note: Comments in the code beginning "Workspace:" describe the
! 159: * minimal amount of workspace needed at that point in the code,
! 160: * as well as the preferred amount for good performance.
! 161: * CWorkspace refers to complex workspace, and RWorkspace refers
! 162: * to real workspace. NB refers to the optimal block size for the
! 163: * immediately following subroutine, as returned by ILAENV.)
! 164: *
! 165: IF( INFO.EQ.0 ) THEN
! 166: MINWRK = 1
! 167: MAXWRK = 1
! 168: IF( MINMN.GT.0 ) THEN
! 169: MM = M
! 170: MNTHR = ILAENV( 6, 'ZGELSS', ' ', M, N, NRHS, -1 )
! 171: IF( M.GE.N .AND. M.GE.MNTHR ) THEN
! 172: *
! 173: * Path 1a - overdetermined, with many more rows than
! 174: * columns
! 175: *
! 176: MM = N
! 177: MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'ZGEQRF', ' ', M,
! 178: $ N, -1, -1 ) )
! 179: MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'ZUNMQR', 'LC',
! 180: $ M, NRHS, N, -1 ) )
! 181: END IF
! 182: IF( M.GE.N ) THEN
! 183: *
! 184: * Path 1 - overdetermined or exactly determined
! 185: *
! 186: MAXWRK = MAX( MAXWRK, 2*N + ( MM + N )*ILAENV( 1,
! 187: $ 'ZGEBRD', ' ', MM, N, -1, -1 ) )
! 188: MAXWRK = MAX( MAXWRK, 2*N + NRHS*ILAENV( 1, 'ZUNMBR',
! 189: $ 'QLC', MM, NRHS, N, -1 ) )
! 190: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
! 191: $ 'ZUNGBR', 'P', N, N, N, -1 ) )
! 192: MAXWRK = MAX( MAXWRK, N*NRHS )
! 193: MINWRK = 2*N + MAX( NRHS, M )
! 194: END IF
! 195: IF( N.GT.M ) THEN
! 196: MINWRK = 2*M + MAX( NRHS, N )
! 197: IF( N.GE.MNTHR ) THEN
! 198: *
! 199: * Path 2a - underdetermined, with many more columns
! 200: * than rows
! 201: *
! 202: MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1,
! 203: $ -1 )
! 204: MAXWRK = MAX( MAXWRK, 3*M + M*M + 2*M*ILAENV( 1,
! 205: $ 'ZGEBRD', ' ', M, M, -1, -1 ) )
! 206: MAXWRK = MAX( MAXWRK, 3*M + M*M + NRHS*ILAENV( 1,
! 207: $ 'ZUNMBR', 'QLC', M, NRHS, M, -1 ) )
! 208: MAXWRK = MAX( MAXWRK, 3*M + M*M + ( M - 1 )*ILAENV( 1,
! 209: $ 'ZUNGBR', 'P', M, M, M, -1 ) )
! 210: IF( NRHS.GT.1 ) THEN
! 211: MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
! 212: ELSE
! 213: MAXWRK = MAX( MAXWRK, M*M + 2*M )
! 214: END IF
! 215: MAXWRK = MAX( MAXWRK, M + NRHS*ILAENV( 1, 'ZUNMLQ',
! 216: $ 'LC', N, NRHS, M, -1 ) )
! 217: ELSE
! 218: *
! 219: * Path 2 - underdetermined
! 220: *
! 221: MAXWRK = 2*M + ( N + M )*ILAENV( 1, 'ZGEBRD', ' ', M,
! 222: $ N, -1, -1 )
! 223: MAXWRK = MAX( MAXWRK, 2*M + NRHS*ILAENV( 1, 'ZUNMBR',
! 224: $ 'QLC', M, NRHS, M, -1 ) )
! 225: MAXWRK = MAX( MAXWRK, 2*M + M*ILAENV( 1, 'ZUNGBR',
! 226: $ 'P', M, N, M, -1 ) )
! 227: MAXWRK = MAX( MAXWRK, N*NRHS )
! 228: END IF
! 229: END IF
! 230: MAXWRK = MAX( MINWRK, MAXWRK )
! 231: END IF
! 232: WORK( 1 ) = MAXWRK
! 233: *
! 234: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
! 235: $ INFO = -12
! 236: END IF
! 237: *
! 238: IF( INFO.NE.0 ) THEN
! 239: CALL XERBLA( 'ZGELSS', -INFO )
! 240: RETURN
! 241: ELSE IF( LQUERY ) THEN
! 242: RETURN
! 243: END IF
! 244: *
! 245: * Quick return if possible
! 246: *
! 247: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
! 248: RANK = 0
! 249: RETURN
! 250: END IF
! 251: *
! 252: * Get machine parameters
! 253: *
! 254: EPS = DLAMCH( 'P' )
! 255: SFMIN = DLAMCH( 'S' )
! 256: SMLNUM = SFMIN / EPS
! 257: BIGNUM = ONE / SMLNUM
! 258: CALL DLABAD( SMLNUM, BIGNUM )
! 259: *
! 260: * Scale A if max element outside range [SMLNUM,BIGNUM]
! 261: *
! 262: ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
! 263: IASCL = 0
! 264: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 265: *
! 266: * Scale matrix norm up to SMLNUM
! 267: *
! 268: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
! 269: IASCL = 1
! 270: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 271: *
! 272: * Scale matrix norm down to BIGNUM
! 273: *
! 274: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
! 275: IASCL = 2
! 276: ELSE IF( ANRM.EQ.ZERO ) THEN
! 277: *
! 278: * Matrix all zero. Return zero solution.
! 279: *
! 280: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
! 281: CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN )
! 282: RANK = 0
! 283: GO TO 70
! 284: END IF
! 285: *
! 286: * Scale B if max element outside range [SMLNUM,BIGNUM]
! 287: *
! 288: BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
! 289: IBSCL = 0
! 290: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
! 291: *
! 292: * Scale matrix norm up to SMLNUM
! 293: *
! 294: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
! 295: IBSCL = 1
! 296: ELSE IF( BNRM.GT.BIGNUM ) THEN
! 297: *
! 298: * Scale matrix norm down to BIGNUM
! 299: *
! 300: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
! 301: IBSCL = 2
! 302: END IF
! 303: *
! 304: * Overdetermined case
! 305: *
! 306: IF( M.GE.N ) THEN
! 307: *
! 308: * Path 1 - overdetermined or exactly determined
! 309: *
! 310: MM = M
! 311: IF( M.GE.MNTHR ) THEN
! 312: *
! 313: * Path 1a - overdetermined, with many more rows than columns
! 314: *
! 315: MM = N
! 316: ITAU = 1
! 317: IWORK = ITAU + N
! 318: *
! 319: * Compute A=Q*R
! 320: * (CWorkspace: need 2*N, prefer N+N*NB)
! 321: * (RWorkspace: none)
! 322: *
! 323: CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
! 324: $ LWORK-IWORK+1, INFO )
! 325: *
! 326: * Multiply B by transpose(Q)
! 327: * (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
! 328: * (RWorkspace: none)
! 329: *
! 330: CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
! 331: $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
! 332: *
! 333: * Zero out below R
! 334: *
! 335: IF( N.GT.1 )
! 336: $ CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
! 337: $ LDA )
! 338: END IF
! 339: *
! 340: IE = 1
! 341: ITAUQ = 1
! 342: ITAUP = ITAUQ + N
! 343: IWORK = ITAUP + N
! 344: *
! 345: * Bidiagonalize R in A
! 346: * (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
! 347: * (RWorkspace: need N)
! 348: *
! 349: CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 350: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
! 351: $ INFO )
! 352: *
! 353: * Multiply B by transpose of left bidiagonalizing vectors of R
! 354: * (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
! 355: * (RWorkspace: none)
! 356: *
! 357: CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
! 358: $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
! 359: *
! 360: * Generate right bidiagonalizing vectors of R in A
! 361: * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
! 362: * (RWorkspace: none)
! 363: *
! 364: CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
! 365: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
! 366: IRWORK = IE + N
! 367: *
! 368: * Perform bidiagonal QR iteration
! 369: * multiply B by transpose of left singular vectors
! 370: * compute right singular vectors in A
! 371: * (CWorkspace: none)
! 372: * (RWorkspace: need BDSPAC)
! 373: *
! 374: CALL ZBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM,
! 375: $ 1, B, LDB, RWORK( IRWORK ), INFO )
! 376: IF( INFO.NE.0 )
! 377: $ GO TO 70
! 378: *
! 379: * Multiply B by reciprocals of singular values
! 380: *
! 381: THR = MAX( RCOND*S( 1 ), SFMIN )
! 382: IF( RCOND.LT.ZERO )
! 383: $ THR = MAX( EPS*S( 1 ), SFMIN )
! 384: RANK = 0
! 385: DO 10 I = 1, N
! 386: IF( S( I ).GT.THR ) THEN
! 387: CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
! 388: RANK = RANK + 1
! 389: ELSE
! 390: CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
! 391: END IF
! 392: 10 CONTINUE
! 393: *
! 394: * Multiply B by right singular vectors
! 395: * (CWorkspace: need N, prefer N*NRHS)
! 396: * (RWorkspace: none)
! 397: *
! 398: IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
! 399: CALL ZGEMM( 'C', 'N', N, NRHS, N, CONE, A, LDA, B, LDB,
! 400: $ CZERO, WORK, LDB )
! 401: CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
! 402: ELSE IF( NRHS.GT.1 ) THEN
! 403: CHUNK = LWORK / N
! 404: DO 20 I = 1, NRHS, CHUNK
! 405: BL = MIN( NRHS-I+1, CHUNK )
! 406: CALL ZGEMM( 'C', 'N', N, BL, N, CONE, A, LDA, B( 1, I ),
! 407: $ LDB, CZERO, WORK, N )
! 408: CALL ZLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
! 409: 20 CONTINUE
! 410: ELSE
! 411: CALL ZGEMV( 'C', N, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
! 412: CALL ZCOPY( N, WORK, 1, B, 1 )
! 413: END IF
! 414: *
! 415: ELSE IF( N.GE.MNTHR .AND. LWORK.GE.3*M+M*M+MAX( M, NRHS, N-2*M ) )
! 416: $ THEN
! 417: *
! 418: * Underdetermined case, M much less than N
! 419: *
! 420: * Path 2a - underdetermined, with many more columns than rows
! 421: * and sufficient workspace for an efficient algorithm
! 422: *
! 423: LDWORK = M
! 424: IF( LWORK.GE.3*M+M*LDA+MAX( M, NRHS, N-2*M ) )
! 425: $ LDWORK = LDA
! 426: ITAU = 1
! 427: IWORK = M + 1
! 428: *
! 429: * Compute A=L*Q
! 430: * (CWorkspace: need 2*M, prefer M+M*NB)
! 431: * (RWorkspace: none)
! 432: *
! 433: CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
! 434: $ LWORK-IWORK+1, INFO )
! 435: IL = IWORK
! 436: *
! 437: * Copy L to WORK(IL), zeroing out above it
! 438: *
! 439: CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
! 440: CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
! 441: $ LDWORK )
! 442: IE = 1
! 443: ITAUQ = IL + LDWORK*M
! 444: ITAUP = ITAUQ + M
! 445: IWORK = ITAUP + M
! 446: *
! 447: * Bidiagonalize L in WORK(IL)
! 448: * (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
! 449: * (RWorkspace: need M)
! 450: *
! 451: CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
! 452: $ WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
! 453: $ LWORK-IWORK+1, INFO )
! 454: *
! 455: * Multiply B by transpose of left bidiagonalizing vectors of L
! 456: * (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
! 457: * (RWorkspace: none)
! 458: *
! 459: CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
! 460: $ WORK( ITAUQ ), B, LDB, WORK( IWORK ),
! 461: $ LWORK-IWORK+1, INFO )
! 462: *
! 463: * Generate right bidiagonalizing vectors of R in WORK(IL)
! 464: * (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
! 465: * (RWorkspace: none)
! 466: *
! 467: CALL ZUNGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
! 468: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
! 469: IRWORK = IE + M
! 470: *
! 471: * Perform bidiagonal QR iteration, computing right singular
! 472: * vectors of L in WORK(IL) and multiplying B by transpose of
! 473: * left singular vectors
! 474: * (CWorkspace: need M*M)
! 475: * (RWorkspace: need BDSPAC)
! 476: *
! 477: CALL ZBDSQR( 'U', M, M, 0, NRHS, S, RWORK( IE ), WORK( IL ),
! 478: $ LDWORK, A, LDA, B, LDB, RWORK( IRWORK ), INFO )
! 479: IF( INFO.NE.0 )
! 480: $ GO TO 70
! 481: *
! 482: * Multiply B by reciprocals of singular values
! 483: *
! 484: THR = MAX( RCOND*S( 1 ), SFMIN )
! 485: IF( RCOND.LT.ZERO )
! 486: $ THR = MAX( EPS*S( 1 ), SFMIN )
! 487: RANK = 0
! 488: DO 30 I = 1, M
! 489: IF( S( I ).GT.THR ) THEN
! 490: CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
! 491: RANK = RANK + 1
! 492: ELSE
! 493: CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
! 494: END IF
! 495: 30 CONTINUE
! 496: IWORK = IL + M*LDWORK
! 497: *
! 498: * Multiply B by right singular vectors of L in WORK(IL)
! 499: * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
! 500: * (RWorkspace: none)
! 501: *
! 502: IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
! 503: CALL ZGEMM( 'C', 'N', M, NRHS, M, CONE, WORK( IL ), LDWORK,
! 504: $ B, LDB, CZERO, WORK( IWORK ), LDB )
! 505: CALL ZLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
! 506: ELSE IF( NRHS.GT.1 ) THEN
! 507: CHUNK = ( LWORK-IWORK+1 ) / M
! 508: DO 40 I = 1, NRHS, CHUNK
! 509: BL = MIN( NRHS-I+1, CHUNK )
! 510: CALL ZGEMM( 'C', 'N', M, BL, M, CONE, WORK( IL ), LDWORK,
! 511: $ B( 1, I ), LDB, CZERO, WORK( IWORK ), M )
! 512: CALL ZLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
! 513: $ LDB )
! 514: 40 CONTINUE
! 515: ELSE
! 516: CALL ZGEMV( 'C', M, M, CONE, WORK( IL ), LDWORK, B( 1, 1 ),
! 517: $ 1, CZERO, WORK( IWORK ), 1 )
! 518: CALL ZCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
! 519: END IF
! 520: *
! 521: * Zero out below first M rows of B
! 522: *
! 523: CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
! 524: IWORK = ITAU + M
! 525: *
! 526: * Multiply transpose(Q) by B
! 527: * (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
! 528: * (RWorkspace: none)
! 529: *
! 530: CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
! 531: $ LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
! 532: *
! 533: ELSE
! 534: *
! 535: * Path 2 - remaining underdetermined cases
! 536: *
! 537: IE = 1
! 538: ITAUQ = 1
! 539: ITAUP = ITAUQ + M
! 540: IWORK = ITAUP + M
! 541: *
! 542: * Bidiagonalize A
! 543: * (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
! 544: * (RWorkspace: need N)
! 545: *
! 546: CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
! 547: $ WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
! 548: $ INFO )
! 549: *
! 550: * Multiply B by transpose of left bidiagonalizing vectors
! 551: * (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
! 552: * (RWorkspace: none)
! 553: *
! 554: CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
! 555: $ B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
! 556: *
! 557: * Generate right bidiagonalizing vectors in A
! 558: * (CWorkspace: need 3*M, prefer 2*M+M*NB)
! 559: * (RWorkspace: none)
! 560: *
! 561: CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
! 562: $ WORK( IWORK ), LWORK-IWORK+1, INFO )
! 563: IRWORK = IE + M
! 564: *
! 565: * Perform bidiagonal QR iteration,
! 566: * computing right singular vectors of A in A and
! 567: * multiplying B by transpose of left singular vectors
! 568: * (CWorkspace: none)
! 569: * (RWorkspace: need BDSPAC)
! 570: *
! 571: CALL ZBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, VDUM,
! 572: $ 1, B, LDB, RWORK( IRWORK ), INFO )
! 573: IF( INFO.NE.0 )
! 574: $ GO TO 70
! 575: *
! 576: * Multiply B by reciprocals of singular values
! 577: *
! 578: THR = MAX( RCOND*S( 1 ), SFMIN )
! 579: IF( RCOND.LT.ZERO )
! 580: $ THR = MAX( EPS*S( 1 ), SFMIN )
! 581: RANK = 0
! 582: DO 50 I = 1, M
! 583: IF( S( I ).GT.THR ) THEN
! 584: CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
! 585: RANK = RANK + 1
! 586: ELSE
! 587: CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
! 588: END IF
! 589: 50 CONTINUE
! 590: *
! 591: * Multiply B by right singular vectors of A
! 592: * (CWorkspace: need N, prefer N*NRHS)
! 593: * (RWorkspace: none)
! 594: *
! 595: IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
! 596: CALL ZGEMM( 'C', 'N', N, NRHS, M, CONE, A, LDA, B, LDB,
! 597: $ CZERO, WORK, LDB )
! 598: CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
! 599: ELSE IF( NRHS.GT.1 ) THEN
! 600: CHUNK = LWORK / N
! 601: DO 60 I = 1, NRHS, CHUNK
! 602: BL = MIN( NRHS-I+1, CHUNK )
! 603: CALL ZGEMM( 'C', 'N', N, BL, M, CONE, A, LDA, B( 1, I ),
! 604: $ LDB, CZERO, WORK, N )
! 605: CALL ZLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
! 606: 60 CONTINUE
! 607: ELSE
! 608: CALL ZGEMV( 'C', M, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
! 609: CALL ZCOPY( N, WORK, 1, B, 1 )
! 610: END IF
! 611: END IF
! 612: *
! 613: * Undo scaling
! 614: *
! 615: IF( IASCL.EQ.1 ) THEN
! 616: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
! 617: CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
! 618: $ INFO )
! 619: ELSE IF( IASCL.EQ.2 ) THEN
! 620: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
! 621: CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
! 622: $ INFO )
! 623: END IF
! 624: IF( IBSCL.EQ.1 ) THEN
! 625: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
! 626: ELSE IF( IBSCL.EQ.2 ) THEN
! 627: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
! 628: END IF
! 629: 70 CONTINUE
! 630: WORK( 1 ) = MAXWRK
! 631: RETURN
! 632: *
! 633: * End of ZGELSS
! 634: *
! 635: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgelss.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack