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 complexOTHEReigen
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      EXTERNAL           ILAENV
217*     ..
218*     .. Intrinsic Functions ..
219      INTRINSIC          INT, MAX, MIN
220*     ..
221*     .. Executable Statements ..
222*
223*     Test the input parameters
224*
225      INFO = 0
226      NP = MIN( N, P )
227      LQUERY = ( LWORK.EQ.-1 )
228      IF( N.LT.0 ) THEN
229         INFO = -1
230      ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
231         INFO = -2
232      ELSE IF( P.LT.0 .OR. P.LT.N-M ) THEN
233         INFO = -3
234      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
235         INFO = -5
236      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
237         INFO = -7
238      END IF
239*
240*     Calculate workspace
241*
242      IF( INFO.EQ.0) THEN
243         IF( N.EQ.0 ) THEN
244            LWKMIN = 1
245            LWKOPT = 1
246         ELSE
247            NB1 = ILAENV( 1, 'CGEQRF', ' ', N, M, -1, -1 )
248            NB2 = ILAENV( 1, 'CGERQF', ' ', N, M, -1, -1 )
249            NB3 = ILAENV( 1, 'CUNMQR', ' ', N, M, P, -1 )
250            NB4 = ILAENV( 1, 'CUNMRQ', ' ', N, M, P, -1 )
251            NB = MAX( NB1, NB2, NB3, NB4 )
252            LWKMIN = M + N + P
253            LWKOPT = M + NP + MAX( N, P )*NB
254         END IF
255         WORK( 1 ) = LWKOPT
256*
257         IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
258            INFO = -12
259         END IF
260      END IF
261*
262      IF( INFO.NE.0 ) THEN
263         CALL XERBLA( 'CGGGLM', -INFO )
264         RETURN
265      ELSE IF( LQUERY ) THEN
266         RETURN
267      END IF
268*
269*     Quick return if possible
270*
271      IF( N.EQ.0 ) THEN
272         DO I = 1, M
273            X(I) = CZERO
274         END DO
275         DO I = 1, P
276            Y(I) = CZERO
277         END DO
278         RETURN
279      END IF
280*
281*     Compute the GQR factorization of matrices A and B:
282*
283*          Q**H*A = ( R11 ) M,    Q**H*B*Z**H = ( T11   T12 ) M
284*                   (  0  ) N-M                 (  0    T22 ) N-M
285*                      M                         M+P-N  N-M
286*
287*     where R11 and T22 are upper triangular, and Q and Z are
288*     unitary.
289*
290      CALL CGGQRF( N, M, P, A, LDA, WORK, B, LDB, WORK( M+1 ),
291     $             WORK( M+NP+1 ), LWORK-M-NP, INFO )
292      LOPT = REAL( WORK( M+NP+1 ) )
293*
294*     Update left-hand-side vector d = Q**H*d = ( d1 ) M
295*                                               ( d2 ) N-M
296*
297      CALL CUNMQR( 'Left', 'Conjugate transpose', N, 1, M, A, LDA, WORK,
298     $             D, MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
299      LOPT = MAX( LOPT, INT( WORK( M+NP+1 ) ) )
300*
301*     Solve T22*y2 = d2 for y2
302*
303      IF( N.GT.M ) THEN
304         CALL CTRTRS( 'Upper', 'No transpose', 'Non unit', N-M, 1,
305     $                B( M+1, M+P-N+1 ), LDB, D( M+1 ), N-M, INFO )
306*
307         IF( INFO.GT.0 ) THEN
308            INFO = 1
309            RETURN
310         END IF
311*
312         CALL CCOPY( N-M, D( M+1 ), 1, Y( M+P-N+1 ), 1 )
313      END IF
314*
315*     Set y1 = 0
316*
317      DO 10 I = 1, M + P - N
318         Y( I ) = CZERO
319   10 CONTINUE
320*
321*     Update d1 = d1 - T12*y2
322*
323      CALL CGEMV( 'No transpose', M, N-M, -CONE, B( 1, M+P-N+1 ), LDB,
324     $            Y( M+P-N+1 ), 1, CONE, D, 1 )
325*
326*     Solve triangular system: R11*x = d1
327*
328      IF( M.GT.0 ) THEN
329         CALL CTRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, LDA,
330     $                D, M, INFO )
331*
332         IF( INFO.GT.0 ) THEN
333            INFO = 2
334            RETURN
335         END IF
336*
337*        Copy D to X
338*
339         CALL CCOPY( M, D, 1, X, 1 )
340      END IF
341*
342*     Backward transformation y = Z**H *y
343*
344      CALL CUNMRQ( 'Left', 'Conjugate transpose', P, 1, NP,
345     $             B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
346     $             MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
347      WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) )
348*
349      RETURN
350*
351*     End of CGGGLM
352*
353      END
354