LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cggglm.f
Go to the documentation of this file.
1*> \brief \b CGGGLM
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGGGLM + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggglm.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggglm.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggglm.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK,
22* INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INFO, LDA, LDB, LWORK, M, N, P
26* ..
27* .. Array Arguments ..
28* COMPLEX A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
29* $ X( * ), Y( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> CGGGLM solves a general Gauss-Markov linear model (GLM) problem:
39*>
40*> minimize || y ||_2 subject to d = A*x + B*y
41*> x
42*>
43*> where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
44*> given N-vector. It is assumed that M <= N <= M+P, and
45*>
46*> rank(A) = M and rank( A B ) = N.
47*>
48*> Under these assumptions, the constrained equation is always
49*> consistent, and there is a unique solution x and a minimal 2-norm
50*> solution y, which is obtained using a generalized QR factorization
51*> of the matrices (A, B) given by
52*>
53*> A = Q*(R), B = Q*T*Z.
54*> (0)
55*>
56*> In particular, if matrix B is square nonsingular, then the problem
57*> GLM is equivalent to the following weighted linear least squares
58*> problem
59*>
60*> minimize || inv(B)*(d-A*x) ||_2
61*> x
62*>
63*> where inv(B) denotes the inverse of B.
64*> \endverbatim
65*
66* Arguments:
67* ==========
68*
69*> \param[in] N
70*> \verbatim
71*> N is INTEGER
72*> The number of rows of the matrices A and B. N >= 0.
73*> \endverbatim
74*>
75*> \param[in] M
76*> \verbatim
77*> M is INTEGER
78*> The number of columns of the matrix A. 0 <= M <= N.
79*> \endverbatim
80*>
81*> \param[in] P
82*> \verbatim
83*> P is INTEGER
84*> The number of columns of the matrix B. P >= N-M.
85*> \endverbatim
86*>
87*> \param[in,out] A
88*> \verbatim
89*> A is COMPLEX array, dimension (LDA,M)
90*> On entry, the N-by-M matrix A.
91*> On exit, the upper triangular part of the array A contains
92*> the M-by-M upper triangular matrix R.
93*> \endverbatim
94*>
95*> \param[in] LDA
96*> \verbatim
97*> LDA is INTEGER
98*> The leading dimension of the array A. LDA >= max(1,N).
99*> \endverbatim
100*>
101*> \param[in,out] B
102*> \verbatim
103*> B is COMPLEX array, dimension (LDB,P)
104*> On entry, the N-by-P matrix B.
105*> On exit, if N <= P, the upper triangle of the subarray
106*> B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
107*> if N > P, the elements on and above the (N-P)th subdiagonal
108*> contain the N-by-P upper trapezoidal matrix T.
109*> \endverbatim
110*>
111*> \param[in] LDB
112*> \verbatim
113*> LDB is INTEGER
114*> The leading dimension of the array B. LDB >= max(1,N).
115*> \endverbatim
116*>
117*> \param[in,out] D
118*> \verbatim
119*> D is COMPLEX array, dimension (N)
120*> On entry, D is the left hand side of the GLM equation.
121*> On exit, D is destroyed.
122*> \endverbatim
123*>
124*> \param[out] X
125*> \verbatim
126*> X is COMPLEX array, dimension (M)
127*> \endverbatim
128*>
129*> \param[out] Y
130*> \verbatim
131*> Y is COMPLEX array, dimension (P)
132*>
133*> On exit, X and Y are the solutions of the GLM problem.
134*> \endverbatim
135*>
136*> \param[out] WORK
137*> \verbatim
138*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
139*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
140*> \endverbatim
141*>
142*> \param[in] LWORK
143*> \verbatim
144*> LWORK is INTEGER
145*> The dimension of the array WORK. LWORK >= max(1,N+M+P).
146*> For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB,
147*> where NB is an upper bound for the optimal blocksizes for
148*> CGEQRF, CGERQF, CUNMQR and CUNMRQ.
149*>
150*> If LWORK = -1, then a workspace query is assumed; the routine
151*> only calculates the optimal size of the WORK array, returns
152*> this value as the first entry of the WORK array, and no error
153*> message related to LWORK is issued by XERBLA.
154*> \endverbatim
155*>
156*> \param[out] INFO
157*> \verbatim
158*> INFO is INTEGER
159*> = 0: successful exit.
160*> < 0: if INFO = -i, the i-th argument had an illegal value.
161*> = 1: the upper triangular factor R associated with A in the
162*> generalized QR factorization of the pair (A, B) is
163*> singular, so that rank(A) < M; the least squares
164*> solution could not be computed.
165*> = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal
166*> factor T associated with B in the generalized QR
167*> factorization of the pair (A, B) is singular, so that
168*> rank( A B ) < N; the least squares solution could not
169*> be computed.
170*> \endverbatim
171*
172* Authors:
173* ========
174*
175*> \author Univ. of Tennessee
176*> \author Univ. of California Berkeley
177*> \author Univ. of Colorado Denver
178*> \author NAG Ltd.
179*
180*> \ingroup ggglm
181*
182* =====================================================================
183 SUBROUTINE cggglm( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK,
184 $ INFO )
185*
186* -- LAPACK driver routine --
187* -- LAPACK is a software package provided by Univ. of Tennessee, --
188* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189*
190* .. Scalar Arguments ..
191 INTEGER INFO, LDA, LDB, LWORK, M, N, P
192* ..
193* .. Array Arguments ..
194 COMPLEX A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
195 $ x( * ), y( * )
196* ..
197*
198* ===================================================================
199*
200* .. Parameters ..
201 COMPLEX CZERO, CONE
202 parameter( czero = ( 0.0e+0, 0.0e+0 ),
203 $ cone = ( 1.0e+0, 0.0e+0 ) )
204* ..
205* .. Local Scalars ..
206 LOGICAL LQUERY
207 INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
208 $ nb4, np
209* ..
210* .. External Subroutines ..
211 EXTERNAL ccopy, cgemv, cggqrf, ctrtrs, cunmqr, cunmrq,
212 $ xerbla
213* ..
214* .. External Functions ..
215 INTEGER ILAENV
216 REAL SROUNDUP_LWORK
217 EXTERNAL ilaenv, sroundup_lwork
218* ..
219* .. Intrinsic Functions ..
220 INTRINSIC int, max, min
221* ..
222* .. Executable Statements ..
223*
224* Test the input parameters
225*
226 info = 0
227 np = min( n, p )
228 lquery = ( lwork.EQ.-1 )
229 IF( n.LT.0 ) THEN
230 info = -1
231 ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
232 info = -2
233 ELSE IF( p.LT.0 .OR. p.LT.n-m ) THEN
234 info = -3
235 ELSE IF( lda.LT.max( 1, n ) ) THEN
236 info = -5
237 ELSE IF( ldb.LT.max( 1, n ) ) THEN
238 info = -7
239 END IF
240*
241* Calculate workspace
242*
243 IF( info.EQ.0) THEN
244 IF( n.EQ.0 ) THEN
245 lwkmin = 1
246 lwkopt = 1
247 ELSE
248 nb1 = ilaenv( 1, 'CGEQRF', ' ', n, m, -1, -1 )
249 nb2 = ilaenv( 1, 'CGERQF', ' ', n, m, -1, -1 )
250 nb3 = ilaenv( 1, 'CUNMQR', ' ', n, m, p, -1 )
251 nb4 = ilaenv( 1, 'CUNMRQ', ' ', n, m, p, -1 )
252 nb = max( nb1, nb2, nb3, nb4 )
253 lwkmin = m + n + p
254 lwkopt = m + np + max( n, p )*nb
255 END IF
256 work( 1 ) = sroundup_lwork(lwkopt)
257*
258 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
259 info = -12
260 END IF
261 END IF
262*
263 IF( info.NE.0 ) THEN
264 CALL xerbla( 'CGGGLM', -info )
265 RETURN
266 ELSE IF( lquery ) THEN
267 RETURN
268 END IF
269*
270* Quick return if possible
271*
272 IF( n.EQ.0 ) THEN
273 DO i = 1, m
274 x(i) = czero
275 END DO
276 DO i = 1, p
277 y(i) = czero
278 END DO
279 RETURN
280 END IF
281*
282* Compute the GQR factorization of matrices A and B:
283*
284* Q**H*A = ( R11 ) M, Q**H*B*Z**H = ( T11 T12 ) M
285* ( 0 ) N-M ( 0 T22 ) N-M
286* M M+P-N N-M
287*
288* where R11 and T22 are upper triangular, and Q and Z are
289* unitary.
290*
291 CALL cggqrf( n, m, p, a, lda, work, b, ldb, work( m+1 ),
292 $ work( m+np+1 ), lwork-m-np, info )
293 lopt = int( work( m+np+1 ) )
294*
295* Update left-hand-side vector d = Q**H*d = ( d1 ) M
296* ( d2 ) N-M
297*
298 CALL cunmqr( 'Left', 'Conjugate transpose', n, 1, m, a, lda, work,
299 $ d, max( 1, n ), work( m+np+1 ), lwork-m-np, info )
300 lopt = max( lopt, int( work( m+np+1 ) ) )
301*
302* Solve T22*y2 = d2 for y2
303*
304 IF( n.GT.m ) THEN
305 CALL ctrtrs( 'Upper', 'No transpose', 'Non unit', n-m, 1,
306 $ b( m+1, m+p-n+1 ), ldb, d( m+1 ), n-m, info )
307*
308 IF( info.GT.0 ) THEN
309 info = 1
310 RETURN
311 END IF
312*
313 CALL ccopy( n-m, d( m+1 ), 1, y( m+p-n+1 ), 1 )
314 END IF
315*
316* Set y1 = 0
317*
318 DO 10 i = 1, m + p - n
319 y( i ) = czero
320 10 CONTINUE
321*
322* Update d1 = d1 - T12*y2
323*
324 CALL cgemv( 'No transpose', m, n-m, -cone, b( 1, m+p-n+1 ), ldb,
325 $ y( m+p-n+1 ), 1, cone, d, 1 )
326*
327* Solve triangular system: R11*x = d1
328*
329 IF( m.GT.0 ) THEN
330 CALL ctrtrs( 'Upper', 'No Transpose', 'Non unit', m, 1, a, lda,
331 $ d, m, info )
332*
333 IF( info.GT.0 ) THEN
334 info = 2
335 RETURN
336 END IF
337*
338* Copy D to X
339*
340 CALL ccopy( m, d, 1, x, 1 )
341 END IF
342*
343* Backward transformation y = Z**H *y
344*
345 CALL cunmrq( 'Left', 'Conjugate transpose', p, 1, np,
346 $ b( max( 1, n-p+1 ), 1 ), ldb, work( m+1 ), y,
347 $ max( 1, p ), work( m+np+1 ), lwork-m-np, info )
348 work( 1 ) = m + np + max( lopt, int( work( m+np+1 ) ) )
349*
350 RETURN
351*
352* End of CGGGLM
353*
354 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
subroutine cggglm(n, m, p, a, lda, b, ldb, d, x, y, work, lwork, info)
CGGGLM
Definition cggglm.f:185
subroutine cggqrf(n, m, p, a, lda, taua, b, ldb, taub, work, lwork, info)
CGGQRF
Definition cggqrf.f:215
subroutine ctrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
CTRTRS
Definition ctrtrs.f:140
subroutine cunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMQR
Definition cunmqr.f:168
subroutine cunmrq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMRQ
Definition cunmrq.f:168