![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZLARFT 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZLARFT + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarft.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarft.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarft.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER DIRECT, STOREV 00025 * INTEGER K, LDT, LDV, N 00026 * .. 00027 * .. Array Arguments .. 00028 * COMPLEX*16 T( LDT, * ), TAU( * ), V( LDV, * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> ZLARFT forms the triangular factor T of a complex block reflector H 00038 *> of order n, which is defined as a product of k elementary reflectors. 00039 *> 00040 *> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; 00041 *> 00042 *> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. 00043 *> 00044 *> If STOREV = 'C', the vector which defines the elementary reflector 00045 *> H(i) is stored in the i-th column of the array V, and 00046 *> 00047 *> H = I - V * T * V**H 00048 *> 00049 *> If STOREV = 'R', the vector which defines the elementary reflector 00050 *> H(i) is stored in the i-th row of the array V, and 00051 *> 00052 *> H = I - V**H * T * V 00053 *> \endverbatim 00054 * 00055 * Arguments: 00056 * ========== 00057 * 00058 *> \param[in] DIRECT 00059 *> \verbatim 00060 *> DIRECT is CHARACTER*1 00061 *> Specifies the order in which the elementary reflectors are 00062 *> multiplied to form the block reflector: 00063 *> = 'F': H = H(1) H(2) . . . H(k) (Forward) 00064 *> = 'B': H = H(k) . . . H(2) H(1) (Backward) 00065 *> \endverbatim 00066 *> 00067 *> \param[in] STOREV 00068 *> \verbatim 00069 *> STOREV is CHARACTER*1 00070 *> Specifies how the vectors which define the elementary 00071 *> reflectors are stored (see also Further Details): 00072 *> = 'C': columnwise 00073 *> = 'R': rowwise 00074 *> \endverbatim 00075 *> 00076 *> \param[in] N 00077 *> \verbatim 00078 *> N is INTEGER 00079 *> The order of the block reflector H. N >= 0. 00080 *> \endverbatim 00081 *> 00082 *> \param[in] K 00083 *> \verbatim 00084 *> K is INTEGER 00085 *> The order of the triangular factor T (= the number of 00086 *> elementary reflectors). K >= 1. 00087 *> \endverbatim 00088 *> 00089 *> \param[in] V 00090 *> \verbatim 00091 *> V is COMPLEX*16 array, dimension 00092 *> (LDV,K) if STOREV = 'C' 00093 *> (LDV,N) if STOREV = 'R' 00094 *> The matrix V. See further details. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] LDV 00098 *> \verbatim 00099 *> LDV is INTEGER 00100 *> The leading dimension of the array V. 00101 *> If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. 00102 *> \endverbatim 00103 *> 00104 *> \param[in] TAU 00105 *> \verbatim 00106 *> TAU is COMPLEX*16 array, dimension (K) 00107 *> TAU(i) must contain the scalar factor of the elementary 00108 *> reflector H(i). 00109 *> \endverbatim 00110 *> 00111 *> \param[out] T 00112 *> \verbatim 00113 *> T is COMPLEX*16 array, dimension (LDT,K) 00114 *> The k by k triangular factor T of the block reflector. 00115 *> If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is 00116 *> lower triangular. The rest of the array is not used. 00117 *> \endverbatim 00118 *> 00119 *> \param[in] LDT 00120 *> \verbatim 00121 *> LDT is INTEGER 00122 *> The leading dimension of the array T. LDT >= K. 00123 *> \endverbatim 00124 * 00125 * Authors: 00126 * ======== 00127 * 00128 *> \author Univ. of Tennessee 00129 *> \author Univ. of California Berkeley 00130 *> \author Univ. of Colorado Denver 00131 *> \author NAG Ltd. 00132 * 00133 *> \date April 2012 00134 * 00135 *> \ingroup complex16OTHERauxiliary 00136 * 00137 *> \par Further Details: 00138 * ===================== 00139 *> 00140 *> \verbatim 00141 *> 00142 *> The shape of the matrix V and the storage of the vectors which define 00143 *> the H(i) is best illustrated by the following example with n = 5 and 00144 *> k = 3. The elements equal to 1 are not stored. 00145 *> 00146 *> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': 00147 *> 00148 *> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) 00149 *> ( v1 1 ) ( 1 v2 v2 v2 ) 00150 *> ( v1 v2 1 ) ( 1 v3 v3 ) 00151 *> ( v1 v2 v3 ) 00152 *> ( v1 v2 v3 ) 00153 *> 00154 *> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': 00155 *> 00156 *> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) 00157 *> ( v1 v2 v3 ) ( v2 v2 v2 1 ) 00158 *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) 00159 *> ( 1 v3 ) 00160 *> ( 1 ) 00161 *> \endverbatim 00162 *> 00163 * ===================================================================== 00164 SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) 00165 * 00166 * -- LAPACK auxiliary routine (version 3.4.1) -- 00167 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00168 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00169 * April 2012 00170 * 00171 * .. Scalar Arguments .. 00172 CHARACTER DIRECT, STOREV 00173 INTEGER K, LDT, LDV, N 00174 * .. 00175 * .. Array Arguments .. 00176 COMPLEX*16 T( LDT, * ), TAU( * ), V( LDV, * ) 00177 * .. 00178 * 00179 * ===================================================================== 00180 * 00181 * .. Parameters .. 00182 COMPLEX*16 ONE, ZERO 00183 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), 00184 $ ZERO = ( 0.0D+0, 0.0D+0 ) ) 00185 * .. 00186 * .. Local Scalars .. 00187 INTEGER I, J, PREVLASTV, LASTV 00188 * .. 00189 * .. External Subroutines .. 00190 EXTERNAL ZGEMV, ZLACGV, ZTRMV 00191 * .. 00192 * .. External Functions .. 00193 LOGICAL LSAME 00194 EXTERNAL LSAME 00195 * .. 00196 * .. Executable Statements .. 00197 * 00198 * Quick return if possible 00199 * 00200 IF( N.EQ.0 ) 00201 $ RETURN 00202 * 00203 IF( LSAME( DIRECT, 'F' ) ) THEN 00204 PREVLASTV = N 00205 DO I = 1, K 00206 PREVLASTV = MAX( PREVLASTV, I ) 00207 IF( TAU( I ).EQ.ZERO ) THEN 00208 * 00209 * H(i) = I 00210 * 00211 DO J = 1, I 00212 T( J, I ) = ZERO 00213 END DO 00214 ELSE 00215 * 00216 * general case 00217 * 00218 IF( LSAME( STOREV, 'C' ) ) THEN 00219 * Skip any trailing zeros. 00220 DO LASTV = N, I+1, -1 00221 IF( V( LASTV, I ).NE.ZERO ) EXIT 00222 END DO 00223 DO J = 1, I-1 00224 T( J, I ) = -TAU( I ) * CONJG( V( I , J ) ) 00225 END DO 00226 J = MIN( LASTV, PREVLASTV ) 00227 * 00228 * T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i) 00229 * 00230 CALL ZGEMV( 'Conjugate transpose', J-I, I-1, 00231 $ -TAU( I ), V( I+1, 1 ), LDV, 00232 $ V( I+1, I ), 1, ONE, T( 1, I ), 1 ) 00233 ELSE 00234 * Skip any trailing zeros. 00235 DO LASTV = N, I+1, -1 00236 IF( V( I, LASTV ).NE.ZERO ) EXIT 00237 END DO 00238 DO J = 1, I-1 00239 T( J, I ) = -TAU( I ) * V( J , I ) 00240 END DO 00241 J = MIN( LASTV, PREVLASTV ) 00242 * 00243 * T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H 00244 * 00245 CALL ZGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ), 00246 $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, 00247 $ ONE, T( 1, I ), LDT ) 00248 END IF 00249 * 00250 * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) 00251 * 00252 CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, 00253 $ LDT, T( 1, I ), 1 ) 00254 T( I, I ) = TAU( I ) 00255 IF( I.GT.1 ) THEN 00256 PREVLASTV = MAX( PREVLASTV, LASTV ) 00257 ELSE 00258 PREVLASTV = LASTV 00259 END IF 00260 END IF 00261 END DO 00262 ELSE 00263 PREVLASTV = 1 00264 DO I = K, 1, -1 00265 IF( TAU( I ).EQ.ZERO ) THEN 00266 * 00267 * H(i) = I 00268 * 00269 DO J = I, K 00270 T( J, I ) = ZERO 00271 END DO 00272 ELSE 00273 * 00274 * general case 00275 * 00276 IF( I.LT.K ) THEN 00277 IF( LSAME( STOREV, 'C' ) ) THEN 00278 * Skip any leading zeros. 00279 DO LASTV = 1, I-1 00280 IF( V( LASTV, I ).NE.ZERO ) EXIT 00281 END DO 00282 DO J = I+1, K 00283 T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) ) 00284 END DO 00285 J = MAX( LASTV, PREVLASTV ) 00286 * 00287 * T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i) 00288 * 00289 CALL ZGEMV( 'Conjugate transpose', N-K+I-J, K-I, 00290 $ -TAU( I ), V( J, I+1 ), LDV, V( J, I ), 00291 $ 1, ONE, T( I+1, I ), 1 ) 00292 ELSE 00293 * Skip any leading zeros. 00294 DO LASTV = 1, I-1 00295 IF( V( I, LASTV ).NE.ZERO ) EXIT 00296 END DO 00297 DO J = I+1, K 00298 T( J, I ) = -TAU( I ) * V( J, N-K+I ) 00299 END DO 00300 J = MAX( LASTV, PREVLASTV ) 00301 * 00302 * T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H 00303 * 00304 CALL ZGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ), 00305 $ V( I+1, J ), LDV, V( I, J ), LDV, 00306 $ ONE, T( I+1, I ), LDT ) 00307 END IF 00308 * 00309 * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) 00310 * 00311 CALL ZTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, 00312 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) 00313 IF( I.GT.1 ) THEN 00314 PREVLASTV = MIN( PREVLASTV, LASTV ) 00315 ELSE 00316 PREVLASTV = LASTV 00317 END IF 00318 END IF 00319 T( I, I ) = TAU( I ) 00320 END IF 00321 END DO 00322 END IF 00323 RETURN 00324 * 00325 * End of ZLARFT 00326 * 00327 END