zlarft.f revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba
1*> \brief \b ZLARFT 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZLARFT + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarft.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarft.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarft.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) 22* 23* .. Scalar Arguments .. 24* CHARACTER DIRECT, STOREV 25* INTEGER K, LDT, LDV, N 26* .. 27* .. Array Arguments .. 28* COMPLEX*16 T( LDT, * ), TAU( * ), V( LDV, * ) 29* .. 30* 31* 32*> \par Purpose: 33* ============= 34*> 35*> \verbatim 36*> 37*> ZLARFT forms the triangular factor T of a complex block reflector H 38*> of order n, which is defined as a product of k elementary reflectors. 39*> 40*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; 41*> 42*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. 43*> 44*> If STOREV = 'C', the vector which defines the elementary reflector 45*> H(i) is stored in the i-th column of the array V, and 46*> 47*> H = I - V * T * V**H 48*> 49*> If STOREV = 'R', the vector which defines the elementary reflector 50*> H(i) is stored in the i-th row of the array V, and 51*> 52*> H = I - V**H * T * V 53*> \endverbatim 54* 55* Arguments: 56* ========== 57* 58*> \param[in] DIRECT 59*> \verbatim 60*> DIRECT is CHARACTER*1 61*> Specifies the order in which the elementary reflectors are 62*> multiplied to form the block reflector: 63*> = 'F': H = H(1) H(2) . . . H(k) (Forward) 64*> = 'B': H = H(k) . . . H(2) H(1) (Backward) 65*> \endverbatim 66*> 67*> \param[in] STOREV 68*> \verbatim 69*> STOREV is CHARACTER*1 70*> Specifies how the vectors which define the elementary 71*> reflectors are stored (see also Further Details): 72*> = 'C': columnwise 73*> = 'R': rowwise 74*> \endverbatim 75*> 76*> \param[in] N 77*> \verbatim 78*> N is INTEGER 79*> The order of the block reflector H. N >= 0. 80*> \endverbatim 81*> 82*> \param[in] K 83*> \verbatim 84*> K is INTEGER 85*> The order of the triangular factor T (= the number of 86*> elementary reflectors). K >= 1. 87*> \endverbatim 88*> 89*> \param[in] V 90*> \verbatim 91*> V is COMPLEX*16 array, dimension 92*> (LDV,K) if STOREV = 'C' 93*> (LDV,N) if STOREV = 'R' 94*> The matrix V. See further details. 95*> \endverbatim 96*> 97*> \param[in] LDV 98*> \verbatim 99*> LDV is INTEGER 100*> The leading dimension of the array V. 101*> If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. 102*> \endverbatim 103*> 104*> \param[in] TAU 105*> \verbatim 106*> TAU is COMPLEX*16 array, dimension (K) 107*> TAU(i) must contain the scalar factor of the elementary 108*> reflector H(i). 109*> \endverbatim 110*> 111*> \param[out] T 112*> \verbatim 113*> T is COMPLEX*16 array, dimension (LDT,K) 114*> The k by k triangular factor T of the block reflector. 115*> If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is 116*> lower triangular. The rest of the array is not used. 117*> \endverbatim 118*> 119*> \param[in] LDT 120*> \verbatim 121*> LDT is INTEGER 122*> The leading dimension of the array T. LDT >= K. 123*> \endverbatim 124* 125* Authors: 126* ======== 127* 128*> \author Univ. of Tennessee 129*> \author Univ. of California Berkeley 130*> \author Univ. of Colorado Denver 131*> \author NAG Ltd. 132* 133*> \date April 2012 134* 135*> \ingroup complex16OTHERauxiliary 136* 137*> \par Further Details: 138* ===================== 139*> 140*> \verbatim 141*> 142*> The shape of the matrix V and the storage of the vectors which define 143*> the H(i) is best illustrated by the following example with n = 5 and 144*> k = 3. The elements equal to 1 are not stored. 145*> 146*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': 147*> 148*> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) 149*> ( v1 1 ) ( 1 v2 v2 v2 ) 150*> ( v1 v2 1 ) ( 1 v3 v3 ) 151*> ( v1 v2 v3 ) 152*> ( v1 v2 v3 ) 153*> 154*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': 155*> 156*> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) 157*> ( v1 v2 v3 ) ( v2 v2 v2 1 ) 158*> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) 159*> ( 1 v3 ) 160*> ( 1 ) 161*> \endverbatim 162*> 163* ===================================================================== 164 SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) 165* 166* -- LAPACK auxiliary routine (version 3.4.1) -- 167* -- LAPACK is a software package provided by Univ. of Tennessee, -- 168* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 169* April 2012 170* 171* .. Scalar Arguments .. 172 CHARACTER DIRECT, STOREV 173 INTEGER K, LDT, LDV, N 174* .. 175* .. Array Arguments .. 176 COMPLEX*16 T( LDT, * ), TAU( * ), V( LDV, * ) 177* .. 178* 179* ===================================================================== 180* 181* .. Parameters .. 182 COMPLEX*16 ONE, ZERO 183 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), 184 $ ZERO = ( 0.0D+0, 0.0D+0 ) ) 185* .. 186* .. Local Scalars .. 187 INTEGER I, J, PREVLASTV, LASTV 188* .. 189* .. External Subroutines .. 190 EXTERNAL ZGEMV, ZLACGV, ZTRMV 191* .. 192* .. External Functions .. 193 LOGICAL LSAME 194 EXTERNAL LSAME 195* .. 196* .. Executable Statements .. 197* 198* Quick return if possible 199* 200 IF( N.EQ.0 ) 201 $ RETURN 202* 203 IF( LSAME( DIRECT, 'F' ) ) THEN 204 PREVLASTV = N 205 DO I = 1, K 206 PREVLASTV = MAX( PREVLASTV, I ) 207 IF( TAU( I ).EQ.ZERO ) THEN 208* 209* H(i) = I 210* 211 DO J = 1, I 212 T( J, I ) = ZERO 213 END DO 214 ELSE 215* 216* general case 217* 218 IF( LSAME( STOREV, 'C' ) ) THEN 219* Skip any trailing zeros. 220 DO LASTV = N, I+1, -1 221 IF( V( LASTV, I ).NE.ZERO ) EXIT 222 END DO 223 DO J = 1, I-1 224 T( J, I ) = -TAU( I ) * CONJG( V( I , J ) ) 225 END DO 226 J = MIN( LASTV, PREVLASTV ) 227* 228* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i) 229* 230 CALL ZGEMV( 'Conjugate transpose', J-I, I-1, 231 $ -TAU( I ), V( I+1, 1 ), LDV, 232 $ V( I+1, I ), 1, ONE, T( 1, I ), 1 ) 233 ELSE 234* Skip any trailing zeros. 235 DO LASTV = N, I+1, -1 236 IF( V( I, LASTV ).NE.ZERO ) EXIT 237 END DO 238 DO J = 1, I-1 239 T( J, I ) = -TAU( I ) * V( J , I ) 240 END DO 241 J = MIN( LASTV, PREVLASTV ) 242* 243* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H 244* 245 CALL ZGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ), 246 $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, 247 $ ONE, T( 1, I ), LDT ) 248 END IF 249* 250* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) 251* 252 CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, 253 $ LDT, T( 1, I ), 1 ) 254 T( I, I ) = TAU( I ) 255 IF( I.GT.1 ) THEN 256 PREVLASTV = MAX( PREVLASTV, LASTV ) 257 ELSE 258 PREVLASTV = LASTV 259 END IF 260 END IF 261 END DO 262 ELSE 263 PREVLASTV = 1 264 DO I = K, 1, -1 265 IF( TAU( I ).EQ.ZERO ) THEN 266* 267* H(i) = I 268* 269 DO J = I, K 270 T( J, I ) = ZERO 271 END DO 272 ELSE 273* 274* general case 275* 276 IF( I.LT.K ) THEN 277 IF( LSAME( STOREV, 'C' ) ) THEN 278* Skip any leading zeros. 279 DO LASTV = 1, I-1 280 IF( V( LASTV, I ).NE.ZERO ) EXIT 281 END DO 282 DO J = I+1, K 283 T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) ) 284 END DO 285 J = MAX( LASTV, PREVLASTV ) 286* 287* T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i) 288* 289 CALL ZGEMV( 'Conjugate transpose', N-K+I-J, K-I, 290 $ -TAU( I ), V( J, I+1 ), LDV, V( J, I ), 291 $ 1, ONE, T( I+1, I ), 1 ) 292 ELSE 293* Skip any leading zeros. 294 DO LASTV = 1, I-1 295 IF( V( I, LASTV ).NE.ZERO ) EXIT 296 END DO 297 DO J = I+1, K 298 T( J, I ) = -TAU( I ) * V( J, N-K+I ) 299 END DO 300 J = MAX( LASTV, PREVLASTV ) 301* 302* T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H 303* 304 CALL ZGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ), 305 $ V( I+1, J ), LDV, V( I, J ), LDV, 306 $ ONE, T( I+1, I ), LDT ) 307 END IF 308* 309* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) 310* 311 CALL ZTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, 312 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) 313 IF( I.GT.1 ) THEN 314 PREVLASTV = MIN( PREVLASTV, LASTV ) 315 ELSE 316 PREVLASTV = LASTV 317 END IF 318 END IF 319 T( I, I ) = TAU( I ) 320 END IF 321 END DO 322 END IF 323 RETURN 324* 325* End of ZLARFT 326* 327 END 328