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
