 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, 00002 $ INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDA, LDB, LWORK, M, N, P 00011 * .. 00012 * .. Array Arguments .. 00013 COMPLEX A( LDA, * ), B( LDB, * ), C( * ), D( * ), 00014 $ WORK( * ), X( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * CGGLSE solves the linear equality-constrained least squares (LSE) 00021 * problem: 00022 * 00023 * minimize || c - A*x ||_2 subject to B*x = d 00024 * 00025 * where A is an M-by-N matrix, B is a P-by-N matrix, c is a given 00026 * M-vector, and d is a given P-vector. It is assumed that 00027 * P <= N <= M+P, and 00028 * 00029 * rank(B) = P and rank( (A) ) = N. 00030 * ( (B) ) 00031 * 00032 * These conditions ensure that the LSE problem has a unique solution, 00033 * which is obtained using a generalized RQ factorization of the 00034 * matrices (B, A) given by 00035 * 00036 * B = (0 R)*Q, A = Z*T*Q. 00037 * 00038 * Arguments 00039 * ========= 00040 * 00041 * M (input) INTEGER 00042 * The number of rows of the matrix A. M >= 0. 00043 * 00044 * N (input) INTEGER 00045 * The number of columns of the matrices A and B. N >= 0. 00046 * 00047 * P (input) INTEGER 00048 * The number of rows of the matrix B. 0 <= P <= N <= M+P. 00049 * 00050 * A (input/output) COMPLEX array, dimension (LDA,N) 00051 * On entry, the M-by-N matrix A. 00052 * On exit, the elements on and above the diagonal of the array 00053 * contain the min(M,N)-by-N upper trapezoidal matrix T. 00054 * 00055 * LDA (input) INTEGER 00056 * The leading dimension of the array A. LDA >= max(1,M). 00057 * 00058 * B (input/output) COMPLEX array, dimension (LDB,N) 00059 * On entry, the P-by-N matrix B. 00060 * On exit, the upper triangle of the subarray B(1:P,N-P+1:N) 00061 * contains the P-by-P upper triangular matrix R. 00062 * 00063 * LDB (input) INTEGER 00064 * The leading dimension of the array B. LDB >= max(1,P). 00065 * 00066 * C (input/output) COMPLEX array, dimension (M) 00067 * On entry, C contains the right hand side vector for the 00068 * least squares part of the LSE problem. 00069 * On exit, the residual sum of squares for the solution 00070 * is given by the sum of squares of elements N-P+1 to M of 00071 * vector C. 00072 * 00073 * D (input/output) COMPLEX array, dimension (P) 00074 * On entry, D contains the right hand side vector for the 00075 * constrained equation. 00076 * On exit, D is destroyed. 00077 * 00078 * X (output) COMPLEX array, dimension (N) 00079 * On exit, X is the solution of the LSE problem. 00080 * 00081 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00082 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00083 * 00084 * LWORK (input) INTEGER 00085 * The dimension of the array WORK. LWORK >= max(1,M+N+P). 00086 * For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB, 00087 * where NB is an upper bound for the optimal blocksizes for 00088 * CGEQRF, CGERQF, CUNMQR and CUNMRQ. 00089 * 00090 * If LWORK = -1, then a workspace query is assumed; the routine 00091 * only calculates the optimal size of the WORK array, returns 00092 * this value as the first entry of the WORK array, and no error 00093 * message related to LWORK is issued by XERBLA. 00094 * 00095 * INFO (output) INTEGER 00096 * = 0: successful exit. 00097 * < 0: if INFO = -i, the i-th argument had an illegal value. 00098 * = 1: the upper triangular factor R associated with B in the 00099 * generalized RQ factorization of the pair (B, A) is 00100 * singular, so that rank(B) < P; the least squares 00101 * solution could not be computed. 00102 * = 2: the (N-P) by (N-P) part of the upper trapezoidal factor 00103 * T associated with A in the generalized RQ factorization 00104 * of the pair (B, A) is singular, so that 00105 * rank( (A) ) < N; the least squares solution could not 00106 * ( (B) ) 00107 * be computed. 00108 * 00109 * ===================================================================== 00110 * 00111 * .. Parameters .. 00112 COMPLEX CONE 00113 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) 00114 * .. 00115 * .. Local Scalars .. 00116 LOGICAL LQUERY 00117 INTEGER LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3, 00118 $ NB4, NR 00119 * .. 00120 * .. External Subroutines .. 00121 EXTERNAL CAXPY, CCOPY, CGEMV, CGGRQF, CTRMV, CTRTRS, 00122 $ CUNMQR, CUNMRQ, XERBLA 00123 * .. 00124 * .. External Functions .. 00125 INTEGER ILAENV 00126 EXTERNAL ILAENV 00127 * .. 00128 * .. Intrinsic Functions .. 00129 INTRINSIC INT, MAX, MIN 00130 * .. 00131 * .. Executable Statements .. 00132 * 00133 * Test the input parameters 00134 * 00135 INFO = 0 00136 MN = MIN( M, N ) 00137 LQUERY = ( LWORK.EQ.-1 ) 00138 IF( M.LT.0 ) THEN 00139 INFO = -1 00140 ELSE IF( N.LT.0 ) THEN 00141 INFO = -2 00142 ELSE IF( P.LT.0 .OR. P.GT.N .OR. P.LT.N-M ) THEN 00143 INFO = -3 00144 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00145 INFO = -5 00146 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00147 INFO = -7 00148 END IF 00149 * 00150 * Calculate workspace 00151 * 00152 IF( INFO.EQ.0) THEN 00153 IF( N.EQ.0 ) THEN 00154 LWKMIN = 1 00155 LWKOPT = 1 00156 ELSE 00157 NB1 = ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 ) 00158 NB2 = ILAENV( 1, 'CGERQF', ' ', M, N, -1, -1 ) 00159 NB3 = ILAENV( 1, 'CUNMQR', ' ', M, N, P, -1 ) 00160 NB4 = ILAENV( 1, 'CUNMRQ', ' ', M, N, P, -1 ) 00161 NB = MAX( NB1, NB2, NB3, NB4 ) 00162 LWKMIN = M + N + P 00163 LWKOPT = P + MN + MAX( M, N )*NB 00164 END IF 00165 WORK( 1 ) = LWKOPT 00166 * 00167 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN 00168 INFO = -12 00169 END IF 00170 END IF 00171 * 00172 IF( INFO.NE.0 ) THEN 00173 CALL XERBLA( 'CGGLSE', -INFO ) 00174 RETURN 00175 ELSE IF( LQUERY ) THEN 00176 RETURN 00177 END IF 00178 * 00179 * Quick return if possible 00180 * 00181 IF( N.EQ.0 ) 00182 $ RETURN 00183 * 00184 * Compute the GRQ factorization of matrices B and A: 00185 * 00186 * B*Q**H = ( 0 T12 ) P Z**H*A*Q**H = ( R11 R12 ) N-P 00187 * N-P P ( 0 R22 ) M+P-N 00188 * N-P P 00189 * 00190 * where T12 and R11 are upper triangular, and Q and Z are 00191 * unitary. 00192 * 00193 CALL CGGRQF( P, M, N, B, LDB, WORK, A, LDA, WORK( P+1 ), 00194 $ WORK( P+MN+1 ), LWORK-P-MN, INFO ) 00195 LOPT = WORK( P+MN+1 ) 00196 * 00197 * Update c = Z**H *c = ( c1 ) N-P 00198 * ( c2 ) M+P-N 00199 * 00200 CALL CUNMQR( 'Left', 'Conjugate Transpose', M, 1, MN, A, LDA, 00201 $ WORK( P+1 ), C, MAX( 1, M ), WORK( P+MN+1 ), 00202 $ LWORK-P-MN, INFO ) 00203 LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) ) 00204 * 00205 * Solve T12*x2 = d for x2 00206 * 00207 IF( P.GT.0 ) THEN 00208 CALL CTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1, 00209 $ B( 1, N-P+1 ), LDB, D, P, INFO ) 00210 * 00211 IF( INFO.GT.0 ) THEN 00212 INFO = 1 00213 RETURN 00214 END IF 00215 * 00216 * Put the solution in X 00217 * 00218 CALL CCOPY( P, D, 1, X( N-P+1 ), 1 ) 00219 * 00220 * Update c1 00221 * 00222 CALL CGEMV( 'No transpose', N-P, P, -CONE, A( 1, N-P+1 ), LDA, 00223 $ D, 1, CONE, C, 1 ) 00224 END IF 00225 * 00226 * Solve R11*x1 = c1 for x1 00227 * 00228 IF( N.GT.P ) THEN 00229 CALL CTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1, 00230 $ A, LDA, C, N-P, INFO ) 00231 * 00232 IF( INFO.GT.0 ) THEN 00233 INFO = 2 00234 RETURN 00235 END IF 00236 * 00237 * Put the solutions in X 00238 * 00239 CALL CCOPY( N-P, C, 1, X, 1 ) 00240 END IF 00241 * 00242 * Compute the residual vector: 00243 * 00244 IF( M.LT.N ) THEN 00245 NR = M + P - N 00246 IF( NR.GT.0 ) 00247 $ CALL CGEMV( 'No transpose', NR, N-M, -CONE, A( N-P+1, M+1 ), 00248 $ LDA, D( NR+1 ), 1, CONE, C( N-P+1 ), 1 ) 00249 ELSE 00250 NR = P 00251 END IF 00252 IF( NR.GT.0 ) THEN 00253 CALL CTRMV( 'Upper', 'No transpose', 'Non unit', NR, 00254 $ A( N-P+1, N-P+1 ), LDA, D, 1 ) 00255 CALL CAXPY( NR, -CONE, D, 1, C( N-P+1 ), 1 ) 00256 END IF 00257 * 00258 * Backward transformation x = Q**H*x 00259 * 00260 CALL CUNMRQ( 'Left', 'Conjugate Transpose', N, 1, P, B, LDB, 00261 $ WORK( 1 ), X, N, WORK( P+MN+1 ), LWORK-P-MN, INFO ) 00262 WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) 00263 * 00264 RETURN 00265 * 00266 * End of CGGLSE 00267 * 00268 END
1.7.3 
