1*> \brief \b DLARFT 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DLARFT + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarft.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarft.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarft.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DLARFT( 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* DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) 29* .. 30* 31* 32*> \par Purpose: 33* ============= 34*> 35*> \verbatim 36*> 37*> DLARFT forms the triangular factor T of a real 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**T 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**T * 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleOTHERauxiliary 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 DLARFT( 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 DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) 177* .. 178* 179* ===================================================================== 180* 181* .. Parameters .. 182 DOUBLE PRECISION ONE, ZERO 183 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 184* .. 185* .. Local Scalars .. 186 INTEGER I, J, PREVLASTV, LASTV 187* .. 188* .. External Subroutines .. 189 EXTERNAL DGEMV, DTRMV 190* .. 191* .. External Functions .. 192 LOGICAL LSAME 193 EXTERNAL LSAME 194* .. 195* .. Executable Statements .. 196* 197* Quick return if possible 198* 199 IF( N.EQ.0 ) 200 $ RETURN 201* 202 IF( LSAME( DIRECT, 'F' ) ) THEN 203 PREVLASTV = N 204 DO I = 1, K 205 PREVLASTV = MAX( I, PREVLASTV ) 206 IF( TAU( I ).EQ.ZERO ) THEN 207* 208* H(i) = I 209* 210 DO J = 1, I 211 T( J, I ) = ZERO 212 END DO 213 ELSE 214* 215* general case 216* 217 IF( LSAME( STOREV, 'C' ) ) THEN 218* Skip any trailing zeros. 219 DO LASTV = N, I+1, -1 220 IF( V( LASTV, I ).NE.ZERO ) EXIT 221 END DO 222 DO J = 1, I-1 223 T( J, I ) = -TAU( I ) * V( I , J ) 224 END DO 225 J = MIN( LASTV, PREVLASTV ) 226* 227* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i) 228* 229 CALL DGEMV( 'Transpose', J-I, I-1, -TAU( I ), 230 $ V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE, 231 $ T( 1, I ), 1 ) 232 ELSE 233* Skip any trailing zeros. 234 DO LASTV = N, I+1, -1 235 IF( V( I, LASTV ).NE.ZERO ) EXIT 236 END DO 237 DO J = 1, I-1 238 T( J, I ) = -TAU( I ) * V( J , I ) 239 END DO 240 J = MIN( LASTV, PREVLASTV ) 241* 242* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T 243* 244 CALL DGEMV( 'No transpose', I-1, J-I, -TAU( I ), 245 $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, ONE, 246 $ T( 1, I ), 1 ) 247 END IF 248* 249* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) 250* 251 CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, 252 $ LDT, T( 1, I ), 1 ) 253 T( I, I ) = TAU( I ) 254 IF( I.GT.1 ) THEN 255 PREVLASTV = MAX( PREVLASTV, LASTV ) 256 ELSE 257 PREVLASTV = LASTV 258 END IF 259 END IF 260 END DO 261 ELSE 262 PREVLASTV = 1 263 DO I = K, 1, -1 264 IF( TAU( I ).EQ.ZERO ) THEN 265* 266* H(i) = I 267* 268 DO J = I, K 269 T( J, I ) = ZERO 270 END DO 271 ELSE 272* 273* general case 274* 275 IF( I.LT.K ) THEN 276 IF( LSAME( STOREV, 'C' ) ) THEN 277* Skip any leading zeros. 278 DO LASTV = 1, I-1 279 IF( V( LASTV, I ).NE.ZERO ) EXIT 280 END DO 281 DO J = I+1, K 282 T( J, I ) = -TAU( I ) * V( N-K+I , J ) 283 END DO 284 J = MAX( LASTV, PREVLASTV ) 285* 286* T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i) 287* 288 CALL DGEMV( 'Transpose', N-K+I-J, K-I, -TAU( I ), 289 $ V( J, I+1 ), LDV, V( J, I ), 1, ONE, 290 $ T( I+1, I ), 1 ) 291 ELSE 292* Skip any leading zeros. 293 DO LASTV = 1, I-1 294 IF( V( I, LASTV ).NE.ZERO ) EXIT 295 END DO 296 DO J = I+1, K 297 T( J, I ) = -TAU( I ) * V( J, N-K+I ) 298 END DO 299 J = MAX( LASTV, PREVLASTV ) 300* 301* T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T 302* 303 CALL DGEMV( 'No transpose', K-I, N-K+I-J, 304 $ -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV, 305 $ ONE, T( I+1, I ), 1 ) 306 END IF 307* 308* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) 309* 310 CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, 311 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) 312 IF( I.GT.1 ) THEN 313 PREVLASTV = MIN( PREVLASTV, LASTV ) 314 ELSE 315 PREVLASTV = LASTV 316 END IF 317 END IF 318 T( I, I ) = TAU( I ) 319 END IF 320 END DO 321 END IF 322 RETURN 323* 324* End of DLARFT 325* 326 END 327