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