[136] | 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
|
---|