1*> \brief \b CTPQRT2 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 CTPQRT2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctpqrt2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctpqrt2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctpqrt2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CTPQRT2( 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   A( LDA, * ), B( LDB, * ), T( LDT, * )
28*       ..
29*
30*
31*> \par Purpose:
32*  =============
33*>
34*> \verbatim
35*>
36*> CTPQRT2 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 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 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 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*> \date September 2012
123*
124*> \ingroup complexOTHERcomputational
125*
126*> \par Further Details:
127*  =====================
128*>
129*> \verbatim
130*>
131*>  The input matrix C is a (N+M)-by-N matrix
132*>
133*>               C = [ A ]
134*>                   [ B ]
135*>
136*>  where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
137*>  matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
138*>  upper trapezoidal matrix B2:
139*>
140*>               B = [ B1 ]  <- (M-L)-by-N rectangular
141*>                   [ B2 ]  <-     L-by-N upper trapezoidal.
142*>
143*>  The upper trapezoidal matrix B2 consists of the first L rows of a
144*>  N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N).  If L=0,
145*>  B is rectangular M-by-N; if M=L=N, B is upper triangular.
146*>
147*>  The matrix W stores the elementary reflectors H(i) in the i-th column
148*>  below the diagonal (of A) in the (N+M)-by-N input matrix C
149*>
150*>               C = [ A ]  <- upper triangular N-by-N
151*>                   [ B ]  <- M-by-N pentagonal
152*>
153*>  so that W can be represented as
154*>
155*>               W = [ I ]  <- identity, N-by-N
156*>                   [ V ]  <- M-by-N, same form as B.
157*>
158*>  Thus, all of information needed for W is contained on exit in B, which
159*>  we call V above.  Note that V has the same form as B; that is,
160*>
161*>               V = [ V1 ] <- (M-L)-by-N rectangular
162*>                   [ V2 ] <-     L-by-N upper trapezoidal.
163*>
164*>  The columns of V represent the vectors which define the H(i)'s.
165*>  The (M+N)-by-(M+N) block reflector H is then given by
166*>
167*>               H = I - W * T * W**H
168*>
169*>  where W**H is the conjugate transpose of W and T is the upper triangular
170*>  factor of the block reflector.
171*> \endverbatim
172*>
173*  =====================================================================
174      SUBROUTINE CTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
175*
176*  -- LAPACK computational routine (version 3.4.2) --
177*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
178*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
179*     September 2012
180*
181*     .. Scalar Arguments ..
182      INTEGER   INFO, LDA, LDB, LDT, N, M, L
183*     ..
184*     .. Array Arguments ..
185      COMPLEX   A( LDA, * ), B( LDB, * ), T( LDT, * )
186*     ..
187*
188*  =====================================================================
189*
190*     .. Parameters ..
191      COMPLEX  ONE, ZERO
192      PARAMETER( ONE = (1.0,0.0), ZERO = (0.0,0.0) )
193*     ..
194*     .. Local Scalars ..
195      INTEGER   I, J, P, MP, NP
196      COMPLEX   ALPHA
197*     ..
198*     .. External Subroutines ..
199      EXTERNAL  CLARFG, CGEMV, CGERC, CTRMV, XERBLA
200*     ..
201*     .. Intrinsic Functions ..
202      INTRINSIC MAX, MIN
203*     ..
204*     .. Executable Statements ..
205*
206*     Test the input arguments
207*
208      INFO = 0
209      IF( M.LT.0 ) THEN
210         INFO = -1
211      ELSE IF( N.LT.0 ) THEN
212         INFO = -2
213      ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
214         INFO = -3
215      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
216         INFO = -5
217      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
218         INFO = -7
219      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
220         INFO = -9
221      END IF
222      IF( INFO.NE.0 ) THEN
223         CALL XERBLA( 'CTPQRT2', -INFO )
224         RETURN
225      END IF
226*
227*     Quick return if possible
228*
229      IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
230*
231      DO I = 1, N
232*
233*        Generate elementary reflector H(I) to annihilate B(:,I)
234*
235         P = M-L+MIN( L, I )
236         CALL CLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) )
237         IF( I.LT.N ) THEN
238*
239*           W(1:N-I) := C(I:M,I+1:N)**H * C(I:M,I) [use W = T(:,N)]
240*
241            DO J = 1, N-I
242               T( J, N ) = CONJG(A( I, I+J ))
243            END DO
244            CALL CGEMV( 'C', P, N-I, ONE, B( 1, I+1 ), LDB,
245     $                  B( 1, I ), 1, ONE, T( 1, N ), 1 )
246*
247*           C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)**H
248*
249            ALPHA = -CONJG(T( I, 1 ))
250            DO J = 1, N-I
251               A( I, I+J ) = A( I, I+J ) + ALPHA*CONJG(T( J, N ))
252            END DO
253            CALL CGERC( P, N-I, ALPHA, B( 1, I ), 1,
254     $           T( 1, N ), 1, B( 1, I+1 ), LDB )
255         END IF
256      END DO
257*
258      DO I = 2, N
259*
260*        T(1:I-1,I) := C(I:M,1:I-1)**H * (alpha * C(I:M,I))
261*
262         ALPHA = -T( I, 1 )
263
264         DO J = 1, I-1
265            T( J, I ) = ZERO
266         END DO
267         P = MIN( I-1, L )
268         MP = MIN( M-L+1, M )
269         NP = MIN( P+1, N )
270*
271*        Triangular part of B2
272*
273         DO J = 1, P
274            T( J, I ) = ALPHA*B( M-L+J, I )
275         END DO
276         CALL CTRMV( 'U', 'C', 'N', P, B( MP, 1 ), LDB,
277     $               T( 1, I ), 1 )
278*
279*        Rectangular part of B2
280*
281         CALL CGEMV( 'C', L, I-1-P, ALPHA, B( MP, NP ), LDB,
282     $               B( MP, I ), 1, ZERO, T( NP, I ), 1 )
283*
284*        B1
285*
286         CALL CGEMV( 'C', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1,
287     $               ONE, T( 1, I ), 1 )
288*
289*        T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
290*
291         CALL CTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
292*
293*        T(I,I) = tau(I)
294*
295         T( I, I ) = T( I, 1 )
296         T( I, 1 ) = ZERO
297      END DO
298
299*
300*     End of CTPQRT2
301*
302      END
303