![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> DGGLSE solves overdetermined or underdetermined systems for OTHER matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DGGLSE + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgglse.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgglse.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgglse.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, 00022 * INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDA, LDB, LWORK, M, N, P 00026 * .. 00027 * .. Array Arguments .. 00028 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( * ), D( * ), 00029 * $ WORK( * ), X( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> DGGLSE solves the linear equality-constrained least squares (LSE) 00039 *> problem: 00040 *> 00041 *> minimize || c - A*x ||_2 subject to B*x = d 00042 *> 00043 *> where A is an M-by-N matrix, B is a P-by-N matrix, c is a given 00044 *> M-vector, and d is a given P-vector. It is assumed that 00045 *> P <= N <= M+P, and 00046 *> 00047 *> rank(B) = P and rank( (A) ) = N. 00048 *> ( (B) ) 00049 *> 00050 *> These conditions ensure that the LSE problem has a unique solution, 00051 *> which is obtained using a generalized RQ factorization of the 00052 *> matrices (B, A) given by 00053 *> 00054 *> B = (0 R)*Q, A = Z*T*Q. 00055 *> \endverbatim 00056 * 00057 * Arguments: 00058 * ========== 00059 * 00060 *> \param[in] M 00061 *> \verbatim 00062 *> M is INTEGER 00063 *> The number of rows of the matrix A. M >= 0. 00064 *> \endverbatim 00065 *> 00066 *> \param[in] N 00067 *> \verbatim 00068 *> N is INTEGER 00069 *> The number of columns of the matrices A and B. N >= 0. 00070 *> \endverbatim 00071 *> 00072 *> \param[in] P 00073 *> \verbatim 00074 *> P is INTEGER 00075 *> The number of rows of the matrix B. 0 <= P <= N <= M+P. 00076 *> \endverbatim 00077 *> 00078 *> \param[in,out] A 00079 *> \verbatim 00080 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00081 *> On entry, the M-by-N matrix A. 00082 *> On exit, the elements on and above the diagonal of the array 00083 *> contain the min(M,N)-by-N upper trapezoidal matrix T. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] LDA 00087 *> \verbatim 00088 *> LDA is INTEGER 00089 *> The leading dimension of the array A. LDA >= max(1,M). 00090 *> \endverbatim 00091 *> 00092 *> \param[in,out] B 00093 *> \verbatim 00094 *> B is DOUBLE PRECISION array, dimension (LDB,N) 00095 *> On entry, the P-by-N matrix B. 00096 *> On exit, the upper triangle of the subarray B(1:P,N-P+1:N) 00097 *> contains the P-by-P upper triangular matrix R. 00098 *> \endverbatim 00099 *> 00100 *> \param[in] LDB 00101 *> \verbatim 00102 *> LDB is INTEGER 00103 *> The leading dimension of the array B. LDB >= max(1,P). 00104 *> \endverbatim 00105 *> 00106 *> \param[in,out] C 00107 *> \verbatim 00108 *> C is DOUBLE PRECISION array, dimension (M) 00109 *> On entry, C contains the right hand side vector for the 00110 *> least squares part of the LSE problem. 00111 *> On exit, the residual sum of squares for the solution 00112 *> is given by the sum of squares of elements N-P+1 to M of 00113 *> vector C. 00114 *> \endverbatim 00115 *> 00116 *> \param[in,out] D 00117 *> \verbatim 00118 *> D is DOUBLE PRECISION array, dimension (P) 00119 *> On entry, D contains the right hand side vector for the 00120 *> constrained equation. 00121 *> On exit, D is destroyed. 00122 *> \endverbatim 00123 *> 00124 *> \param[out] X 00125 *> \verbatim 00126 *> X is DOUBLE PRECISION array, dimension (N) 00127 *> On exit, X is the solution of the LSE problem. 00128 *> \endverbatim 00129 *> 00130 *> \param[out] WORK 00131 *> \verbatim 00132 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00133 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00134 *> \endverbatim 00135 *> 00136 *> \param[in] LWORK 00137 *> \verbatim 00138 *> LWORK is INTEGER 00139 *> The dimension of the array WORK. LWORK >= max(1,M+N+P). 00140 *> For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB, 00141 *> where NB is an upper bound for the optimal blocksizes for 00142 *> DGEQRF, SGERQF, DORMQR and SORMRQ. 00143 *> 00144 *> If LWORK = -1, then a workspace query is assumed; the routine 00145 *> only calculates the optimal size of the WORK array, returns 00146 *> this value as the first entry of the WORK array, and no error 00147 *> message related to LWORK is issued by XERBLA. 00148 *> \endverbatim 00149 *> 00150 *> \param[out] INFO 00151 *> \verbatim 00152 *> INFO is INTEGER 00153 *> = 0: successful exit. 00154 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00155 *> = 1: the upper triangular factor R associated with B in the 00156 *> generalized RQ factorization of the pair (B, A) is 00157 *> singular, so that rank(B) < P; the least squares 00158 *> solution could not be computed. 00159 *> = 2: the (N-P) by (N-P) part of the upper trapezoidal factor 00160 *> T associated with A in the generalized RQ factorization 00161 *> of the pair (B, A) is singular, so that 00162 *> rank( (A) ) < N; the least squares solution could not 00163 *> ( (B) ) 00164 *> be computed. 00165 *> \endverbatim 00166 * 00167 * Authors: 00168 * ======== 00169 * 00170 *> \author Univ. of Tennessee 00171 *> \author Univ. of California Berkeley 00172 *> \author Univ. of Colorado Denver 00173 *> \author NAG Ltd. 00174 * 00175 *> \date November 2011 00176 * 00177 *> \ingroup doubleOTHERsolve 00178 * 00179 * ===================================================================== 00180 SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, 00181 $ INFO ) 00182 * 00183 * -- LAPACK driver routine (version 3.4.0) -- 00184 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00185 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00186 * November 2011 00187 * 00188 * .. Scalar Arguments .. 00189 INTEGER INFO, LDA, LDB, LWORK, M, N, P 00190 * .. 00191 * .. Array Arguments .. 00192 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( * ), D( * ), 00193 $ WORK( * ), X( * ) 00194 * .. 00195 * 00196 * ===================================================================== 00197 * 00198 * .. Parameters .. 00199 DOUBLE PRECISION ONE 00200 PARAMETER ( ONE = 1.0D+0 ) 00201 * .. 00202 * .. Local Scalars .. 00203 LOGICAL LQUERY 00204 INTEGER LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3, 00205 $ NB4, NR 00206 * .. 00207 * .. External Subroutines .. 00208 EXTERNAL DAXPY, DCOPY, DGEMV, DGGRQF, DORMQR, DORMRQ, 00209 $ DTRMV, DTRTRS, XERBLA 00210 * .. 00211 * .. External Functions .. 00212 INTEGER ILAENV 00213 EXTERNAL ILAENV 00214 * .. 00215 * .. Intrinsic Functions .. 00216 INTRINSIC INT, MAX, MIN 00217 * .. 00218 * .. Executable Statements .. 00219 * 00220 * Test the input parameters 00221 * 00222 INFO = 0 00223 MN = MIN( M, N ) 00224 LQUERY = ( LWORK.EQ.-1 ) 00225 IF( M.LT.0 ) THEN 00226 INFO = -1 00227 ELSE IF( N.LT.0 ) THEN 00228 INFO = -2 00229 ELSE IF( P.LT.0 .OR. P.GT.N .OR. P.LT.N-M ) THEN 00230 INFO = -3 00231 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00232 INFO = -5 00233 ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 00234 INFO = -7 00235 END IF 00236 * 00237 * Calculate workspace 00238 * 00239 IF( INFO.EQ.0) THEN 00240 IF( N.EQ.0 ) THEN 00241 LWKMIN = 1 00242 LWKOPT = 1 00243 ELSE 00244 NB1 = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 ) 00245 NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 ) 00246 NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, P, -1 ) 00247 NB4 = ILAENV( 1, 'DORMRQ', ' ', M, N, P, -1 ) 00248 NB = MAX( NB1, NB2, NB3, NB4 ) 00249 LWKMIN = M + N + P 00250 LWKOPT = P + MN + MAX( M, N )*NB 00251 END IF 00252 WORK( 1 ) = LWKOPT 00253 * 00254 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN 00255 INFO = -12 00256 END IF 00257 END IF 00258 * 00259 IF( INFO.NE.0 ) THEN 00260 CALL XERBLA( 'DGGLSE', -INFO ) 00261 RETURN 00262 ELSE IF( LQUERY ) THEN 00263 RETURN 00264 END IF 00265 * 00266 * Quick return if possible 00267 * 00268 IF( N.EQ.0 ) 00269 $ RETURN 00270 * 00271 * Compute the GRQ factorization of matrices B and A: 00272 * 00273 * B*Q**T = ( 0 T12 ) P Z**T*A*Q**T = ( R11 R12 ) N-P 00274 * N-P P ( 0 R22 ) M+P-N 00275 * N-P P 00276 * 00277 * where T12 and R11 are upper triangular, and Q and Z are 00278 * orthogonal. 00279 * 00280 CALL DGGRQF( P, M, N, B, LDB, WORK, A, LDA, WORK( P+1 ), 00281 $ WORK( P+MN+1 ), LWORK-P-MN, INFO ) 00282 LOPT = WORK( P+MN+1 ) 00283 * 00284 * Update c = Z**T *c = ( c1 ) N-P 00285 * ( c2 ) M+P-N 00286 * 00287 CALL DORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, WORK( P+1 ), 00288 $ C, MAX( 1, M ), WORK( P+MN+1 ), LWORK-P-MN, INFO ) 00289 LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) ) 00290 * 00291 * Solve T12*x2 = d for x2 00292 * 00293 IF( P.GT.0 ) THEN 00294 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1, 00295 $ B( 1, N-P+1 ), LDB, D, P, INFO ) 00296 * 00297 IF( INFO.GT.0 ) THEN 00298 INFO = 1 00299 RETURN 00300 END IF 00301 * 00302 * Put the solution in X 00303 * 00304 CALL DCOPY( P, D, 1, X( N-P+1 ), 1 ) 00305 * 00306 * Update c1 00307 * 00308 CALL DGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), LDA, 00309 $ D, 1, ONE, C, 1 ) 00310 END IF 00311 * 00312 * Solve R11*x1 = c1 for x1 00313 * 00314 IF( N.GT.P ) THEN 00315 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1, 00316 $ A, LDA, C, N-P, INFO ) 00317 * 00318 IF( INFO.GT.0 ) THEN 00319 INFO = 2 00320 RETURN 00321 END IF 00322 * 00323 * Put the solutions in X 00324 * 00325 CALL DCOPY( N-P, C, 1, X, 1 ) 00326 END IF 00327 * 00328 * Compute the residual vector: 00329 * 00330 IF( M.LT.N ) THEN 00331 NR = M + P - N 00332 IF( NR.GT.0 ) 00333 $ CALL DGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, M+1 ), 00334 $ LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 ) 00335 ELSE 00336 NR = P 00337 END IF 00338 IF( NR.GT.0 ) THEN 00339 CALL DTRMV( 'Upper', 'No transpose', 'Non unit', NR, 00340 $ A( N-P+1, N-P+1 ), LDA, D, 1 ) 00341 CALL DAXPY( NR, -ONE, D, 1, C( N-P+1 ), 1 ) 00342 END IF 00343 * 00344 * Backward transformation x = Q**T*x 00345 * 00346 CALL DORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), X, 00347 $ N, WORK( P+MN+1 ), LWORK-P-MN, INFO ) 00348 WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) 00349 * 00350 RETURN 00351 * 00352 * End of DGGLSE 00353 * 00354 END