1*> \brief \b CGGQRF
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGGQRF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggqrf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggqrf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggqrf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CGGQRF( N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK,
22*                          LWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       INTEGER            INFO, LDA, LDB, LWORK, M, N, P
26*       ..
27*       .. Array Arguments ..
28*       COMPLEX            A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
29*      $                   WORK( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> CGGQRF computes a generalized QR factorization of an N-by-M matrix A
39*> and an N-by-P matrix B:
40*>
41*>             A = Q*R,        B = Q*T*Z,
42*>
43*> where Q is an N-by-N unitary matrix, Z is a P-by-P unitary matrix,
44*> and R and T assume one of the forms:
45*>
46*> if N >= M,  R = ( R11 ) M  ,   or if N < M,  R = ( R11  R12 ) N,
47*>                 (  0  ) N-M                         N   M-N
48*>                    M
49*>
50*> where R11 is upper triangular, and
51*>
52*> if N <= P,  T = ( 0  T12 ) N,   or if N > P,  T = ( T11 ) N-P,
53*>                  P-N  N                           ( T21 ) P
54*>                                                      P
55*>
56*> where T12 or T21 is upper triangular.
57*>
58*> In particular, if B is square and nonsingular, the GQR factorization
59*> of A and B implicitly gives the QR factorization of inv(B)*A:
60*>
61*>              inv(B)*A = Z**H * (inv(T)*R)
62*>
63*> where inv(B) denotes the inverse of the matrix B, and Z' denotes the
64*> conjugate transpose of matrix Z.
65*> \endverbatim
66*
67*  Arguments:
68*  ==========
69*
70*> \param[in] N
71*> \verbatim
72*>          N is INTEGER
73*>          The number of rows of the matrices A and B. N >= 0.
74*> \endverbatim
75*>
76*> \param[in] M
77*> \verbatim
78*>          M is INTEGER
79*>          The number of columns of the matrix A.  M >= 0.
80*> \endverbatim
81*>
82*> \param[in] P
83*> \verbatim
84*>          P is INTEGER
85*>          The number of columns of the matrix B.  P >= 0.
86*> \endverbatim
87*>
88*> \param[in,out] A
89*> \verbatim
90*>          A is COMPLEX array, dimension (LDA,M)
91*>          On entry, the N-by-M matrix A.
92*>          On exit, the elements on and above the diagonal of the array
93*>          contain the min(N,M)-by-M upper trapezoidal matrix R (R is
94*>          upper triangular if N >= M); the elements below the diagonal,
95*>          with the array TAUA, represent the unitary matrix Q as a
96*>          product of min(N,M) elementary reflectors (see Further
97*>          Details).
98*> \endverbatim
99*>
100*> \param[in] LDA
101*> \verbatim
102*>          LDA is INTEGER
103*>          The leading dimension of the array A. LDA >= max(1,N).
104*> \endverbatim
105*>
106*> \param[out] TAUA
107*> \verbatim
108*>          TAUA is COMPLEX array, dimension (min(N,M))
109*>          The scalar factors of the elementary reflectors which
110*>          represent the unitary matrix Q (see Further Details).
111*> \endverbatim
112*>
113*> \param[in,out] B
114*> \verbatim
115*>          B is COMPLEX array, dimension (LDB,P)
116*>          On entry, the N-by-P matrix B.
117*>          On exit, if N <= P, the upper triangle of the subarray
118*>          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
119*>          if N > P, the elements on and above the (N-P)-th subdiagonal
120*>          contain the N-by-P upper trapezoidal matrix T; the remaining
121*>          elements, with the array TAUB, represent the unitary
122*>          matrix Z as a product of elementary reflectors (see Further
123*>          Details).
124*> \endverbatim
125*>
126*> \param[in] LDB
127*> \verbatim
128*>          LDB is INTEGER
129*>          The leading dimension of the array B. LDB >= max(1,N).
130*> \endverbatim
131*>
132*> \param[out] TAUB
133*> \verbatim
134*>          TAUB is COMPLEX array, dimension (min(N,P))
135*>          The scalar factors of the elementary reflectors which
136*>          represent the unitary matrix Z (see Further Details).
137*> \endverbatim
138*>
139*> \param[out] WORK
140*> \verbatim
141*>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
142*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
143*> \endverbatim
144*>
145*> \param[in] LWORK
146*> \verbatim
147*>          LWORK is INTEGER
148*>          The dimension of the array WORK. LWORK >= max(1,N,M,P).
149*>          For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3),
150*>          where NB1 is the optimal blocksize for the QR factorization
151*>          of an N-by-M matrix, NB2 is the optimal blocksize for the
152*>          RQ factorization of an N-by-P matrix, and NB3 is the optimal
153*>          blocksize for a call of CUNMQR.
154*>
155*>          If LWORK = -1, then a workspace query is assumed; the routine
156*>          only calculates the optimal size of the WORK array, returns
157*>          this value as the first entry of the WORK array, and no error
158*>          message related to LWORK is issued by XERBLA.
159*> \endverbatim
160*>
161*> \param[out] INFO
162*> \verbatim
163*>          INFO is INTEGER
164*>           = 0:  successful exit
165*>           < 0:  if INFO = -i, the i-th argument had an illegal value.
166*> \endverbatim
167*
168*  Authors:
169*  ========
170*
171*> \author Univ. of Tennessee
172*> \author Univ. of California Berkeley
173*> \author Univ. of Colorado Denver
174*> \author NAG Ltd.
175*
176*> \date November 2011
177*
178*> \ingroup complexOTHERcomputational
179*
180*> \par Further Details:
181*  =====================
182*>
183*> \verbatim
184*>
185*>  The matrix Q is represented as a product of elementary reflectors
186*>
187*>     Q = H(1) H(2) . . . H(k), where k = min(n,m).
188*>
189*>  Each H(i) has the form
190*>
191*>     H(i) = I - taua * v * v**H
192*>
193*>  where taua is a complex scalar, and v is a complex vector with
194*>  v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
195*>  and taua in TAUA(i).
196*>  To form Q explicitly, use LAPACK subroutine CUNGQR.
197*>  To use Q to update another matrix, use LAPACK subroutine CUNMQR.
198*>
199*>  The matrix Z is represented as a product of elementary reflectors
200*>
201*>     Z = H(1) H(2) . . . H(k), where k = min(n,p).
202*>
203*>  Each H(i) has the form
204*>
205*>     H(i) = I - taub * v * v**H
206*>
207*>  where taub is a complex scalar, and v is a complex vector with
208*>  v(p-k+i+1:p) = 0 and v(p-k+i) = 1; v(1:p-k+i-1) is stored on exit in
209*>  B(n-k+i,1:p-k+i-1), and taub in TAUB(i).
210*>  To form Z explicitly, use LAPACK subroutine CUNGRQ.
211*>  To use Z to update another matrix, use LAPACK subroutine CUNMRQ.
212*> \endverbatim
213*>
214*  =====================================================================
215      SUBROUTINE CGGQRF( N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK,
216     $                   LWORK, INFO )
217*
218*  -- LAPACK computational routine (version 3.4.0) --
219*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
220*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
221*     November 2011
222*
223*     .. Scalar Arguments ..
224      INTEGER            INFO, LDA, LDB, LWORK, M, N, P
225*     ..
226*     .. Array Arguments ..
227      COMPLEX            A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
228     $                   WORK( * )
229*     ..
230*
231*  =====================================================================
232*
233*     .. Local Scalars ..
234      LOGICAL            LQUERY
235      INTEGER            LOPT, LWKOPT, NB, NB1, NB2, NB3
236*     ..
237*     .. External Subroutines ..
238      EXTERNAL           CGEQRF, CGERQF, CUNMQR, XERBLA
239*     ..
240*     .. External Functions ..
241      INTEGER            ILAENV
242      EXTERNAL           ILAENV
243*     ..
244*     .. Intrinsic Functions ..
245      INTRINSIC          INT, MAX, MIN
246*     ..
247*     .. Executable Statements ..
248*
249*     Test the input parameters
250*
251      INFO = 0
252      NB1 = ILAENV( 1, 'CGEQRF', ' ', N, M, -1, -1 )
253      NB2 = ILAENV( 1, 'CGERQF', ' ', N, P, -1, -1 )
254      NB3 = ILAENV( 1, 'CUNMQR', ' ', N, M, P, -1 )
255      NB = MAX( NB1, NB2, NB3 )
256      LWKOPT = MAX( N, M, P)*NB
257      WORK( 1 ) = LWKOPT
258      LQUERY = ( LWORK.EQ.-1 )
259      IF( N.LT.0 ) THEN
260         INFO = -1
261      ELSE IF( M.LT.0 ) THEN
262         INFO = -2
263      ELSE IF( P.LT.0 ) THEN
264         INFO = -3
265      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
266         INFO = -5
267      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
268         INFO = -8
269      ELSE IF( LWORK.LT.MAX( 1, N, M, P ) .AND. .NOT.LQUERY ) THEN
270         INFO = -11
271      END IF
272      IF( INFO.NE.0 ) THEN
273         CALL XERBLA( 'CGGQRF', -INFO )
274         RETURN
275      ELSE IF( LQUERY ) THEN
276         RETURN
277      END IF
278*
279*     QR factorization of N-by-M matrix A: A = Q*R
280*
281      CALL CGEQRF( N, M, A, LDA, TAUA, WORK, LWORK, INFO )
282      LOPT = WORK( 1 )
283*
284*     Update B := Q**H*B.
285*
286      CALL CUNMQR( 'Left', 'Conjugate Transpose', N, P, MIN( N, M ), A,
287     $             LDA, TAUA, B, LDB, WORK, LWORK, INFO )
288      LOPT = MAX( LOPT, INT( WORK( 1 ) ) )
289*
290*     RQ factorization of N-by-P matrix B: B = T*Z.
291*
292      CALL CGERQF( N, P, B, LDB, TAUB, WORK, LWORK, INFO )
293      WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) )
294*
295      RETURN
296*
297*     End of CGGQRF
298*
299      END
300