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