![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZGGGLM( 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: COMPLEX*16 A( LDA, * ), B( LDB, * ), D( * ), WORK( * ), 14: $ X( * ), Y( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * ZGGGLM 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 array, dimension (M) 82: * Y (output) COMPLEX*16 array, dimension (P) 83: * On exit, X and Y are the solutions of the GLM problem. 84: * 85: * WORK (workspace/output) COMPLEX*16 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: * ZGEQRF, ZGERQF, ZUNMQR and ZUNMRQ. 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: COMPLEX*16 CZERO, CONE 116: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 117: $ CONE = ( 1.0D+0, 0.0D+0 ) ) 118: * .. 119: * .. Local Scalars .. 120: LOGICAL LQUERY 121: INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3, 122: $ NB4, NP 123: * .. 124: * .. External Subroutines .. 125: EXTERNAL XERBLA, ZCOPY, ZGEMV, ZGGQRF, ZTRTRS, ZUNMQR, 126: $ ZUNMRQ 127: * .. 128: * .. External Functions .. 129: INTEGER ILAENV 130: EXTERNAL ILAENV 131: * .. 132: * .. Intrinsic Functions .. 133: INTRINSIC INT, MAX, MIN 134: * .. 135: * .. Executable Statements .. 136: * 137: * Test the input parameters 138: * 139: INFO = 0 140: NP = MIN( N, P ) 141: LQUERY = ( LWORK.EQ.-1 ) 142: IF( N.LT.0 ) THEN 143: INFO = -1 144: ELSE IF( M.LT.0 .OR. M.GT.N ) THEN 145: INFO = -2 146: ELSE IF( P.LT.0 .OR. P.LT.N-M ) THEN 147: INFO = -3 148: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 149: INFO = -5 150: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 151: INFO = -7 152: END IF 153: * 154: * Calculate workspace 155: * 156: IF( INFO.EQ.0) THEN 157: IF( N.EQ.0 ) THEN 158: LWKMIN = 1 159: LWKOPT = 1 160: ELSE 161: NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, M, -1, -1 ) 162: NB2 = ILAENV( 1, 'ZGERQF', ' ', N, M, -1, -1 ) 163: NB3 = ILAENV( 1, 'ZUNMQR', ' ', N, M, P, -1 ) 164: NB4 = ILAENV( 1, 'ZUNMRQ', ' ', N, M, P, -1 ) 165: NB = MAX( NB1, NB2, NB3, NB4 ) 166: LWKMIN = M + N + P 167: LWKOPT = M + NP + MAX( N, P )*NB 168: END IF 169: WORK( 1 ) = LWKOPT 170: * 171: IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN 172: INFO = -12 173: END IF 174: END IF 175: * 176: IF( INFO.NE.0 ) THEN 177: CALL XERBLA( 'ZGGGLM', -INFO ) 178: RETURN 179: ELSE IF( LQUERY ) THEN 180: RETURN 181: END IF 182: * 183: * Quick return if possible 184: * 185: IF( N.EQ.0 ) 186: $ RETURN 187: * 188: * Compute the GQR factorization of matrices A and B: 189: * 190: * Q'*A = ( R11 ) M, Q'*B*Z' = ( T11 T12 ) M 191: * ( 0 ) N-M ( 0 T22 ) N-M 192: * M M+P-N N-M 193: * 194: * where R11 and T22 are upper triangular, and Q and Z are 195: * unitary. 196: * 197: CALL ZGGQRF( N, M, P, A, LDA, WORK, B, LDB, WORK( M+1 ), 198: $ WORK( M+NP+1 ), LWORK-M-NP, INFO ) 199: LOPT = WORK( M+NP+1 ) 200: * 201: * Update left-hand-side vector d = Q'*d = ( d1 ) M 202: * ( d2 ) N-M 203: * 204: CALL ZUNMQR( 'Left', 'Conjugate transpose', N, 1, M, A, LDA, WORK, 205: $ D, MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) 206: LOPT = MAX( LOPT, INT( WORK( M+NP+1 ) ) ) 207: * 208: * Solve T22*y2 = d2 for y2 209: * 210: IF( N.GT.M ) THEN 211: CALL ZTRTRS( 'Upper', 'No transpose', 'Non unit', N-M, 1, 212: $ B( M+1, M+P-N+1 ), LDB, D( M+1 ), N-M, INFO ) 213: * 214: IF( INFO.GT.0 ) THEN 215: INFO = 1 216: RETURN 217: END IF 218: * 219: CALL ZCOPY( N-M, D( M+1 ), 1, Y( M+P-N+1 ), 1 ) 220: END IF 221: * 222: * Set y1 = 0 223: * 224: DO 10 I = 1, M + P - N 225: Y( I ) = CZERO 226: 10 CONTINUE 227: * 228: * Update d1 = d1 - T12*y2 229: * 230: CALL ZGEMV( 'No transpose', M, N-M, -CONE, B( 1, M+P-N+1 ), LDB, 231: $ Y( M+P-N+1 ), 1, CONE, D, 1 ) 232: * 233: * Solve triangular system: R11*x = d1 234: * 235: IF( M.GT.0 ) THEN 236: CALL ZTRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, LDA, 237: $ D, M, INFO ) 238: * 239: IF( INFO.GT.0 ) THEN 240: INFO = 2 241: RETURN 242: END IF 243: * 244: * Copy D to X 245: * 246: CALL ZCOPY( M, D, 1, X, 1 ) 247: END IF 248: * 249: * Backward transformation y = Z'*y 250: * 251: CALL ZUNMRQ( 'Left', 'Conjugate transpose', P, 1, NP, 252: $ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y, 253: $ MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) 254: WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) ) 255: * 256: RETURN 257: * 258: * End of ZGGGLM 259: * 260: END