![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, 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: INTEGER INFO, LDA, LDB, LWORK, M, N, P 11: * .. 12: * .. Array Arguments .. 13: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), D( * ), WORK( * ), 14: $ X( * ), Y( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * DGGGLM solves a general Gauss-Markov linear model (GLM) problem: 21: * 22: * minimize || y ||_2 subject to d = A*x + B*y 23: * x 24: * 25: * where A is an N-by-M matrix, B is an N-by-P matrix, and d is a 26: * given N-vector. It is assumed that M <= N <= M+P, and 27: * 28: * rank(A) = M and rank( A B ) = N. 29: * 30: * Under these assumptions, the constrained equation is always 31: * consistent, and there is a unique solution x and a minimal 2-norm 32: * solution y, which is obtained using a generalized QR factorization 33: * of the matrices (A, B) given by 34: * 35: * A = Q*(R), B = Q*T*Z. 36: * (0) 37: * 38: * In particular, if matrix B is square nonsingular, then the problem 39: * GLM is equivalent to the following weighted linear least squares 40: * problem 41: * 42: * minimize || inv(B)*(d-A*x) ||_2 43: * x 44: * 45: * where inv(B) denotes the inverse of B. 46: * 47: * Arguments 48: * ========= 49: * 50: * N (input) INTEGER 51: * The number of rows of the matrices A and B. N >= 0. 52: * 53: * M (input) INTEGER 54: * The number of columns of the matrix A. 0 <= M <= N. 55: * 56: * P (input) INTEGER 57: * The number of columns of the matrix B. P >= N-M. 58: * 59: * A (input/output) DOUBLE PRECISION array, dimension (LDA,M) 60: * On entry, the N-by-M matrix A. 61: * On exit, the upper triangular part of the array A contains 62: * the M-by-M upper triangular matrix R. 63: * 64: * LDA (input) INTEGER 65: * The leading dimension of the array A. LDA >= max(1,N). 66: * 67: * B (input/output) DOUBLE PRECISION array, dimension (LDB,P) 68: * On entry, the N-by-P matrix B. 69: * On exit, if N <= P, the upper triangle of the subarray 70: * B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T; 71: * if N > P, the elements on and above the (N-P)th subdiagonal 72: * contain the N-by-P upper trapezoidal matrix T. 73: * 74: * LDB (input) INTEGER 75: * The leading dimension of the array B. LDB >= max(1,N). 76: * 77: * D (input/output) DOUBLE PRECISION array, dimension (N) 78: * On entry, D is the left hand side of the GLM equation. 79: * On exit, D is destroyed. 80: * 81: * X (output) DOUBLE PRECISION array, dimension (M) 82: * Y (output) DOUBLE PRECISION array, dimension (P) 83: * On exit, X and Y are the solutions of the GLM problem. 84: * 85: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 86: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 87: * 88: * LWORK (input) INTEGER 89: * The dimension of the array WORK. LWORK >= max(1,N+M+P). 90: * For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB, 91: * where NB is an upper bound for the optimal blocksizes for 92: * DGEQRF, SGERQF, DORMQR and SORMRQ. 93: * 94: * If LWORK = -1, then a workspace query is assumed; the routine 95: * only calculates the optimal size of the WORK array, returns 96: * this value as the first entry of the WORK array, and no error 97: * message related to LWORK is issued by XERBLA. 98: * 99: * INFO (output) INTEGER 100: * = 0: successful exit. 101: * < 0: if INFO = -i, the i-th argument had an illegal value. 102: * = 1: the upper triangular factor R associated with A in the 103: * generalized QR factorization of the pair (A, B) is 104: * singular, so that rank(A) < M; the least squares 105: * solution could not be computed. 106: * = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal 107: * factor T associated with B in the generalized QR 108: * factorization of the pair (A, B) is singular, so that 109: * rank( A B ) < N; the least squares solution could not 110: * be computed. 111: * 112: * =================================================================== 113: * 114: * .. Parameters .. 115: DOUBLE PRECISION ZERO, ONE 116: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 117: * .. 118: * .. Local Scalars .. 119: LOGICAL LQUERY 120: INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3, 121: $ NB4, NP 122: * .. 123: * .. External Subroutines .. 124: EXTERNAL DCOPY, DGEMV, DGGQRF, DORMQR, DORMRQ, DTRTRS, 125: $ XERBLA 126: * .. 127: * .. External Functions .. 128: INTEGER ILAENV 129: EXTERNAL ILAENV 130: * .. 131: * .. Intrinsic Functions .. 132: INTRINSIC INT, MAX, MIN 133: * .. 134: * .. Executable Statements .. 135: * 136: * Test the input parameters 137: * 138: INFO = 0 139: NP = MIN( N, P ) 140: LQUERY = ( LWORK.EQ.-1 ) 141: IF( N.LT.0 ) THEN 142: INFO = -1 143: ELSE IF( M.LT.0 .OR. M.GT.N ) THEN 144: INFO = -2 145: ELSE IF( P.LT.0 .OR. P.LT.N-M ) THEN 146: INFO = -3 147: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 148: INFO = -5 149: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 150: INFO = -7 151: END IF 152: * 153: * Calculate workspace 154: * 155: IF( INFO.EQ.0) THEN 156: IF( N.EQ.0 ) THEN 157: LWKMIN = 1 158: LWKOPT = 1 159: ELSE 160: NB1 = ILAENV( 1, 'DGEQRF', ' ', N, M, -1, -1 ) 161: NB2 = ILAENV( 1, 'DGERQF', ' ', N, M, -1, -1 ) 162: NB3 = ILAENV( 1, 'DORMQR', ' ', N, M, P, -1 ) 163: NB4 = ILAENV( 1, 'DORMRQ', ' ', N, M, P, -1 ) 164: NB = MAX( NB1, NB2, NB3, NB4 ) 165: LWKMIN = M + N + P 166: LWKOPT = M + NP + MAX( N, P )*NB 167: END IF 168: WORK( 1 ) = LWKOPT 169: * 170: IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN 171: INFO = -12 172: END IF 173: END IF 174: * 175: IF( INFO.NE.0 ) THEN 176: CALL XERBLA( 'DGGGLM', -INFO ) 177: RETURN 178: ELSE IF( LQUERY ) THEN 179: RETURN 180: END IF 181: * 182: * Quick return if possible 183: * 184: IF( N.EQ.0 ) 185: $ RETURN 186: * 187: * Compute the GQR factorization of matrices A and B: 188: * 189: * Q'*A = ( R11 ) M, Q'*B*Z' = ( T11 T12 ) M 190: * ( 0 ) N-M ( 0 T22 ) N-M 191: * M M+P-N N-M 192: * 193: * where R11 and T22 are upper triangular, and Q and Z are 194: * orthogonal. 195: * 196: CALL DGGQRF( N, M, P, A, LDA, WORK, B, LDB, WORK( M+1 ), 197: $ WORK( M+NP+1 ), LWORK-M-NP, INFO ) 198: LOPT = WORK( M+NP+1 ) 199: * 200: * Update left-hand-side vector d = Q'*d = ( d1 ) M 201: * ( d2 ) N-M 202: * 203: CALL DORMQR( 'Left', 'Transpose', N, 1, M, A, LDA, WORK, D, 204: $ MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) 205: LOPT = MAX( LOPT, INT( WORK( M+NP+1 ) ) ) 206: * 207: * Solve T22*y2 = d2 for y2 208: * 209: IF( N.GT.M ) THEN 210: CALL DTRTRS( 'Upper', 'No transpose', 'Non unit', N-M, 1, 211: $ B( M+1, M+P-N+1 ), LDB, D( M+1 ), N-M, INFO ) 212: * 213: IF( INFO.GT.0 ) THEN 214: INFO = 1 215: RETURN 216: END IF 217: * 218: CALL DCOPY( N-M, D( M+1 ), 1, Y( M+P-N+1 ), 1 ) 219: END IF 220: * 221: * Set y1 = 0 222: * 223: DO 10 I = 1, M + P - N 224: Y( I ) = ZERO 225: 10 CONTINUE 226: * 227: * Update d1 = d1 - T12*y2 228: * 229: CALL DGEMV( 'No transpose', M, N-M, -ONE, B( 1, M+P-N+1 ), LDB, 230: $ Y( M+P-N+1 ), 1, ONE, D, 1 ) 231: * 232: * Solve triangular system: R11*x = d1 233: * 234: IF( M.GT.0 ) THEN 235: CALL DTRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, LDA, 236: $ D, M, INFO ) 237: * 238: IF( INFO.GT.0 ) THEN 239: INFO = 2 240: RETURN 241: END IF 242: * 243: * Copy D to X 244: * 245: CALL DCOPY( M, D, 1, X, 1 ) 246: END IF 247: * 248: * Backward transformation y = Z'*y 249: * 250: CALL DORMRQ( 'Left', 'Transpose', P, 1, NP, 251: $ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y, 252: $ MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) 253: WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) ) 254: * 255: RETURN 256: * 257: * End of DGGGLM 258: * 259: END