1*> \brief \b ZGGQRF
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZGGQRF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggqrf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggqrf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggqrf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZGGQRF( 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*16         A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
29*      $                   WORK( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> ZGGQRF 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**H 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*16 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*16 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*16 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*16 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*16 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 ZUNMQR.
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*> \ingroup complex16OTHERcomputational
177*
178*> \par Further Details:
179*  =====================
180*>
181*> \verbatim
182*>
183*>  The matrix Q is represented as a product of elementary reflectors
184*>
185*>     Q = H(1) H(2) . . . H(k), where k = min(n,m).
186*>
187*>  Each H(i) has the form
188*>
189*>     H(i) = I - taua * v * v**H
190*>
191*>  where taua is a complex scalar, and v is a complex vector with
192*>  v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
193*>  and taua in TAUA(i).
194*>  To form Q explicitly, use LAPACK subroutine ZUNGQR.
195*>  To use Q to update another matrix, use LAPACK subroutine ZUNMQR.
196*>
197*>  The matrix Z is represented as a product of elementary reflectors
198*>
199*>     Z = H(1) H(2) . . . H(k), where k = min(n,p).
200*>
201*>  Each H(i) has the form
202*>
203*>     H(i) = I - taub * v * v**H
204*>
205*>  where taub is a complex scalar, and v is a complex vector with
206*>  v(p-k+i+1:p) = 0 and v(p-k+i) = 1; v(1:p-k+i-1) is stored on exit in
207*>  B(n-k+i,1:p-k+i-1), and taub in TAUB(i).
208*>  To form Z explicitly, use LAPACK subroutine ZUNGRQ.
209*>  To use Z to update another matrix, use LAPACK subroutine ZUNMRQ.
210*> \endverbatim
211*>
212*  =====================================================================
213      SUBROUTINE ZGGQRF( N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK,
214     $                   LWORK, INFO )
215*
216*  -- LAPACK computational routine --
217*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
218*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
219*
220*     .. Scalar Arguments ..
221      INTEGER            INFO, LDA, LDB, LWORK, M, N, P
222*     ..
223*     .. Array Arguments ..
224      COMPLEX*16         A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ),
225     $                   WORK( * )
226*     ..
227*
228*  =====================================================================
229*
230*     .. Local Scalars ..
231      LOGICAL            LQUERY
232      INTEGER            LOPT, LWKOPT, NB, NB1, NB2, NB3
233*     ..
234*     .. External Subroutines ..
235      EXTERNAL           XERBLA, ZGEQRF, ZGERQF, ZUNMQR
236*     ..
237*     .. External Functions ..
238      INTEGER            ILAENV
239      EXTERNAL           ILAENV
240*     ..
241*     .. Intrinsic Functions ..
242      INTRINSIC          INT, MAX, MIN
243*     ..
244*     .. Executable Statements ..
245*
246*     Test the input parameters
247*
248      INFO = 0
249      NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, M, -1, -1 )
250      NB2 = ILAENV( 1, 'ZGERQF', ' ', N, P, -1, -1 )
251      NB3 = ILAENV( 1, 'ZUNMQR', ' ', N, M, P, -1 )
252      NB = MAX( NB1, NB2, NB3 )
253      LWKOPT = MAX( N, M, P )*NB
254      WORK( 1 ) = LWKOPT
255      LQUERY = ( LWORK.EQ.-1 )
256      IF( N.LT.0 ) THEN
257         INFO = -1
258      ELSE IF( M.LT.0 ) THEN
259         INFO = -2
260      ELSE IF( P.LT.0 ) THEN
261         INFO = -3
262      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
263         INFO = -5
264      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
265         INFO = -8
266      ELSE IF( LWORK.LT.MAX( 1, N, M, P ) .AND. .NOT.LQUERY ) THEN
267         INFO = -11
268      END IF
269      IF( INFO.NE.0 ) THEN
270         CALL XERBLA( 'ZGGQRF', -INFO )
271         RETURN
272      ELSE IF( LQUERY ) THEN
273         RETURN
274      END IF
275*
276*     QR factorization of N-by-M matrix A: A = Q*R
277*
278      CALL ZGEQRF( N, M, A, LDA, TAUA, WORK, LWORK, INFO )
279      LOPT = DBLE( WORK( 1 ) )
280*
281*     Update B := Q**H*B.
282*
283      CALL ZUNMQR( 'Left', 'Conjugate Transpose', N, P, MIN( N, M ), A,
284     $             LDA, TAUA, B, LDB, WORK, LWORK, INFO )
285      LOPT = MAX( LOPT, INT( WORK( 1 ) ) )
286*
287*     RQ factorization of N-by-P matrix B: B = T*Z.
288*
289      CALL ZGERQF( N, P, B, LDB, TAUB, WORK, LWORK, INFO )
290      WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) )
291*
292      RETURN
293*
294*     End of ZGGQRF
295*
296      END
297