1 SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK,
2 $ INFO )
3 *
4 * -- LAPACK driver routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDA, LDB, LWORK, M, N, P
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( * ), D( * ),
14 $ WORK( * ), X( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGGLSE 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (N)
79 * On exit, X is the solution of the LSE problem.
80 *
81 * WORK (workspace/output) DOUBLE PRECISION 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 * DGEQRF, SGERQF, DORMQR and SORMRQ.
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 DOUBLE PRECISION ONE
113 PARAMETER ( ONE = 1.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 DAXPY, DCOPY, DGEMV, DGGRQF, DORMQR, DORMRQ,
122 $ DTRMV, DTRTRS, XERBLA
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, 'DGEQRF', ' ', M, N, -1, -1 )
158 NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
159 NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, P, -1 )
160 NB4 = ILAENV( 1, 'DORMRQ', ' ', 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( 'DGGLSE', -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**T = ( 0 T12 ) P Z**T*A*Q**T = ( 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 * orthogonal.
192 *
193 CALL DGGRQF( 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**T *c = ( c1 ) N-P
198 * ( c2 ) M+P-N
199 *
200 CALL DORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, WORK( P+1 ),
201 $ C, MAX( 1, M ), WORK( P+MN+1 ), LWORK-P-MN, INFO )
202 LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
203 *
204 * Solve T12*x2 = d for x2
205 *
206 IF( P.GT.0 ) THEN
207 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1,
208 $ B( 1, N-P+1 ), LDB, D, P, INFO )
209 *
210 IF( INFO.GT.0 ) THEN
211 INFO = 1
212 RETURN
213 END IF
214 *
215 * Put the solution in X
216 *
217 CALL DCOPY( P, D, 1, X( N-P+1 ), 1 )
218 *
219 * Update c1
220 *
221 CALL DGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), LDA,
222 $ D, 1, ONE, C, 1 )
223 END IF
224 *
225 * Solve R11*x1 = c1 for x1
226 *
227 IF( N.GT.P ) THEN
228 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1,
229 $ A, LDA, C, N-P, INFO )
230 *
231 IF( INFO.GT.0 ) THEN
232 INFO = 2
233 RETURN
234 END IF
235 *
236 * Put the solutions in X
237 *
238 CALL DCOPY( N-P, C, 1, X, 1 )
239 END IF
240 *
241 * Compute the residual vector:
242 *
243 IF( M.LT.N ) THEN
244 NR = M + P - N
245 IF( NR.GT.0 )
246 $ CALL DGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, M+1 ),
247 $ LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 )
248 ELSE
249 NR = P
250 END IF
251 IF( NR.GT.0 ) THEN
252 CALL DTRMV( 'Upper', 'No transpose', 'Non unit', NR,
253 $ A( N-P+1, N-P+1 ), LDA, D, 1 )
254 CALL DAXPY( NR, -ONE, D, 1, C( N-P+1 ), 1 )
255 END IF
256 *
257 * Backward transformation x = Q**T*x
258 *
259 CALL DORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), X,
260 $ N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
261 WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
262 *
263 RETURN
264 *
265 * End of DGGLSE
266 *
267 END
2 $ INFO )
3 *
4 * -- LAPACK driver routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 INTEGER INFO, LDA, LDB, LWORK, M, N, P
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( * ), D( * ),
14 $ WORK( * ), X( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGGLSE 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (N)
79 * On exit, X is the solution of the LSE problem.
80 *
81 * WORK (workspace/output) DOUBLE PRECISION 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 * DGEQRF, SGERQF, DORMQR and SORMRQ.
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 DOUBLE PRECISION ONE
113 PARAMETER ( ONE = 1.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 DAXPY, DCOPY, DGEMV, DGGRQF, DORMQR, DORMRQ,
122 $ DTRMV, DTRTRS, XERBLA
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, 'DGEQRF', ' ', M, N, -1, -1 )
158 NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
159 NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, P, -1 )
160 NB4 = ILAENV( 1, 'DORMRQ', ' ', 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( 'DGGLSE', -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**T = ( 0 T12 ) P Z**T*A*Q**T = ( 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 * orthogonal.
192 *
193 CALL DGGRQF( 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**T *c = ( c1 ) N-P
198 * ( c2 ) M+P-N
199 *
200 CALL DORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, WORK( P+1 ),
201 $ C, MAX( 1, M ), WORK( P+MN+1 ), LWORK-P-MN, INFO )
202 LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
203 *
204 * Solve T12*x2 = d for x2
205 *
206 IF( P.GT.0 ) THEN
207 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1,
208 $ B( 1, N-P+1 ), LDB, D, P, INFO )
209 *
210 IF( INFO.GT.0 ) THEN
211 INFO = 1
212 RETURN
213 END IF
214 *
215 * Put the solution in X
216 *
217 CALL DCOPY( P, D, 1, X( N-P+1 ), 1 )
218 *
219 * Update c1
220 *
221 CALL DGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), LDA,
222 $ D, 1, ONE, C, 1 )
223 END IF
224 *
225 * Solve R11*x1 = c1 for x1
226 *
227 IF( N.GT.P ) THEN
228 CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1,
229 $ A, LDA, C, N-P, INFO )
230 *
231 IF( INFO.GT.0 ) THEN
232 INFO = 2
233 RETURN
234 END IF
235 *
236 * Put the solutions in X
237 *
238 CALL DCOPY( N-P, C, 1, X, 1 )
239 END IF
240 *
241 * Compute the residual vector:
242 *
243 IF( M.LT.N ) THEN
244 NR = M + P - N
245 IF( NR.GT.0 )
246 $ CALL DGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, M+1 ),
247 $ LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 )
248 ELSE
249 NR = P
250 END IF
251 IF( NR.GT.0 ) THEN
252 CALL DTRMV( 'Upper', 'No transpose', 'Non unit', NR,
253 $ A( N-P+1, N-P+1 ), LDA, D, 1 )
254 CALL DAXPY( NR, -ONE, D, 1, C( N-P+1 ), 1 )
255 END IF
256 *
257 * Backward transformation x = Q**T*x
258 *
259 CALL DORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), X,
260 $ N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
261 WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
262 *
263 RETURN
264 *
265 * End of DGGLSE
266 *
267 END