Annotation of rpl/lapack/lapack/zgels.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
! 2: $ INFO )
! 3: *
! 4: * -- LAPACK driver routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER TRANS
! 11: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
! 12: * ..
! 13: * .. Array Arguments ..
! 14: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
! 15: * ..
! 16: *
! 17: * Purpose
! 18: * =======
! 19: *
! 20: * ZGELS solves overdetermined or underdetermined complex linear systems
! 21: * involving an M-by-N matrix A, or its conjugate-transpose, using a QR
! 22: * or LQ factorization of A. It is assumed that A has full rank.
! 23: *
! 24: * The following options are provided:
! 25: *
! 26: * 1. If TRANS = 'N' and m >= n: find the least squares solution of
! 27: * an overdetermined system, i.e., solve the least squares problem
! 28: * minimize || B - A*X ||.
! 29: *
! 30: * 2. If TRANS = 'N' and m < n: find the minimum norm solution of
! 31: * an underdetermined system A * X = B.
! 32: *
! 33: * 3. If TRANS = 'C' and m >= n: find the minimum norm solution of
! 34: * an undetermined system A**H * X = B.
! 35: *
! 36: * 4. If TRANS = 'C' and m < n: find the least squares solution of
! 37: * an overdetermined system, i.e., solve the least squares problem
! 38: * minimize || B - A**H * X ||.
! 39: *
! 40: * Several right hand side vectors b and solution vectors x can be
! 41: * handled in a single call; they are stored as the columns of the
! 42: * M-by-NRHS right hand side matrix B and the N-by-NRHS solution
! 43: * matrix X.
! 44: *
! 45: * Arguments
! 46: * =========
! 47: *
! 48: * TRANS (input) CHARACTER*1
! 49: * = 'N': the linear system involves A;
! 50: * = 'C': the linear system involves A**H.
! 51: *
! 52: * M (input) INTEGER
! 53: * The number of rows of the matrix A. M >= 0.
! 54: *
! 55: * N (input) INTEGER
! 56: * The number of columns of the matrix A. N >= 0.
! 57: *
! 58: * NRHS (input) INTEGER
! 59: * The number of right hand sides, i.e., the number of
! 60: * columns of the matrices B and X. NRHS >= 0.
! 61: *
! 62: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
! 63: * On entry, the M-by-N matrix A.
! 64: * if M >= N, A is overwritten by details of its QR
! 65: * factorization as returned by ZGEQRF;
! 66: * if M < N, A is overwritten by details of its LQ
! 67: * factorization as returned by ZGELQF.
! 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 matrix B of right hand side vectors, stored
! 74: * columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
! 75: * if TRANS = 'C'.
! 76: * On exit, if INFO = 0, B is overwritten by the solution
! 77: * vectors, stored columnwise:
! 78: * if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
! 79: * squares solution vectors; the residual sum of squares for the
! 80: * solution in each column is given by the sum of squares of the
! 81: * modulus of elements N+1 to M in that column;
! 82: * if TRANS = 'N' and m < n, rows 1 to N of B contain the
! 83: * minimum norm solution vectors;
! 84: * if TRANS = 'C' and m >= n, rows 1 to M of B contain the
! 85: * minimum norm solution vectors;
! 86: * if TRANS = 'C' and m < n, rows 1 to M of B contain the
! 87: * least squares solution vectors; the residual sum of squares
! 88: * for the solution in each column is given by the sum of
! 89: * squares of the modulus of elements M+1 to N in that column.
! 90: *
! 91: * LDB (input) INTEGER
! 92: * The leading dimension of the array B. LDB >= MAX(1,M,N).
! 93: *
! 94: * WORK (workspace/output) COMPLEX*16 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.
! 99: * LWORK >= max( 1, MN + max( MN, NRHS ) ).
! 100: * For optimal performance,
! 101: * LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
! 102: * where MN = min(M,N) and NB is the optimum block size.
! 103: *
! 104: * If LWORK = -1, then a workspace query is assumed; the routine
! 105: * only calculates the optimal size of the WORK array, returns
! 106: * this value as the first entry of the WORK array, and no error
! 107: * message related to LWORK is issued by XERBLA.
! 108: *
! 109: * INFO (output) INTEGER
! 110: * = 0: successful exit
! 111: * < 0: if INFO = -i, the i-th argument had an illegal value
! 112: * > 0: if INFO = i, the i-th diagonal element of the
! 113: * triangular factor of A is zero, so that A does not have
! 114: * full rank; the least squares solution could not be
! 115: * computed.
! 116: *
! 117: * =====================================================================
! 118: *
! 119: * .. Parameters ..
! 120: DOUBLE PRECISION ZERO, ONE
! 121: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 122: COMPLEX*16 CZERO
! 123: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
! 124: * ..
! 125: * .. Local Scalars ..
! 126: LOGICAL LQUERY, TPSD
! 127: INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
! 128: DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
! 129: * ..
! 130: * .. Local Arrays ..
! 131: DOUBLE PRECISION RWORK( 1 )
! 132: * ..
! 133: * .. External Functions ..
! 134: LOGICAL LSAME
! 135: INTEGER ILAENV
! 136: DOUBLE PRECISION DLAMCH, ZLANGE
! 137: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
! 138: * ..
! 139: * .. External Subroutines ..
! 140: EXTERNAL DLABAD, XERBLA, ZGELQF, ZGEQRF, ZLASCL, ZLASET,
! 141: $ ZTRTRS, ZUNMLQ, ZUNMQR
! 142: * ..
! 143: * .. Intrinsic Functions ..
! 144: INTRINSIC DBLE, MAX, MIN
! 145: * ..
! 146: * .. Executable Statements ..
! 147: *
! 148: * Test the input arguments.
! 149: *
! 150: INFO = 0
! 151: MN = MIN( M, N )
! 152: LQUERY = ( LWORK.EQ.-1 )
! 153: IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'C' ) ) ) THEN
! 154: INFO = -1
! 155: ELSE IF( M.LT.0 ) THEN
! 156: INFO = -2
! 157: ELSE IF( N.LT.0 ) THEN
! 158: INFO = -3
! 159: ELSE IF( NRHS.LT.0 ) THEN
! 160: INFO = -4
! 161: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 162: INFO = -6
! 163: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
! 164: INFO = -8
! 165: ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY )
! 166: $ THEN
! 167: INFO = -10
! 168: END IF
! 169: *
! 170: * Figure out optimal block size
! 171: *
! 172: IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
! 173: *
! 174: TPSD = .TRUE.
! 175: IF( LSAME( TRANS, 'N' ) )
! 176: $ TPSD = .FALSE.
! 177: *
! 178: IF( M.GE.N ) THEN
! 179: NB = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 180: IF( TPSD ) THEN
! 181: NB = MAX( NB, ILAENV( 1, 'ZUNMQR', 'LN', M, NRHS, N,
! 182: $ -1 ) )
! 183: ELSE
! 184: NB = MAX( NB, ILAENV( 1, 'ZUNMQR', 'LC', M, NRHS, N,
! 185: $ -1 ) )
! 186: END IF
! 187: ELSE
! 188: NB = ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
! 189: IF( TPSD ) THEN
! 190: NB = MAX( NB, ILAENV( 1, 'ZUNMLQ', 'LC', N, NRHS, M,
! 191: $ -1 ) )
! 192: ELSE
! 193: NB = MAX( NB, ILAENV( 1, 'ZUNMLQ', 'LN', N, NRHS, M,
! 194: $ -1 ) )
! 195: END IF
! 196: END IF
! 197: *
! 198: WSIZE = MAX( 1, MN+MAX( MN, NRHS )*NB )
! 199: WORK( 1 ) = DBLE( WSIZE )
! 200: *
! 201: END IF
! 202: *
! 203: IF( INFO.NE.0 ) THEN
! 204: CALL XERBLA( 'ZGELS ', -INFO )
! 205: RETURN
! 206: ELSE IF( LQUERY ) THEN
! 207: RETURN
! 208: END IF
! 209: *
! 210: * Quick return if possible
! 211: *
! 212: IF( MIN( M, N, NRHS ).EQ.0 ) THEN
! 213: CALL ZLASET( 'Full', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
! 214: RETURN
! 215: END IF
! 216: *
! 217: * Get machine parameters
! 218: *
! 219: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
! 220: BIGNUM = ONE / SMLNUM
! 221: CALL DLABAD( SMLNUM, BIGNUM )
! 222: *
! 223: * Scale A, B if max element outside range [SMLNUM,BIGNUM]
! 224: *
! 225: ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
! 226: IASCL = 0
! 227: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 228: *
! 229: * Scale matrix norm up to SMLNUM
! 230: *
! 231: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
! 232: IASCL = 1
! 233: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 234: *
! 235: * Scale matrix norm down to BIGNUM
! 236: *
! 237: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
! 238: IASCL = 2
! 239: ELSE IF( ANRM.EQ.ZERO ) THEN
! 240: *
! 241: * Matrix all zero. Return zero solution.
! 242: *
! 243: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
! 244: GO TO 50
! 245: END IF
! 246: *
! 247: BROW = M
! 248: IF( TPSD )
! 249: $ BROW = N
! 250: BNRM = ZLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
! 251: IBSCL = 0
! 252: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
! 253: *
! 254: * Scale matrix norm up to SMLNUM
! 255: *
! 256: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
! 257: $ INFO )
! 258: IBSCL = 1
! 259: ELSE IF( BNRM.GT.BIGNUM ) THEN
! 260: *
! 261: * Scale matrix norm down to BIGNUM
! 262: *
! 263: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
! 264: $ INFO )
! 265: IBSCL = 2
! 266: END IF
! 267: *
! 268: IF( M.GE.N ) THEN
! 269: *
! 270: * compute QR factorization of A
! 271: *
! 272: CALL ZGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
! 273: $ INFO )
! 274: *
! 275: * workspace at least N, optimally N*NB
! 276: *
! 277: IF( .NOT.TPSD ) THEN
! 278: *
! 279: * Least-Squares Problem min || A * X - B ||
! 280: *
! 281: * B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
! 282: *
! 283: CALL ZUNMQR( 'Left', 'Conjugate transpose', M, NRHS, N, A,
! 284: $ LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
! 285: $ INFO )
! 286: *
! 287: * workspace at least NRHS, optimally NRHS*NB
! 288: *
! 289: * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
! 290: *
! 291: CALL ZTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS,
! 292: $ A, LDA, B, LDB, INFO )
! 293: *
! 294: IF( INFO.GT.0 ) THEN
! 295: RETURN
! 296: END IF
! 297: *
! 298: SCLLEN = N
! 299: *
! 300: ELSE
! 301: *
! 302: * Overdetermined system of equations A' * X = B
! 303: *
! 304: * B(1:N,1:NRHS) := inv(R') * B(1:N,1:NRHS)
! 305: *
! 306: CALL ZTRTRS( 'Upper', 'Conjugate transpose','Non-unit',
! 307: $ N, NRHS, A, LDA, B, LDB, INFO )
! 308: *
! 309: IF( INFO.GT.0 ) THEN
! 310: RETURN
! 311: END IF
! 312: *
! 313: * B(N+1:M,1:NRHS) = ZERO
! 314: *
! 315: DO 20 J = 1, NRHS
! 316: DO 10 I = N + 1, M
! 317: B( I, J ) = CZERO
! 318: 10 CONTINUE
! 319: 20 CONTINUE
! 320: *
! 321: * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
! 322: *
! 323: CALL ZUNMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA,
! 324: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
! 325: $ INFO )
! 326: *
! 327: * workspace at least NRHS, optimally NRHS*NB
! 328: *
! 329: SCLLEN = M
! 330: *
! 331: END IF
! 332: *
! 333: ELSE
! 334: *
! 335: * Compute LQ factorization of A
! 336: *
! 337: CALL ZGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
! 338: $ INFO )
! 339: *
! 340: * workspace at least M, optimally M*NB.
! 341: *
! 342: IF( .NOT.TPSD ) THEN
! 343: *
! 344: * underdetermined system of equations A * X = B
! 345: *
! 346: * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
! 347: *
! 348: CALL ZTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS,
! 349: $ A, LDA, B, LDB, INFO )
! 350: *
! 351: IF( INFO.GT.0 ) THEN
! 352: RETURN
! 353: END IF
! 354: *
! 355: * B(M+1:N,1:NRHS) = 0
! 356: *
! 357: DO 40 J = 1, NRHS
! 358: DO 30 I = M + 1, N
! 359: B( I, J ) = CZERO
! 360: 30 CONTINUE
! 361: 40 CONTINUE
! 362: *
! 363: * B(1:N,1:NRHS) := Q(1:N,:)' * B(1:M,1:NRHS)
! 364: *
! 365: CALL ZUNMLQ( 'Left', 'Conjugate transpose', N, NRHS, M, A,
! 366: $ LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
! 367: $ INFO )
! 368: *
! 369: * workspace at least NRHS, optimally NRHS*NB
! 370: *
! 371: SCLLEN = N
! 372: *
! 373: ELSE
! 374: *
! 375: * overdetermined system min || A' * X - B ||
! 376: *
! 377: * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
! 378: *
! 379: CALL ZUNMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA,
! 380: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
! 381: $ INFO )
! 382: *
! 383: * workspace at least NRHS, optimally NRHS*NB
! 384: *
! 385: * B(1:M,1:NRHS) := inv(L') * B(1:M,1:NRHS)
! 386: *
! 387: CALL ZTRTRS( 'Lower', 'Conjugate transpose', 'Non-unit',
! 388: $ M, NRHS, A, LDA, B, LDB, INFO )
! 389: *
! 390: IF( INFO.GT.0 ) THEN
! 391: RETURN
! 392: END IF
! 393: *
! 394: SCLLEN = M
! 395: *
! 396: END IF
! 397: *
! 398: END IF
! 399: *
! 400: * Undo scaling
! 401: *
! 402: IF( IASCL.EQ.1 ) THEN
! 403: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
! 404: $ INFO )
! 405: ELSE IF( IASCL.EQ.2 ) THEN
! 406: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
! 407: $ INFO )
! 408: END IF
! 409: IF( IBSCL.EQ.1 ) THEN
! 410: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
! 411: $ INFO )
! 412: ELSE IF( IBSCL.EQ.2 ) THEN
! 413: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
! 414: $ INFO )
! 415: END IF
! 416: *
! 417: 50 CONTINUE
! 418: WORK( 1 ) = DBLE( WSIZE )
! 419: *
! 420: RETURN
! 421: *
! 422: * End of ZGELS
! 423: *
! 424: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgels.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack