1 *> \brief \b CLARFT 2 * 3 * =========== DOCUMENTATION =========== 4 * 5 * Online html documentation available at 6 * http://www.netlib.org/lapack/explore-html/ 7 * 8 *> \htmlonly 9 *> Download CLARFT + dependencies 10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarft.f"> 11 *> [TGZ]</a> 12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarft.f"> 13 *> [ZIP]</a> 14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarft.f"> 15 *> [TXT]</a> 16 *> \endhtmlonly 17 * 18 * Definition: 19 * =========== 20 * 21 * SUBROUTINE CLARFT( 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 T( LDT, * ), TAU( * ), V( LDV, * ) 29 * .. 30 * 31 * 32 *> \par Purpose: 33 * ============= 34 *> 35 *> \verbatim 36 *> 37 *> CLARFT 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 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 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 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 complexOTHERauxiliary 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 CLARFT( 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 T( LDT, * ), TAU( * ), V( LDV, * ) 177 * .. 178 * 179 * ===================================================================== 180 * 181 * .. Parameters .. 182 COMPLEX ONE, ZERO 183 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), 184 $ ZERO = ( 0.0E+0, 0.0E+0 ) ) 185 * .. 186 * .. Local Scalars .. 187 INTEGER I, J, PREVLASTV, LASTV 188 * .. 189 * .. External Subroutines .. 190 EXTERNAL CGEMV, CLACGV, CTRMV 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 CGEMV( 'Conjugate transpose', J-I, I-1, 231 $ -TAU( I ), V( I+1, 1 ), LDV, 232 $ V( I+1, I ), 1, 233 $ ONE, T( 1, I ), 1 ) 234 ELSE 235 * Skip any trailing zeros. 236 DO LASTV = N, I+1, -1 237 IF( V( I, LASTV ).NE.ZERO ) EXIT 238 END DO 239 DO J = 1, I-1 240 T( J, I ) = -TAU( I ) * V( J , I ) 241 END DO 242 J = MIN( LASTV, PREVLASTV ) 243 * 244 * T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H 245 * 246 CALL CGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ), 247 $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, 248 $ ONE, T( 1, I ), LDT ) 249 END IF 250 * 251 * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) 252 * 253 CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, 254 $ LDT, T( 1, I ), 1 ) 255 T( I, I ) = TAU( I ) 256 IF( I.GT.1 ) THEN 257 PREVLASTV = MAX( PREVLASTV, LASTV ) 258 ELSE 259 PREVLASTV = LASTV 260 END IF 261 END IF 262 END DO 263 ELSE 264 PREVLASTV = 1 265 DO I = K, 1, -1 266 IF( TAU( I ).EQ.ZERO ) THEN 267 * 268 * H(i) = I 269 * 270 DO J = I, K 271 T( J, I ) = ZERO 272 END DO 273 ELSE 274 * 275 * general case 276 * 277 IF( I.LT.K ) THEN 278 IF( LSAME( STOREV, 'C' ) ) THEN 279 * Skip any leading zeros. 280 DO LASTV = 1, I-1 281 IF( V( LASTV, I ).NE.ZERO ) EXIT 282 END DO 283 DO J = I+1, K 284 T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) ) 285 END DO 286 J = MAX( LASTV, PREVLASTV ) 287 * 288 * T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i) 289 * 290 CALL CGEMV( 'Conjugate transpose', N-K+I-J, K-I, 291 $ -TAU( I ), V( J, I+1 ), LDV, V( J, I ), 292 $ 1, ONE, T( I+1, I ), 1 ) 293 ELSE 294 * Skip any leading zeros. 295 DO LASTV = 1, I-1 296 IF( V( I, LASTV ).NE.ZERO ) EXIT 297 END DO 298 DO J = I+1, K 299 T( J, I ) = -TAU( I ) * V( J, N-K+I ) 300 END DO 301 J = MAX( LASTV, PREVLASTV ) 302 * 303 * T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H 304 * 305 CALL CGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ), 306 $ V( I+1, J ), LDV, V( I, J ), LDV, 307 $ ONE, T( I+1, I ), LDT ) 308 END IF 309 * 310 * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) 311 * 312 CALL CTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, 313 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) 314 IF( I.GT.1 ) THEN 315 PREVLASTV = MIN( PREVLASTV, LASTV ) 316 ELSE 317 PREVLASTV = LASTV 318 END IF 319 END IF 320 T( I, I ) = TAU( I ) 321 END IF 322 END DO 323 END IF 324 RETURN 325 * 326 * End of CLARFT 327 * 328 END 329
