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