![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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