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
