1*> \brief \b CGEQRT3 recursively computes a QR factorization of a general real or complex matrix using the compact WY representation of Q.
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*> \date September 2012
104*
105*> \ingroup complexGEcomputational
106*
107*> \par Further Details:
108*  =====================
109*>
110*> \verbatim
111*>
112*>  The matrix V stores the elementary reflectors H(i) in the i-th column
113*>  below the diagonal. For example, if M=5 and N=3, the matrix V is
114*>
115*>               V = (  1       )
116*>                   ( v1  1    )
117*>                   ( v1 v2  1 )
118*>                   ( v1 v2 v3 )
119*>                   ( v1 v2 v3 )
120*>
121*>  where the vi's represent the vectors which define H(i), which are returned
122*>  in the matrix A.  The 1's along the diagonal of V are not stored in A.  The
123*>  block reflector H is then given by
124*>
125*>               H = I - V * T * V**H
126*>
127*>  where V**H is the conjugate transpose of V.
128*>
129*>  For details of the algorithm, see Elmroth and Gustavson (cited above).
130*> \endverbatim
131*>
132*  =====================================================================
133      RECURSIVE SUBROUTINE CGEQRT3( M, N, A, LDA, T, LDT, INFO )
134*
135*  -- LAPACK computational routine (version 3.4.2) --
136*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
137*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
138*     September 2012
139*
140*     .. Scalar Arguments ..
141      INTEGER   INFO, LDA, M, N, LDT
142*     ..
143*     .. Array Arguments ..
144      COMPLEX   A( LDA, * ), T( LDT, * )
145*     ..
146*
147*  =====================================================================
148*
149*     .. Parameters ..
150      COMPLEX   ONE
151      PARAMETER ( ONE = (1.0,0.0) )
152*     ..
153*     .. Local Scalars ..
154      INTEGER   I, I1, J, J1, N1, N2, IINFO
155*     ..
156*     .. External Subroutines ..
157      EXTERNAL  CLARFG, CTRMM, CGEMM, XERBLA
158*     ..
159*     .. Executable Statements ..
160*
161      INFO = 0
162      IF( N .LT. 0 ) THEN
163         INFO = -2
164      ELSE IF( M .LT. N ) THEN
165         INFO = -1
166      ELSE IF( LDA .LT. MAX( 1, M ) ) THEN
167         INFO = -4
168      ELSE IF( LDT .LT. MAX( 1, N ) ) THEN
169         INFO = -6
170      END IF
171      IF( INFO.NE.0 ) THEN
172         CALL XERBLA( 'CGEQRT3', -INFO )
173         RETURN
174      END IF
175*
176      IF( N.EQ.1 ) THEN
177*
178*        Compute Householder transform when N=1
179*
180         CALL CLARFG( M, A, A( MIN( 2, M ), 1 ), 1, T )
181*
182      ELSE
183*
184*        Otherwise, split A into blocks...
185*
186         N1 = N/2
187         N2 = N-N1
188         J1 = MIN( N1+1, N )
189         I1 = MIN( N+1, M )
190*
191*        Compute A(1:M,1:N1) <- (Y1,R1,T1), where Q1 = I - Y1 T1 Y1**H
192*
193         CALL CGEQRT3( M, N1, A, LDA, T, LDT, IINFO )
194*
195*        Compute A(1:M,J1:N) = Q1**H A(1:M,J1:N) [workspace: T(1:N1,J1:N)]
196*
197         DO J=1,N2
198            DO I=1,N1
199               T( I, J+N1 ) = A( I, J+N1 )
200            END DO
201         END DO
202         CALL CTRMM( 'L', 'L', 'C', 'U', N1, N2, ONE,
203     &               A, LDA, T( 1, J1 ), LDT )
204*
205         CALL CGEMM( 'C', 'N', N1, N2, M-N1, ONE, A( J1, 1 ), LDA,
206     &               A( J1, J1 ), LDA, ONE, T( 1, J1 ), LDT)
207*
208         CALL CTRMM( 'L', 'U', 'C', 'N', N1, N2, ONE,
209     &               T, LDT, T( 1, J1 ), LDT )
210*
211         CALL CGEMM( 'N', 'N', M-N1, N2, N1, -ONE, A( J1, 1 ), LDA,
212     &               T( 1, J1 ), LDT, ONE, A( J1, J1 ), LDA )
213*
214         CALL CTRMM( 'L', 'L', 'N', 'U', N1, N2, ONE,
215     &               A, LDA, T( 1, J1 ), LDT )
216*
217         DO J=1,N2
218            DO I=1,N1
219               A( I, J+N1 ) = A( I, J+N1 ) - T( I, J+N1 )
220            END DO
221         END DO
222*
223*        Compute A(J1:M,J1:N) <- (Y2,R2,T2) where Q2 = I - Y2 T2 Y2**H
224*
225         CALL CGEQRT3( M-N1, N2, A( J1, J1 ), LDA,
226     &                T( J1, J1 ), LDT, IINFO )
227*
228*        Compute T3 = T(1:N1,J1:N) = -T1 Y1**H Y2 T2
229*
230         DO I=1,N1
231            DO J=1,N2
232               T( I, J+N1 ) = CONJG(A( J+N1, I ))
233            END DO
234         END DO
235*
236         CALL CTRMM( 'R', 'L', 'N', 'U', N1, N2, ONE,
237     &               A( J1, J1 ), LDA, T( 1, J1 ), LDT )
238*
239         CALL CGEMM( 'C', 'N', N1, N2, M-N, ONE, A( I1, 1 ), LDA,
240     &               A( I1, J1 ), LDA, ONE, T( 1, J1 ), LDT )
241*
242         CALL CTRMM( 'L', 'U', 'N', 'N', N1, N2, -ONE, T, LDT,
243     &               T( 1, J1 ), LDT )
244*
245         CALL CTRMM( 'R', 'U', 'N', 'N', N1, N2, ONE,
246     &               T( J1, J1 ), LDT, T( 1, J1 ), LDT )
247*
248*        Y = (Y1,Y2); R = [ R1  A(1:N1,J1:N) ];  T = [T1 T3]
249*                         [  0        R2     ]       [ 0 T2]
250*
251      END IF
252*
253      RETURN
254*
255*     End of CGEQRT3
256*
257      END
258