![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLARFT 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLARFT + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarft.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarft.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarft.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLARFT( 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 * REAL T( LDT, * ), TAU( * ), V( LDV, * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> SLARFT forms the triangular factor T of a real 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**T 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**T * 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 REAL 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 REAL 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 REAL 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 realOTHERauxiliary 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 SLARFT( 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 REAL T( LDT, * ), TAU( * ), V( LDV, * ) 00177 * .. 00178 * 00179 * ===================================================================== 00180 * 00181 * .. Parameters .. 00182 REAL ONE, ZERO 00183 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00184 * .. 00185 * .. Local Scalars .. 00186 INTEGER I, J, PREVLASTV, LASTV 00187 * .. 00188 * .. External Subroutines .. 00189 EXTERNAL SGEMV, STRMV 00190 * .. 00191 * .. External Functions .. 00192 LOGICAL LSAME 00193 EXTERNAL LSAME 00194 * .. 00195 * .. Executable Statements .. 00196 * 00197 * Quick return if possible 00198 * 00199 IF( N.EQ.0 ) 00200 $ RETURN 00201 * 00202 IF( LSAME( DIRECT, 'F' ) ) THEN 00203 PREVLASTV = N 00204 DO I = 1, K 00205 PREVLASTV = MAX( I, PREVLASTV ) 00206 IF( TAU( I ).EQ.ZERO ) THEN 00207 * 00208 * H(i) = I 00209 * 00210 DO J = 1, I 00211 T( J, I ) = ZERO 00212 END DO 00213 ELSE 00214 * 00215 * general case 00216 * 00217 IF( LSAME( STOREV, 'C' ) ) THEN 00218 * Skip any trailing zeros. 00219 DO LASTV = N, I+1, -1 00220 IF( V( LASTV, I ).NE.ZERO ) EXIT 00221 END DO 00222 DO J = 1, I-1 00223 T( J, I ) = -TAU( I ) * V( I , J ) 00224 END DO 00225 J = MIN( LASTV, PREVLASTV ) 00226 * 00227 * T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i) 00228 * 00229 CALL SGEMV( 'Transpose', J-I, I-1, -TAU( I ), 00230 $ V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE, 00231 $ T( 1, I ), 1 ) 00232 ELSE 00233 * Skip any trailing zeros. 00234 DO LASTV = N, I+1, -1 00235 IF( V( I, LASTV ).NE.ZERO ) EXIT 00236 END DO 00237 DO J = 1, I-1 00238 T( J, I ) = -TAU( I ) * V( J , I ) 00239 END DO 00240 J = MIN( LASTV, PREVLASTV ) 00241 * 00242 * T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T 00243 * 00244 CALL SGEMV( 'No transpose', I-1, J-I, -TAU( I ), 00245 $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, 00246 $ ONE, T( 1, I ), 1 ) 00247 END IF 00248 * 00249 * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) 00250 * 00251 CALL STRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, 00252 $ LDT, T( 1, I ), 1 ) 00253 T( I, I ) = TAU( I ) 00254 IF( I.GT.1 ) THEN 00255 PREVLASTV = MAX( PREVLASTV, LASTV ) 00256 ELSE 00257 PREVLASTV = LASTV 00258 END IF 00259 END IF 00260 END DO 00261 ELSE 00262 PREVLASTV = 1 00263 DO I = K, 1, -1 00264 IF( TAU( I ).EQ.ZERO ) THEN 00265 * 00266 * H(i) = I 00267 * 00268 DO J = I, K 00269 T( J, I ) = ZERO 00270 END DO 00271 ELSE 00272 * 00273 * general case 00274 * 00275 IF( I.LT.K ) THEN 00276 IF( LSAME( STOREV, 'C' ) ) THEN 00277 * Skip any leading zeros. 00278 DO LASTV = 1, I-1 00279 IF( V( LASTV, I ).NE.ZERO ) EXIT 00280 END DO 00281 DO J = I+1, K 00282 T( J, I ) = -TAU( I ) * V( N-K+I , J ) 00283 END DO 00284 J = MAX( LASTV, PREVLASTV ) 00285 * 00286 * T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i) 00287 * 00288 CALL SGEMV( 'Transpose', N-K+I-J, K-I, -TAU( I ), 00289 $ V( J, I+1 ), LDV, V( J, I ), 1, ONE, 00290 $ T( I+1, I ), 1 ) 00291 ELSE 00292 * Skip any leading zeros. 00293 DO LASTV = 1, I-1 00294 IF( V( I, LASTV ).NE.ZERO ) EXIT 00295 END DO 00296 DO J = I+1, K 00297 T( J, I ) = -TAU( I ) * V( J, N-K+I ) 00298 END DO 00299 J = MAX( LASTV, PREVLASTV ) 00300 * 00301 * T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T 00302 * 00303 CALL SGEMV( 'No transpose', K-I, N-K+I-J, 00304 $ -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV, 00305 $ ONE, T( I+1, I ), 1 ) 00306 END IF 00307 * 00308 * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) 00309 * 00310 CALL STRMV( 'Lower', 'No transpose', 'Non-unit', K-I, 00311 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) 00312 IF( I.GT.1 ) THEN 00313 PREVLASTV = MIN( PREVLASTV, LASTV ) 00314 ELSE 00315 PREVLASTV = LASTV 00316 END IF 00317 END IF 00318 T( I, I ) = TAU( I ) 00319 END IF 00320 END DO 00321 END IF 00322 RETURN 00323 * 00324 * End of SLARFT 00325 * 00326 END