1 SUBROUTINE ZGGLSE( 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 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**H = ( 0 T12 ) P Z**H*A*Q**H = ( 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**H *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**H*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
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 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**H = ( 0 T12 ) P Z**H*A*Q**H = ( 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**H *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**H*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