1*> \brief <b> CGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact WY representation of Q. </b>
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGEQRT3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgeqrt3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgeqrt3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgeqrt3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       RECURSIVE SUBROUTINE CGEQRT3( M, N, A, LDA, T, LDT, INFO )
22*
23*       .. Scalar Arguments ..
24*       INTEGER   INFO, LDA, M, N, LDT
25*       ..
26*       .. Array Arguments ..
27*       COMPLEX   A( LDA, * ), T( LDT, * )
28*       ..
29*
30*
31*> \par Purpose:
32*  =============
33*>
34*> \verbatim
35*>
36*> CGEQRT3 recursively computes a QR factorization of a complex M-by-N matrix A,
37*> using the compact WY representation of Q.
38*>
39*> Based on the algorithm of Elmroth and Gustavson,
40*> IBM J. Res. Develop. Vol 44 No. 4 July 2000.
41*> \endverbatim
42*
43*  Arguments:
44*  ==========
45*
46*> \param[in] M
47*> \verbatim
48*>          M is INTEGER
49*>          The number of rows of the matrix A.  M >= N.
50*> \endverbatim
51*>
52*> \param[in] N
53*> \verbatim
54*>          N is INTEGER
55*>          The number of columns of the matrix A.  N >= 0.
56*> \endverbatim
57*>
58*> \param[in,out] A
59*> \verbatim
60*>          A is COMPLEX array, dimension (LDA,N)
61*>          On entry, the complex M-by-N matrix A.  On exit, the elements on and
62*>          above the diagonal contain the N-by-N upper triangular matrix R; the
63*>          elements below the diagonal are the columns of V.  See below for
64*>          further details.
65*> \endverbatim
66*>
67*> \param[in] LDA
68*> \verbatim
69*>          LDA is INTEGER
70*>          The leading dimension of the array A.  LDA >= max(1,M).
71*> \endverbatim
72*>
73*> \param[out] T
74*> \verbatim
75*>          T is COMPLEX array, dimension (LDT,N)
76*>          The N-by-N upper triangular factor of the block reflector.
77*>          The elements on and above the diagonal contain the block
78*>          reflector T; the elements below the diagonal are not used.
79*>          See below for further details.
80*> \endverbatim
81*>
82*> \param[in] LDT
83*> \verbatim
84*>          LDT is INTEGER
85*>          The leading dimension of the array T.  LDT >= max(1,N).
86*> \endverbatim
87*>
88*> \param[out] INFO
89*> \verbatim
90*>          INFO is INTEGER
91*>          = 0: successful exit
92*>          < 0: if INFO = -i, the i-th argument had an illegal value
93*> \endverbatim
94*
95*  Authors:
96*  ========
97*
98*> \author Univ. of Tennessee
99*> \author Univ. of California Berkeley
100*> \author Univ. of Colorado Denver
101*> \author NAG Ltd.
102*
103*> \ingroup complexGEcomputational
104*
105*> \par Further Details:
106*  =====================
107*>
108*> \verbatim
109*>
110*>  The matrix V stores the elementary reflectors H(i) in the i-th column
111*>  below the diagonal. For example, if M=5 and N=3, the matrix V is
112*>
113*>               V = (  1       )
114*>                   ( v1  1    )
115*>                   ( v1 v2  1 )
116*>                   ( v1 v2 v3 )
117*>                   ( v1 v2 v3 )
118*>
119*>  where the vi's represent the vectors which define H(i), which are returned
120*>  in the matrix A.  The 1's along the diagonal of V are not stored in A.  The
121*>  block reflector H is then given by
122*>
123*>               H = I - V * T * V**H
124*>
125*>  where V**H is the conjugate transpose of V.
126*>
127*>  For details of the algorithm, see Elmroth and Gustavson (cited above).
128*> \endverbatim
129*>
130*  =====================================================================
131      RECURSIVE SUBROUTINE CGEQRT3( M, N, A, LDA, T, LDT, INFO )
132*
133*  -- LAPACK computational routine --
134*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
135*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136*
137*     .. Scalar Arguments ..
138      INTEGER   INFO, LDA, M, N, LDT
139*     ..
140*     .. Array Arguments ..
141      COMPLEX   A( LDA, * ), T( LDT, * )
142*     ..
143*
144*  =====================================================================
145*
146*     .. Parameters ..
147      COMPLEX   ONE
148      PARAMETER ( ONE = (1.0,0.0) )
149*     ..
150*     .. Local Scalars ..
151      INTEGER   I, I1, J, J1, N1, N2, IINFO
152*     ..
153*     .. External Subroutines ..
154      EXTERNAL  CLARFG, CTRMM, CGEMM, XERBLA
155*     ..
156*     .. Executable Statements ..
157*
158      INFO = 0
159      IF( N .LT. 0 ) THEN
160         INFO = -2
161      ELSE IF( M .LT. N ) THEN
162         INFO = -1
163      ELSE IF( LDA .LT. MAX( 1, M ) ) THEN
164         INFO = -4
165      ELSE IF( LDT .LT. MAX( 1, N ) ) THEN
166         INFO = -6
167      END IF
168      IF( INFO.NE.0 ) THEN
169         CALL XERBLA( 'CGEQRT3', -INFO )
170         RETURN
171      END IF
172*
173      IF( N.EQ.1 ) THEN
174*
175*        Compute Householder transform when N=1
176*
177         CALL CLARFG( M, A(1,1), A( MIN( 2, M ), 1 ), 1, T(1,1) )
178*
179      ELSE
180*
181*        Otherwise, split A into blocks...
182*
183         N1 = N/2
184         N2 = N-N1
185         J1 = MIN( N1+1, N )
186         I1 = MIN( N+1, M )
187*
188*        Compute A(1:M,1:N1) <- (Y1,R1,T1), where Q1 = I - Y1 T1 Y1**H
189*
190         CALL CGEQRT3( M, N1, A, LDA, T, LDT, IINFO )
191*
192*        Compute A(1:M,J1:N) = Q1**H A(1:M,J1:N) [workspace: T(1:N1,J1:N)]
193*
194         DO J=1,N2
195            DO I=1,N1
196               T( I, J+N1 ) = A( I, J+N1 )
197            END DO
198         END DO
199         CALL CTRMM( 'L', 'L', 'C', 'U', N1, N2, ONE,
200     &               A, LDA, T( 1, J1 ), LDT )
201*
202         CALL CGEMM( 'C', 'N', N1, N2, M-N1, ONE, A( J1, 1 ), LDA,
203     &               A( J1, J1 ), LDA, ONE, T( 1, J1 ), LDT)
204*
205         CALL CTRMM( 'L', 'U', 'C', 'N', N1, N2, ONE,
206     &               T, LDT, T( 1, J1 ), LDT )
207*
208         CALL CGEMM( 'N', 'N', M-N1, N2, N1, -ONE, A( J1, 1 ), LDA,
209     &               T( 1, J1 ), LDT, ONE, A( J1, J1 ), LDA )
210*
211         CALL CTRMM( 'L', 'L', 'N', 'U', N1, N2, ONE,
212     &               A, LDA, T( 1, J1 ), LDT )
213*
214         DO J=1,N2
215            DO I=1,N1
216               A( I, J+N1 ) = A( I, J+N1 ) - T( I, J+N1 )
217            END DO
218         END DO
219*
220*        Compute A(J1:M,J1:N) <- (Y2,R2,T2) where Q2 = I - Y2 T2 Y2**H
221*
222         CALL CGEQRT3( M-N1, N2, A( J1, J1 ), LDA,
223     &                T( J1, J1 ), LDT, IINFO )
224*
225*        Compute T3 = T(1:N1,J1:N) = -T1 Y1**H Y2 T2
226*
227         DO I=1,N1
228            DO J=1,N2
229               T( I, J+N1 ) = CONJG(A( J+N1, I ))
230            END DO
231         END DO
232*
233         CALL CTRMM( 'R', 'L', 'N', 'U', N1, N2, ONE,
234     &               A( J1, J1 ), LDA, T( 1, J1 ), LDT )
235*
236         CALL CGEMM( 'C', 'N', N1, N2, M-N, ONE, A( I1, 1 ), LDA,
237     &               A( I1, J1 ), LDA, ONE, T( 1, J1 ), LDT )
238*
239         CALL CTRMM( 'L', 'U', 'N', 'N', N1, N2, -ONE, T, LDT,
240     &               T( 1, J1 ), LDT )
241*
242         CALL CTRMM( 'R', 'U', 'N', 'N', N1, N2, ONE,
243     &               T( J1, J1 ), LDT, T( 1, J1 ), LDT )
244*
245*        Y = (Y1,Y2); R = [ R1  A(1:N1,J1:N) ];  T = [T1 T3]
246*                         [  0        R2     ]       [ 0 T2]
247*
248      END IF
249*
250      RETURN
251*
252*     End of CGEQRT3
253*
254      END
255