1 SUBROUTINE CUNGLQ( M, N, K, A, LDA, TAU, WORK, LWORK, INFO ) 2* 3* -- LAPACK routine (version 3.0) -- 4* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5* Courant Institute, Argonne National Lab, and Rice University 6* June 30, 1999 7* 8* .. Scalar Arguments .. 9 INTEGER INFO, K, LDA, LWORK, M, N 10* .. 11* .. Array Arguments .. 12 COMPLEX A( LDA, * ), TAU( * ), WORK( * ) 13* .. 14* 15* Purpose 16* ======= 17* 18* CUNGLQ generates an M-by-N complex matrix Q with orthonormal rows, 19* which is defined as the first M rows of a product of K elementary 20* reflectors of order N 21* 22* Q = H(k)' . . . H(2)' H(1)' 23* 24* as returned by CGELQF. 25* 26* Arguments 27* ========= 28* 29* M (input) INTEGER 30* The number of rows of the matrix Q. M >= 0. 31* 32* N (input) INTEGER 33* The number of columns of the matrix Q. N >= M. 34* 35* K (input) INTEGER 36* The number of elementary reflectors whose product defines the 37* matrix Q. M >= K >= 0. 38* 39* A (input/output) COMPLEX array, dimension (LDA,N) 40* On entry, the i-th row must contain the vector which defines 41* the elementary reflector H(i), for i = 1,2,...,k, as returned 42* by CGELQF in the first k rows of its array argument A. 43* On exit, the M-by-N matrix Q. 44* 45* LDA (input) INTEGER 46* The first dimension of the array A. LDA >= max(1,M). 47* 48* TAU (input) COMPLEX array, dimension (K) 49* TAU(i) must contain the scalar factor of the elementary 50* reflector H(i), as returned by CGELQF. 51* 52* WORK (workspace/output) COMPLEX array, dimension (LWORK) 53* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 54* 55* LWORK (input) INTEGER 56* The dimension of the array WORK. LWORK >= max(1,M). 57* For optimum performance LWORK >= M*NB, where NB is 58* the optimal blocksize. 59* 60* If LWORK = -1, then a workspace query is assumed; the routine 61* only calculates the optimal size of the WORK array, returns 62* this value as the first entry of the WORK array, and no error 63* message related to LWORK is issued by XERBLA. 64* 65* INFO (output) INTEGER 66* = 0: successful exit; 67* < 0: if INFO = -i, the i-th argument has an illegal value 68* 69* ===================================================================== 70* 71* .. Parameters .. 72 COMPLEX ZERO 73 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) ) 74* .. 75* .. Local Scalars .. 76 LOGICAL LQUERY 77 INTEGER I, IB, IINFO, IWS, J, KI, KK, L, LDWORK, 78 $ LWKOPT, NB, NBMIN, NX 79* .. 80* .. External Subroutines .. 81 EXTERNAL CLARFB, CLARFT, CUNGL2, XERBLA 82* .. 83* .. Intrinsic Functions .. 84 INTRINSIC MAX, MIN 85* .. 86* .. External Functions .. 87 INTEGER ILAENV 88 EXTERNAL ILAENV 89* .. 90* .. Executable Statements .. 91* 92* Test the input arguments 93* 94 INFO = 0 95 NB = ILAENV( 1, 'CUNGLQ', ' ', M, N, K, -1 ) 96 LWKOPT = MAX( 1, M )*NB 97 WORK( 1 ) = LWKOPT 98 LQUERY = ( LWORK.EQ.-1 ) 99 IF( M.LT.0 ) THEN 100 INFO = -1 101 ELSE IF( N.LT.M ) THEN 102 INFO = -2 103 ELSE IF( K.LT.0 .OR. K.GT.M ) THEN 104 INFO = -3 105 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 106 INFO = -5 107 ELSE IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN 108 INFO = -8 109 END IF 110 IF( INFO.NE.0 ) THEN 111 CALL XERBLA( 'CUNGLQ', -INFO ) 112 RETURN 113 ELSE IF( LQUERY ) THEN 114 RETURN 115 END IF 116* 117* Quick return if possible 118* 119 IF( M.LE.0 ) THEN 120 WORK( 1 ) = 1 121 RETURN 122 END IF 123* 124 NBMIN = 2 125 NX = 0 126 IWS = M 127 IF( NB.GT.1 .AND. NB.LT.K ) THEN 128* 129* Determine when to cross over from blocked to unblocked code. 130* 131 NX = MAX( 0, ILAENV( 3, 'CUNGLQ', ' ', M, N, K, -1 ) ) 132 IF( NX.LT.K ) THEN 133* 134* Determine if workspace is large enough for blocked code. 135* 136 LDWORK = M 137 IWS = LDWORK*NB 138 IF( LWORK.LT.IWS ) THEN 139* 140* Not enough workspace to use optimal NB: reduce NB and 141* determine the minimum value of NB. 142* 143 NB = LWORK / LDWORK 144 NBMIN = MAX( 2, ILAENV( 2, 'CUNGLQ', ' ', M, N, K, -1 ) ) 145 END IF 146 END IF 147 END IF 148* 149 IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN 150* 151* Use blocked code after the last block. 152* The first kk rows are handled by the block method. 153* 154 KI = ( ( K-NX-1 ) / NB )*NB 155 KK = MIN( K, KI+NB ) 156* 157* Set A(kk+1:m,1:kk) to zero. 158* 159 DO 20 J = 1, KK 160 DO 10 I = KK + 1, M 161 A( I, J ) = ZERO 162 10 CONTINUE 163 20 CONTINUE 164 ELSE 165 KK = 0 166 END IF 167* 168* Use unblocked code for the last or only block. 169* 170 IF( KK.LT.M ) 171 $ CALL CUNGL2( M-KK, N-KK, K-KK, A( KK+1, KK+1 ), LDA, 172 $ TAU( KK+1 ), WORK, IINFO ) 173* 174 IF( KK.GT.0 ) THEN 175* 176* Use blocked code 177* 178 DO 50 I = KI + 1, 1, -NB 179 IB = MIN( NB, K-I+1 ) 180 IF( I+IB.LE.M ) THEN 181* 182* Form the triangular factor of the block reflector 183* H = H(i) H(i+1) . . . H(i+ib-1) 184* 185 CALL CLARFT( 'Forward', 'Rowwise', N-I+1, IB, A( I, I ), 186 $ LDA, TAU( I ), WORK, LDWORK ) 187* 188* Apply H' to A(i+ib:m,i:n) from the right 189* 190 CALL CLARFB( 'Right', 'Conjugate transpose', 'Forward', 191 $ 'Rowwise', M-I-IB+1, N-I+1, IB, A( I, I ), 192 $ LDA, WORK, LDWORK, A( I+IB, I ), LDA, 193 $ WORK( IB+1 ), LDWORK ) 194 END IF 195* 196* Apply H' to columns i:n of current block 197* 198 CALL CUNGL2( IB, N-I+1, IB, A( I, I ), LDA, TAU( I ), WORK, 199 $ IINFO ) 200* 201* Set columns 1:i-1 of current block to zero 202* 203 DO 40 J = 1, I - 1 204 DO 30 L = I, I + IB - 1 205 A( L, J ) = ZERO 206 30 CONTINUE 207 40 CONTINUE 208 50 CONTINUE 209 END IF 210* 211 WORK( 1 ) = IWS 212 RETURN 213* 214* End of CUNGLQ 215* 216 END 217