![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLARFT 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CLARFT + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarft.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarft.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarft.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CLARFT( 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 T( LDT, * ), TAU( * ), V( LDV, * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> CLARFT 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 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 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 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 complexOTHERauxiliary 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 CLARFT( 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 T( LDT, * ), TAU( * ), V( LDV, * ) 00177 * .. 00178 * 00179 * ===================================================================== 00180 * 00181 * .. Parameters .. 00182 COMPLEX ONE, ZERO 00183 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), 00184 $ ZERO = ( 0.0E+0, 0.0E+0 ) ) 00185 * .. 00186 * .. Local Scalars .. 00187 INTEGER I, J, PREVLASTV, LASTV 00188 * .. 00189 * .. External Subroutines .. 00190 EXTERNAL CGEMV, CLACGV, CTRMV 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 CGEMV( 'Conjugate transpose', J-I, I-1, 00231 $ -TAU( I ), V( I+1, 1 ), LDV, 00232 $ V( I+1, I ), 1, 00233 $ ONE, T( 1, I ), 1 ) 00234 ELSE 00235 * Skip any trailing zeros. 00236 DO LASTV = N, I+1, -1 00237 IF( V( I, LASTV ).NE.ZERO ) EXIT 00238 END DO 00239 DO J = 1, I-1 00240 T( J, I ) = -TAU( I ) * V( J , I ) 00241 END DO 00242 J = MIN( LASTV, PREVLASTV ) 00243 * 00244 * T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H 00245 * 00246 CALL CGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ), 00247 $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, 00248 $ ONE, T( 1, I ), LDT ) 00249 END IF 00250 * 00251 * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) 00252 * 00253 CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, 00254 $ LDT, T( 1, I ), 1 ) 00255 T( I, I ) = TAU( I ) 00256 IF( I.GT.1 ) THEN 00257 PREVLASTV = MAX( PREVLASTV, LASTV ) 00258 ELSE 00259 PREVLASTV = LASTV 00260 END IF 00261 END IF 00262 END DO 00263 ELSE 00264 PREVLASTV = 1 00265 DO I = K, 1, -1 00266 IF( TAU( I ).EQ.ZERO ) THEN 00267 * 00268 * H(i) = I 00269 * 00270 DO J = I, K 00271 T( J, I ) = ZERO 00272 END DO 00273 ELSE 00274 * 00275 * general case 00276 * 00277 IF( I.LT.K ) THEN 00278 IF( LSAME( STOREV, 'C' ) ) THEN 00279 * Skip any leading zeros. 00280 DO LASTV = 1, I-1 00281 IF( V( LASTV, I ).NE.ZERO ) EXIT 00282 END DO 00283 DO J = I+1, K 00284 T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) ) 00285 END DO 00286 J = MAX( LASTV, PREVLASTV ) 00287 * 00288 * T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i) 00289 * 00290 CALL CGEMV( 'Conjugate transpose', N-K+I-J, K-I, 00291 $ -TAU( I ), V( J, I+1 ), LDV, V( J, I ), 00292 $ 1, ONE, T( I+1, I ), 1 ) 00293 ELSE 00294 * Skip any leading zeros. 00295 DO LASTV = 1, I-1 00296 IF( V( I, LASTV ).NE.ZERO ) EXIT 00297 END DO 00298 DO J = I+1, K 00299 T( J, I ) = -TAU( I ) * V( J, N-K+I ) 00300 END DO 00301 J = MAX( LASTV, PREVLASTV ) 00302 * 00303 * T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H 00304 * 00305 CALL CGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ), 00306 $ V( I+1, J ), LDV, V( I, J ), LDV, 00307 $ ONE, T( I+1, I ), LDT ) 00308 END IF 00309 * 00310 * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) 00311 * 00312 CALL CTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, 00313 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) 00314 IF( I.GT.1 ) THEN 00315 PREVLASTV = MIN( PREVLASTV, LASTV ) 00316 ELSE 00317 PREVLASTV = LASTV 00318 END IF 00319 END IF 00320 T( I, I ) = TAU( I ) 00321 END IF 00322 END DO 00323 END IF 00324 RETURN 00325 * 00326 * End of CLARFT 00327 * 00328 END