![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, 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: INTEGER JPVT( * ) 15: DOUBLE PRECISION RWORK( * ) 16: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) 17: * .. 18: * 19: * Purpose 20: * ======= 21: * 22: * ZGELSY computes the minimum-norm solution to a complex linear least 23: * squares problem: 24: * minimize || A * X - B || 25: * using a complete orthogonal factorization 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 routine first computes a QR factorization with column pivoting: 34: * A * P = Q * [ R11 R12 ] 35: * [ 0 R22 ] 36: * with R11 defined as the largest leading submatrix whose estimated 37: * condition number is less than 1/RCOND. The order of R11, RANK, 38: * is the effective rank of A. 39: * 40: * Then, R22 is considered to be negligible, and R12 is annihilated 41: * by unitary transformations from the right, arriving at the 42: * complete orthogonal factorization: 43: * A * P = Q * [ T11 0 ] * Z 44: * [ 0 0 ] 45: * The minimum-norm solution is then 46: * X = P * Z' [ inv(T11)*Q1'*B ] 47: * [ 0 ] 48: * where Q1 consists of the first RANK columns of Q. 49: * 50: * This routine is basically identical to the original xGELSX except 51: * three differences: 52: * o The permutation of matrix B (the right hand side) is faster and 53: * more simple. 54: * o The call to the subroutine xGEQPF has been substituted by the 55: * the call to the subroutine xGEQP3. This subroutine is a Blas-3 56: * version of the QR factorization with column pivoting. 57: * o Matrix B (the right hand side) is updated with Blas-3. 58: * 59: * Arguments 60: * ========= 61: * 62: * M (input) INTEGER 63: * The number of rows of the matrix A. M >= 0. 64: * 65: * N (input) INTEGER 66: * The number of columns of the matrix A. N >= 0. 67: * 68: * NRHS (input) INTEGER 69: * The number of right hand sides, i.e., the number of 70: * columns of matrices B and X. NRHS >= 0. 71: * 72: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 73: * On entry, the M-by-N matrix A. 74: * On exit, A has been overwritten by details of its 75: * complete orthogonal factorization. 76: * 77: * LDA (input) INTEGER 78: * The leading dimension of the array A. LDA >= max(1,M). 79: * 80: * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) 81: * On entry, the M-by-NRHS right hand side matrix B. 82: * On exit, the N-by-NRHS solution matrix X. 83: * 84: * LDB (input) INTEGER 85: * The leading dimension of the array B. LDB >= max(1,M,N). 86: * 87: * JPVT (input/output) INTEGER array, dimension (N) 88: * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted 89: * to the front of AP, otherwise column i is a free column. 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/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 105: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 106: * 107: * LWORK (input) INTEGER 108: * The dimension of the array WORK. 109: * The unblocked strategy requires that: 110: * LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS ) 111: * where MN = min(M,N). 112: * The block algorithm requires that: 113: * LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS ) 114: * where NB is an upper bound on the blocksize returned 115: * by ILAENV for the routines ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR, 116: * and ZUNMRZ. 117: * 118: * If LWORK = -1, then a workspace query is assumed; the routine 119: * only calculates the optimal size of the WORK array, returns 120: * this value as the first entry of the WORK array, and no error 121: * message related to LWORK is issued by XERBLA. 122: * 123: * RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) 124: * 125: * INFO (output) INTEGER 126: * = 0: successful exit 127: * < 0: if INFO = -i, the i-th argument had an illegal value 128: * 129: * Further Details 130: * =============== 131: * 132: * Based on contributions by 133: * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA 134: * E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain 135: * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain 136: * 137: * ===================================================================== 138: * 139: * .. Parameters .. 140: INTEGER IMAX, IMIN 141: PARAMETER ( IMAX = 1, IMIN = 2 ) 142: DOUBLE PRECISION ZERO, ONE 143: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 144: COMPLEX*16 CZERO, CONE 145: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 146: $ CONE = ( 1.0D+0, 0.0D+0 ) ) 147: * .. 148: * .. Local Scalars .. 149: LOGICAL LQUERY 150: INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN, 151: $ NB, NB1, NB2, NB3, NB4 152: DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR, 153: $ SMLNUM, WSIZE 154: COMPLEX*16 C1, C2, S1, S2 155: * .. 156: * .. External Subroutines .. 157: EXTERNAL DLABAD, XERBLA, ZCOPY, ZGEQP3, ZLAIC1, ZLASCL, 158: $ ZLASET, ZTRSM, ZTZRZF, ZUNMQR, ZUNMRZ 159: * .. 160: * .. External Functions .. 161: INTEGER ILAENV 162: DOUBLE PRECISION DLAMCH, ZLANGE 163: EXTERNAL ILAENV, DLAMCH, ZLANGE 164: * .. 165: * .. Intrinsic Functions .. 166: INTRINSIC ABS, DBLE, DCMPLX, MAX, MIN 167: * .. 168: * .. Executable Statements .. 169: * 170: MN = MIN( M, N ) 171: ISMIN = MN + 1 172: ISMAX = 2*MN + 1 173: * 174: * Test the input arguments. 175: * 176: INFO = 0 177: NB1 = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 ) 178: NB2 = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 ) 179: NB3 = ILAENV( 1, 'ZUNMQR', ' ', M, N, NRHS, -1 ) 180: NB4 = ILAENV( 1, 'ZUNMRQ', ' ', M, N, NRHS, -1 ) 181: NB = MAX( NB1, NB2, NB3, NB4 ) 182: LWKOPT = MAX( 1, MN+2*N+NB*( N+1 ), 2*MN+NB*NRHS ) 183: WORK( 1 ) = DCMPLX( LWKOPT ) 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, M, N ) ) THEN 194: INFO = -7 195: ELSE IF( LWORK.LT.( MN+MAX( 2*MN, N+1, MN+NRHS ) ) .AND. .NOT. 196: $ LQUERY ) THEN 197: INFO = -12 198: END IF 199: * 200: IF( INFO.NE.0 ) THEN 201: CALL XERBLA( 'ZGELSY', -INFO ) 202: RETURN 203: ELSE IF( LQUERY ) THEN 204: RETURN 205: END IF 206: * 207: * Quick return if possible 208: * 209: IF( MIN( M, N, NRHS ).EQ.0 ) THEN 210: RANK = 0 211: RETURN 212: END IF 213: * 214: * Get machine parameters 215: * 216: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) 217: BIGNUM = ONE / SMLNUM 218: CALL DLABAD( SMLNUM, BIGNUM ) 219: * 220: * Scale A, B if max entries outside range [SMLNUM,BIGNUM] 221: * 222: ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK ) 223: IASCL = 0 224: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 225: * 226: * Scale matrix norm up to SMLNUM 227: * 228: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO ) 229: IASCL = 1 230: ELSE IF( ANRM.GT.BIGNUM ) THEN 231: * 232: * Scale matrix norm down to BIGNUM 233: * 234: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO ) 235: IASCL = 2 236: ELSE IF( ANRM.EQ.ZERO ) THEN 237: * 238: * Matrix all zero. Return zero solution. 239: * 240: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 241: RANK = 0 242: GO TO 70 243: END IF 244: * 245: BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK ) 246: IBSCL = 0 247: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN 248: * 249: * Scale matrix norm up to SMLNUM 250: * 251: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO ) 252: IBSCL = 1 253: ELSE IF( BNRM.GT.BIGNUM ) THEN 254: * 255: * Scale matrix norm down to BIGNUM 256: * 257: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO ) 258: IBSCL = 2 259: END IF 260: * 261: * Compute QR factorization with column pivoting of A: 262: * A * P = Q * R 263: * 264: CALL ZGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), 265: $ LWORK-MN, RWORK, INFO ) 266: WSIZE = MN + DBLE( WORK( MN+1 ) ) 267: * 268: * complex workspace: MN+NB*(N+1). real workspace 2*N. 269: * Details of Householder rotations stored in WORK(1:MN). 270: * 271: * Determine RANK using incremental condition estimation 272: * 273: WORK( ISMIN ) = CONE 274: WORK( ISMAX ) = CONE 275: SMAX = ABS( A( 1, 1 ) ) 276: SMIN = SMAX 277: IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN 278: RANK = 0 279: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB ) 280: GO TO 70 281: ELSE 282: RANK = 1 283: END IF 284: * 285: 10 CONTINUE 286: IF( RANK.LT.MN ) THEN 287: I = RANK + 1 288: CALL ZLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ), 289: $ A( I, I ), SMINPR, S1, C1 ) 290: CALL ZLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ), 291: $ A( I, I ), SMAXPR, S2, C2 ) 292: * 293: IF( SMAXPR*RCOND.LE.SMINPR ) THEN 294: DO 20 I = 1, RANK 295: WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 ) 296: WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 ) 297: 20 CONTINUE 298: WORK( ISMIN+RANK ) = C1 299: WORK( ISMAX+RANK ) = C2 300: SMIN = SMINPR 301: SMAX = SMAXPR 302: RANK = RANK + 1 303: GO TO 10 304: END IF 305: END IF 306: * 307: * complex workspace: 3*MN. 308: * 309: * Logically partition R = [ R11 R12 ] 310: * [ 0 R22 ] 311: * where R11 = R(1:RANK,1:RANK) 312: * 313: * [R11,R12] = [ T11, 0 ] * Y 314: * 315: IF( RANK.LT.N ) 316: $ CALL ZTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ), 317: $ LWORK-2*MN, INFO ) 318: * 319: * complex workspace: 2*MN. 320: * Details of Householder rotations stored in WORK(MN+1:2*MN) 321: * 322: * B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) 323: * 324: CALL ZUNMQR( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA, 325: $ WORK( 1 ), B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO ) 326: WSIZE = MAX( WSIZE, 2*MN+DBLE( WORK( 2*MN+1 ) ) ) 327: * 328: * complex workspace: 2*MN+NB*NRHS. 329: * 330: * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) 331: * 332: CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK, 333: $ NRHS, CONE, A, LDA, B, LDB ) 334: * 335: DO 40 J = 1, NRHS 336: DO 30 I = RANK + 1, N 337: B( I, J ) = CZERO 338: 30 CONTINUE 339: 40 CONTINUE 340: * 341: * B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) 342: * 343: IF( RANK.LT.N ) THEN 344: CALL ZUNMRZ( 'Left', 'Conjugate transpose', N, NRHS, RANK, 345: $ N-RANK, A, LDA, WORK( MN+1 ), B, LDB, 346: $ WORK( 2*MN+1 ), LWORK-2*MN, INFO ) 347: END IF 348: * 349: * complex workspace: 2*MN+NRHS. 350: * 351: * B(1:N,1:NRHS) := P * B(1:N,1:NRHS) 352: * 353: DO 60 J = 1, NRHS 354: DO 50 I = 1, N 355: WORK( JPVT( I ) ) = B( I, J ) 356: 50 CONTINUE 357: CALL ZCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 ) 358: 60 CONTINUE 359: * 360: * complex workspace: N. 361: * 362: * Undo scaling 363: * 364: IF( IASCL.EQ.1 ) THEN 365: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO ) 366: CALL ZLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA, 367: $ INFO ) 368: ELSE IF( IASCL.EQ.2 ) THEN 369: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO ) 370: CALL ZLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA, 371: $ INFO ) 372: END IF 373: IF( IBSCL.EQ.1 ) THEN 374: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO ) 375: ELSE IF( IBSCL.EQ.2 ) THEN 376: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO ) 377: END IF 378: * 379: 70 CONTINUE 380: WORK( 1 ) = DCMPLX( LWKOPT ) 381: * 382: RETURN 383: * 384: * End of ZGELSY 385: * 386: END