|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, 00002 $ INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDA, LDB, LWORK, M, N, P 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( * ), D( * ), 00014 $ WORK( * ), X( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DGGLSE 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (N) 00079 * On exit, X is the solution of the LSE problem. 00080 * 00081 * WORK (workspace/output) DOUBLE PRECISION 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 * DGEQRF, SGERQF, DORMQR and SORMRQ. 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 DOUBLE PRECISION ONE 00113 PARAMETER ( ONE = 1.0D+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 DAXPY, DCOPY, DGEMV, DGGRQF, DORMQR, DORMRQ, 00122 $ DTRMV, DTRTRS, 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, 'DGEQRF', ' ', M, N, -1, -1 ) 00158 NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 ) 00159 NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, P, -1 ) 00160 NB4 = ILAENV( 1, 'DORMRQ', ' ', 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( 'DGGLSE', -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' = ( 0 T12 ) P Z'*A*Q' = ( 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 * orthogonal. 00192 * 00193 CALL DGGRQF( 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'*c = ( c1 ) N-P 00198 * ( c2 ) M+P-N 00199 * 00200 CALL DORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, WORK( P+1 ), 00201 $ C, MAX( 1, M ), WORK( P+MN+1 ), LWORK-P-MN, INFO ) 00202 LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) ) 00203 * 00204 * Solve T12*x2 = d for x2 00205 * 00206 IF( P.GT.0 ) THEN 00207 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1, 00208 $ B( 1, N-P+1 ), LDB, D, P, INFO ) 00209 * 00210 IF( INFO.GT.0 ) THEN 00211 INFO = 1 00212 RETURN 00213 END IF 00214 * 00215 * Put the solution in X 00216 * 00217 CALL DCOPY( P, D, 1, X( N-P+1 ), 1 ) 00218 * 00219 * Update c1 00220 * 00221 CALL DGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), LDA, 00222 $ D, 1, ONE, C, 1 ) 00223 END IF 00224 * 00225 * Solve R11*x1 = c1 for x1 00226 * 00227 IF( N.GT.P ) THEN 00228 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1, 00229 $ A, LDA, C, N-P, INFO ) 00230 * 00231 IF( INFO.GT.0 ) THEN 00232 INFO = 2 00233 RETURN 00234 END IF 00235 * 00236 * Put the solutions in X 00237 * 00238 CALL DCOPY( N-P, C, 1, X, 1 ) 00239 END IF 00240 * 00241 * Compute the residual vector: 00242 * 00243 IF( M.LT.N ) THEN 00244 NR = M + P - N 00245 IF( NR.GT.0 ) 00246 $ CALL DGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, M+1 ), 00247 $ LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 ) 00248 ELSE 00249 NR = P 00250 END IF 00251 IF( NR.GT.0 ) THEN 00252 CALL DTRMV( 'Upper', 'No transpose', 'Non unit', NR, 00253 $ A( N-P+1, N-P+1 ), LDA, D, 1 ) 00254 CALL DAXPY( NR, -ONE, D, 1, C( N-P+1 ), 1 ) 00255 END IF 00256 * 00257 * Backward transformation x = Q'*x 00258 * 00259 CALL DORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), X, 00260 $ N, WORK( P+MN+1 ), LWORK-P-MN, INFO ) 00261 WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) 00262 * 00263 RETURN 00264 * 00265 * End of DGGLSE 00266 * 00267 END
1.7.3 
