![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, 2: $ WORK, 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, M, N, NRHS, RANK 11: DOUBLE PRECISION RCOND 12: * .. 13: * .. Array Arguments .. 14: INTEGER JPVT( * ) 15: DOUBLE PRECISION RWORK( * ) 16: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) 17: * .. 18: * 19: * Purpose 20: * ======= 21: * 22: * This routine is deprecated and has been replaced by routine ZGELSY. 23: * 24: * ZGELSX computes the minimum-norm solution to a complex linear least 25: * squares problem: 26: * minimize || A * X - B || 27: * using a complete orthogonal factorization of A. A is an M-by-N 28: * matrix which may be rank-deficient. 29: * 30: * Several right hand side vectors b and solution vectors x can be 31: * handled in a single call; they are stored as the columns of the 32: * M-by-NRHS right hand side matrix B and the N-by-NRHS solution 33: * matrix X. 34: * 35: * The routine first computes a QR factorization with column pivoting: 36: * A * P = Q * [ R11 R12 ] 37: * [ 0 R22 ] 38: * with R11 defined as the largest leading submatrix whose estimated 39: * condition number is less than 1/RCOND. The order of R11, RANK, 40: * is the effective rank of A. 41: * 42: * Then, R22 is considered to be negligible, and R12 is annihilated 43: * by unitary transformations from the right, arriving at the 44: * complete orthogonal factorization: 45: * A * P = Q * [ T11 0 ] * Z 46: * [ 0 0 ] 47: * The minimum-norm solution is then 48: * X = P * Z' [ inv(T11)*Q1'*B ] 49: * [ 0 ] 50: * where Q1 consists of the first RANK columns of Q. 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 63: * columns of matrices B and X. NRHS >= 0. 64: * 65: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 66: * On entry, the M-by-N matrix A. 67: * On exit, A has been overwritten by details of its 68: * complete orthogonal factorization. 69: * 70: * LDA (input) INTEGER 71: * The leading dimension of the array A. LDA >= max(1,M). 72: * 73: * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) 74: * On entry, the M-by-NRHS right hand side matrix B. 75: * On exit, the N-by-NRHS solution matrix X. 76: * If m >= n and RANK = n, the residual sum-of-squares for 77: * the solution in the i-th column is given by the sum of 78: * squares of elements N+1:M in that column. 79: * 80: * LDB (input) INTEGER 81: * The leading dimension of the array B. LDB >= max(1,M,N). 82: * 83: * JPVT (input/output) INTEGER array, dimension (N) 84: * On entry, if JPVT(i) .ne. 0, the i-th column of A is an 85: * initial column, otherwise it is a free column. Before 86: * the QR factorization of A, all initial columns are 87: * permuted to the leading positions; only the remaining 88: * free columns are moved as a result of column pivoting 89: * during the factorization. 90: * On exit, if JPVT(i) = k, then the i-th column of A*P 91: * was the k-th column of A. 92: * 93: * RCOND (input) DOUBLE PRECISION 94: * RCOND is used to determine the effective rank of A, which 95: * is defined as the order of the largest leading triangular 96: * submatrix R11 in the QR factorization with pivoting of A, 97: * whose estimated condition number < 1/RCOND. 98: * 99: * RANK (output) INTEGER 100: * The effective rank of A, i.e., the order of the submatrix 101: * R11. This is the same as the order of the submatrix T11 102: * in the complete orthogonal factorization of A. 103: * 104: * WORK (workspace) COMPLEX*16 array, dimension 105: * (min(M,N) + max( N, 2*min(M,N)+NRHS )), 106: * 107: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) 108: * 109: * INFO (output) INTEGER 110: * = 0: successful exit 111: * < 0: if INFO = -i, the i-th argument had an illegal value 112: * 113: * ===================================================================== 114: * 115: * .. Parameters .. 116: INTEGER IMAX, IMIN 117: PARAMETER ( IMAX = 1, IMIN = 2 ) 118: DOUBLE PRECISION ZERO, ONE, DONE, NTDONE 119: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, DONE = ZERO, 120: $ NTDONE = ONE ) 121: COMPLEX*16 CZERO, CONE 122: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 123: $ CONE = ( 1.0D+0, 0.0D+0 ) ) 124: * .. 125: * .. Local Scalars .. 126: INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN 127: DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR, 128: $ SMLNUM 129: COMPLEX*16 C1, C2, S1, S2, T1, T2 130: * .. 131: * .. External Subroutines .. 132: EXTERNAL XERBLA, ZGEQPF, ZLAIC1, ZLASCL, ZLASET, ZLATZM, 133: $ ZTRSM, ZTZRQF, ZUNM2R 134: * .. 135: * .. External Functions .. 136: DOUBLE PRECISION DLAMCH, ZLANGE 137: EXTERNAL DLAMCH, ZLANGE 138: * .. 139: * .. Intrinsic Functions .. 140: INTRINSIC ABS, DCONJG, MAX, MIN 141: * .. 142: * .. Executable Statements .. 143: * 144: MN = MIN( M, N ) 145: ISMIN = MN + 1 146: ISMAX = 2*MN + 1 147: * 148: * Test the input arguments. 149: * 150: INFO = 0 151: IF( M.LT.0 ) THEN 152: INFO = -1 153: ELSE IF( N.LT.0 ) THEN 154: INFO = -2 155: ELSE IF( NRHS.LT.0 ) THEN 156: INFO = -3 157: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 158: INFO = -5 159: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN 160: INFO = -7 161: END IF 162: * 163: IF( INFO.NE.0 ) THEN 164: CALL XERBLA( 'ZGELSX', -INFO ) 165: RETURN 166: END IF 167: * 168: * Quick return if possible 169: * 170: IF( MIN( M, N, NRHS ).EQ.0 ) THEN 171: RANK = 0 172: RETURN 173: END IF 174: * 175: * Get machine parameters 176: * 177: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) 178: BIGNUM = ONE / SMLNUM 179: CALL DLABAD( SMLNUM, BIGNUM ) 180: * 181: * Scale A, B if max elements outside range [SMLNUM,BIGNUM] 182: * 183: ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK ) 184: IASCL = 0 185: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 186: * 187: * Scale matrix norm up to SMLNUM 188: * 189: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 190: IASCL = 1 191: ELSE IF( ANRM.GT.BIGNUM ) THEN 192: * 193: * Scale matrix norm down to BIGNUM 194: * 195: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 196: IASCL = 2 197: ELSE IF( ANRM.EQ.ZERO ) THEN 198: * 199: * Matrix all zero. Return zero solution. 200: * 201: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 202: RANK = 0 203: GO TO 100 204: END IF 205: * 206: BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK ) 207: IBSCL = 0 208: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 209: * 210: * Scale matrix norm up to SMLNUM 211: * 212: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 213: IBSCL = 1 214: ELSE IF( BNRM.GT.BIGNUM ) THEN 215: * 216: * Scale matrix norm down to BIGNUM 217: * 218: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 219: IBSCL = 2 220: END IF 221: * 222: * Compute QR factorization with column pivoting of A: 223: * A * P = Q * R 224: * 225: CALL ZGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), RWORK, 226: $ INFO ) 227: * 228: * complex workspace MN+N. Real workspace 2*N. Details of Householder 229: * rotations stored in WORK(1:MN). 230: * 231: * Determine RANK using incremental condition estimation 232: * 233: WORK( ISMIN ) = CONE 234: WORK( ISMAX ) = CONE 235: SMAX = ABS( A( 1, 1 ) ) 236: SMIN = SMAX 237: IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN 238: RANK = 0 239: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 240: GO TO 100 241: ELSE 242: RANK = 1 243: END IF 244: * 245: 10 CONTINUE 246: IF( RANK.LT.MN ) THEN 247: I = RANK + 1 248: CALL ZLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ), 249: $ A( I, I ), SMINPR, S1, C1 ) 250: CALL ZLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ), 251: $ A( I, I ), SMAXPR, S2, C2 ) 252: * 253: IF( SMAXPR*RCOND.LE.SMINPR ) THEN 254: DO 20 I = 1, RANK 255: WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 ) 256: WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 ) 257: 20 CONTINUE 258: WORK( ISMIN+RANK ) = C1 259: WORK( ISMAX+RANK ) = C2 260: SMIN = SMINPR 261: SMAX = SMAXPR 262: RANK = RANK + 1 263: GO TO 10 264: END IF 265: END IF 266: * 267: * Logically partition R = [ R11 R12 ] 268: * [ 0 R22 ] 269: * where R11 = R(1:RANK,1:RANK) 270: * 271: * [R11,R12] = [ T11, 0 ] * Y 272: * 273: IF( RANK.LT.N ) 274: $ CALL ZTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO ) 275: * 276: * Details of Householder rotations stored in WORK(MN+1:2*MN) 277: * 278: * B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) 279: * 280: CALL ZUNM2R( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA, 281: $ WORK( 1 ), B, LDB, WORK( 2*MN+1 ), INFO ) 282: * 283: * workspace NRHS 284: * 285: * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) 286: * 287: CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK, 288: $ NRHS, CONE, A, LDA, B, LDB ) 289: * 290: DO 40 I = RANK + 1, N 291: DO 30 J = 1, NRHS 292: B( I, J ) = CZERO 293: 30 CONTINUE 294: 40 CONTINUE 295: * 296: * B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) 297: * 298: IF( RANK.LT.N ) THEN 299: DO 50 I = 1, RANK 300: CALL ZLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA, 301: $ DCONJG( WORK( MN+I ) ), B( I, 1 ), 302: $ B( RANK+1, 1 ), LDB, WORK( 2*MN+1 ) ) 303: 50 CONTINUE 304: END IF 305: * 306: * workspace NRHS 307: * 308: * B(1:N,1:NRHS) := P * B(1:N,1:NRHS) 309: * 310: DO 90 J = 1, NRHS 311: DO 60 I = 1, N 312: WORK( 2*MN+I ) = NTDONE 313: 60 CONTINUE 314: DO 80 I = 1, N 315: IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN 316: IF( JPVT( I ).NE.I ) THEN 317: K = I 318: T1 = B( K, J ) 319: T2 = B( JPVT( K ), J ) 320: 70 CONTINUE 321: B( JPVT( K ), J ) = T1 322: WORK( 2*MN+K ) = DONE 323: T1 = T2 324: K = JPVT( K ) 325: T2 = B( JPVT( K ), J ) 326: IF( JPVT( K ).NE.I ) 327: $ GO TO 70 328: B( I, J ) = T1 329: WORK( 2*MN+K ) = DONE 330: END IF 331: END IF 332: 80 CONTINUE 333: 90 CONTINUE 334: * 335: * Undo scaling 336: * 337: IF( IASCL.EQ.1 ) THEN 338: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 339: CALL ZLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA, 340: $ INFO ) 341: ELSE IF( IASCL.EQ.2 ) THEN 342: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 343: CALL ZLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA, 344: $ INFO ) 345: END IF 346: IF( IBSCL.EQ.1 ) THEN 347: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 348: ELSE IF( IBSCL.EQ.2 ) THEN 349: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 350: END IF 351: * 352: 100 CONTINUE 353: * 354: RETURN 355: * 356: * End of ZGELSX 357: * 358: END