1*> \brief \b ZTPQRT2 computes a QR factorization of a real or complex "triangular-pentagonal" matrix, which is composed of a triangular block and a pentagonal block, using the compact WY representation for Q.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZTPQRT2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztpqrt2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztpqrt2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztpqrt2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
22*
23*       .. Scalar Arguments ..
24*       INTEGER   INFO, LDA, LDB, LDT, N, M, L
25*       ..
26*       .. Array Arguments ..
27*       COMPLEX*16   A( LDA, * ), B( LDB, * ), T( LDT, * )
28*       ..
29*
30*
31*> \par Purpose:
32*  =============
33*>
34*> \verbatim
35*>
36*> ZTPQRT2 computes a QR factorization of a complex "triangular-pentagonal"
37*> matrix C, which is composed of a triangular block A and pentagonal block B,
38*> using the compact WY representation for Q.
39*> \endverbatim
40*
41*  Arguments:
42*  ==========
43*
44*> \param[in] M
45*> \verbatim
46*>          M is INTEGER
47*>          The total number of rows of the matrix B.
48*>          M >= 0.
49*> \endverbatim
50*>
51*> \param[in] N
52*> \verbatim
53*>          N is INTEGER
54*>          The number of columns of the matrix B, and the order of
55*>          the triangular matrix A.
56*>          N >= 0.
57*> \endverbatim
58*>
59*> \param[in] L
60*> \verbatim
61*>          L is INTEGER
62*>          The number of rows of the upper trapezoidal part of B.
63*>          MIN(M,N) >= L >= 0.  See Further Details.
64*> \endverbatim
65*>
66*> \param[in,out] A
67*> \verbatim
68*>          A is COMPLEX*16 array, dimension (LDA,N)
69*>          On entry, the upper triangular N-by-N matrix A.
70*>          On exit, the elements on and above the diagonal of the array
71*>          contain the upper triangular matrix R.
72*> \endverbatim
73*>
74*> \param[in] LDA
75*> \verbatim
76*>          LDA is INTEGER
77*>          The leading dimension of the array A.  LDA >= max(1,N).
78*> \endverbatim
79*>
80*> \param[in,out] B
81*> \verbatim
82*>          B is COMPLEX*16 array, dimension (LDB,N)
83*>          On entry, the pentagonal M-by-N matrix B.  The first M-L rows
84*>          are rectangular, and the last L rows are upper trapezoidal.
85*>          On exit, B contains the pentagonal matrix V.  See Further Details.
86*> \endverbatim
87*>
88*> \param[in] LDB
89*> \verbatim
90*>          LDB is INTEGER
91*>          The leading dimension of the array B.  LDB >= max(1,M).
92*> \endverbatim
93*>
94*> \param[out] T
95*> \verbatim
96*>          T is COMPLEX*16 array, dimension (LDT,N)
97*>          The N-by-N upper triangular factor T of the block reflector.
98*>          See Further Details.
99*> \endverbatim
100*>
101*> \param[in] LDT
102*> \verbatim
103*>          LDT is INTEGER
104*>          The leading dimension of the array T.  LDT >= max(1,N)
105*> \endverbatim
106*>
107*> \param[out] INFO
108*> \verbatim
109*>          INFO is INTEGER
110*>          = 0: successful exit
111*>          < 0: if INFO = -i, the i-th argument had an illegal value
112*> \endverbatim
113*
114*  Authors:
115*  ========
116*
117*> \author Univ. of Tennessee
118*> \author Univ. of California Berkeley
119*> \author Univ. of Colorado Denver
120*> \author NAG Ltd.
121*
122*> \ingroup complex16OTHERcomputational
123*
124*> \par Further Details:
125*  =====================
126*>
127*> \verbatim
128*>
129*>  The input matrix C is a (N+M)-by-N matrix
130*>
131*>               C = [ A ]
132*>                   [ B ]
133*>
134*>  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
135*>  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
136*>  upper trapezoidal matrix B2:
137*>
138*>               B = [ B1 ]  <- (M-L)-by-N rectangular
139*>                   [ B2 ]  <-     L-by-N upper trapezoidal.
140*>
141*>  The upper trapezoidal matrix B2 consists of the first L rows of a
142*>  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0,
143*>  B is rectangular M-by-N; if M=L=N, B is upper triangular.
144*>
145*>  The matrix W stores the elementary reflectors H(i) in the i-th column
146*>  below the diagonal (of A) in the (N+M)-by-N input matrix C
147*>
148*>               C = [ A ]  <- upper triangular N-by-N
149*>                   [ B ]  <- M-by-N pentagonal
150*>
151*>  so that W can be represented as
152*>
153*>               W = [ I ]  <- identity, N-by-N
154*>                   [ V ]  <- M-by-N, same form as B.
155*>
156*>  Thus, all of information needed for W is contained on exit in B, which
157*>  we call V above.  Note that V has the same form as B; that is,
158*>
159*>               V = [ V1 ] <- (M-L)-by-N rectangular
160*>                   [ V2 ] <-     L-by-N upper trapezoidal.
161*>
162*>  The columns of V represent the vectors which define the H(i)'s.
163*>  The (M+N)-by-(M+N) block reflector H is then given by
164*>
165*>               H = I - W * T * W**H
166*>
167*>  where W**H is the conjugate transpose of W and T is the upper triangular
168*>  factor of the block reflector.
169*> \endverbatim
170*>
171*  =====================================================================
172      SUBROUTINE ZTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
173*
174*  -- LAPACK computational routine --
175*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
176*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
177*
178*     .. Scalar Arguments ..
179      INTEGER   INFO, LDA, LDB, LDT, N, M, L
180*     ..
181*     .. Array Arguments ..
182      COMPLEX*16   A( LDA, * ), B( LDB, * ), T( LDT, * )
183*     ..
184*
185*  =====================================================================
186*
187*     .. Parameters ..
188      COMPLEX*16  ONE, ZERO
189      PARAMETER( ONE = (1.0,0.0), ZERO = (0.0,0.0) )
190*     ..
191*     .. Local Scalars ..
192      INTEGER   I, J, P, MP, NP
193      COMPLEX*16   ALPHA
194*     ..
195*     .. External Subroutines ..
196      EXTERNAL  ZLARFG, ZGEMV, ZGERC, ZTRMV, XERBLA
197*     ..
198*     .. Intrinsic Functions ..
199      INTRINSIC MAX, MIN
200*     ..
201*     .. Executable Statements ..
202*
203*     Test the input arguments
204*
205      INFO = 0
206      IF( M.LT.0 ) THEN
207         INFO = -1
208      ELSE IF( N.LT.0 ) THEN
209         INFO = -2
210      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
211         INFO = -3
212      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
213         INFO = -5
214      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
215         INFO = -7
216      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
217         INFO = -9
218      END IF
219      IF( INFO.NE.0 ) THEN
220         CALL XERBLA( 'ZTPQRT2', -INFO )
221         RETURN
222      END IF
223*
224*     Quick return if possible
225*
226      IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
227*
228      DO I = 1, N
229*
230*        Generate elementary reflector H(I) to annihilate B(:,I)
231*
232         P = M-L+MIN( L, I )
233         CALL ZLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) )
234         IF( I.LT.N ) THEN
235*
236*           W(1:N-I) := C(I:M,I+1:N)**H * C(I:M,I) [use W = T(:,N)]
237*
238            DO J = 1, N-I
239               T( J, N ) = CONJG(A( I, I+J ))
240            END DO
241            CALL ZGEMV( 'C', P, N-I, ONE, B( 1, I+1 ), LDB,
242     $                  B( 1, I ), 1, ONE, T( 1, N ), 1 )
243*
244*           C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)**H
245*
246            ALPHA = -CONJG(T( I, 1 ))
247            DO J = 1, N-I
248               A( I, I+J ) = A( I, I+J ) + ALPHA*CONJG(T( J, N ))
249            END DO
250            CALL ZGERC( P, N-I, ALPHA, B( 1, I ), 1,
251     $           T( 1, N ), 1, B( 1, I+1 ), LDB )
252         END IF
253      END DO
254*
255      DO I = 2, N
256*
257*        T(1:I-1,I) := C(I:M,1:I-1)**H * (alpha * C(I:M,I))
258*
259         ALPHA = -T( I, 1 )
260
261         DO J = 1, I-1
262            T( J, I ) = ZERO
263         END DO
264         P = MIN( I-1, L )
265         MP = MIN( M-L+1, M )
266         NP = MIN( P+1, N )
267*
268*        Triangular part of B2
269*
270         DO J = 1, P
271            T( J, I ) = ALPHA*B( M-L+J, I )
272         END DO
273         CALL ZTRMV( 'U', 'C', 'N', P, B( MP, 1 ), LDB,
274     $               T( 1, I ), 1 )
275*
276*        Rectangular part of B2
277*
278         CALL ZGEMV( 'C', L, I-1-P, ALPHA, B( MP, NP ), LDB,
279     $               B( MP, I ), 1, ZERO, T( NP, I ), 1 )
280*
281*        B1
282*
283         CALL ZGEMV( 'C', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1,
284     $               ONE, T( 1, I ), 1 )
285*
286*        T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
287*
288         CALL ZTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
289*
290*        T(I,I) = tau(I)
291*
292         T( I, I ) = T( I, 1 )
293         T( I, 1 ) = ZERO
294      END DO
295
296*
297*     End of ZTPQRT2
298*
299      END
300