 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) 00002 IMPLICIT NONE 00003 * 00004 * -- LAPACK auxiliary routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER DIRECT, STOREV 00011 INTEGER K, LDT, LDV, N 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DLARFT forms the triangular factor T of a real block reflector H 00021 * of order n, which is defined as a product of k elementary reflectors. 00022 * 00023 * If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; 00024 * 00025 * If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. 00026 * 00027 * If STOREV = 'C', the vector which defines the elementary reflector 00028 * H(i) is stored in the i-th column of the array V, and 00029 * 00030 * H = I - V * T * V**T 00031 * 00032 * If STOREV = 'R', the vector which defines the elementary reflector 00033 * H(i) is stored in the i-th row of the array V, and 00034 * 00035 * H = I - V**T * T * V 00036 * 00037 * Arguments 00038 * ========= 00039 * 00040 * DIRECT (input) CHARACTER*1 00041 * Specifies the order in which the elementary reflectors are 00042 * multiplied to form the block reflector: 00043 * = 'F': H = H(1) H(2) . . . H(k) (Forward) 00044 * = 'B': H = H(k) . . . H(2) H(1) (Backward) 00045 * 00046 * STOREV (input) CHARACTER*1 00047 * Specifies how the vectors which define the elementary 00048 * reflectors are stored (see also Further Details): 00049 * = 'C': columnwise 00050 * = 'R': rowwise 00051 * 00052 * N (input) INTEGER 00053 * The order of the block reflector H. N >= 0. 00054 * 00055 * K (input) INTEGER 00056 * The order of the triangular factor T (= the number of 00057 * elementary reflectors). K >= 1. 00058 * 00059 * V (input/output) DOUBLE PRECISION array, dimension 00060 * (LDV,K) if STOREV = 'C' 00061 * (LDV,N) if STOREV = 'R' 00062 * The matrix V. See further details. 00063 * 00064 * LDV (input) INTEGER 00065 * The leading dimension of the array V. 00066 * If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. 00067 * 00068 * TAU (input) DOUBLE PRECISION array, dimension (K) 00069 * TAU(i) must contain the scalar factor of the elementary 00070 * reflector H(i). 00071 * 00072 * T (output) DOUBLE PRECISION array, dimension (LDT,K) 00073 * The k by k triangular factor T of the block reflector. 00074 * If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is 00075 * lower triangular. The rest of the array is not used. 00076 * 00077 * LDT (input) INTEGER 00078 * The leading dimension of the array T. LDT >= K. 00079 * 00080 * Further Details 00081 * =============== 00082 * 00083 * The shape of the matrix V and the storage of the vectors which define 00084 * the H(i) is best illustrated by the following example with n = 5 and 00085 * k = 3. The elements equal to 1 are not stored; the corresponding 00086 * array elements are modified but restored on exit. The rest of the 00087 * array is not used. 00088 * 00089 * DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': 00090 * 00091 * V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) 00092 * ( v1 1 ) ( 1 v2 v2 v2 ) 00093 * ( v1 v2 1 ) ( 1 v3 v3 ) 00094 * ( v1 v2 v3 ) 00095 * ( v1 v2 v3 ) 00096 * 00097 * DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': 00098 * 00099 * V = ( v1 v2 v3 ) V = ( v1 v1 1 ) 00100 * ( v1 v2 v3 ) ( v2 v2 v2 1 ) 00101 * ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) 00102 * ( 1 v3 ) 00103 * ( 1 ) 00104 * 00105 * ===================================================================== 00106 * 00107 * .. Parameters .. 00108 DOUBLE PRECISION ONE, ZERO 00109 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00110 * .. 00111 * .. Local Scalars .. 00112 INTEGER I, J, PREVLASTV, LASTV 00113 DOUBLE PRECISION VII 00114 * .. 00115 * .. External Subroutines .. 00116 EXTERNAL DGEMV, DTRMV 00117 * .. 00118 * .. External Functions .. 00119 LOGICAL LSAME 00120 EXTERNAL LSAME 00121 * .. 00122 * .. Executable Statements .. 00123 * 00124 * Quick return if possible 00125 * 00126 IF( N.EQ.0 ) 00127 $ RETURN 00128 * 00129 IF( LSAME( DIRECT, 'F' ) ) THEN 00130 PREVLASTV = N 00131 DO 20 I = 1, K 00132 PREVLASTV = MAX( I, PREVLASTV ) 00133 IF( TAU( I ).EQ.ZERO ) THEN 00134 * 00135 * H(i) = I 00136 * 00137 DO 10 J = 1, I 00138 T( J, I ) = ZERO 00139 10 CONTINUE 00140 ELSE 00141 * 00142 * general case 00143 * 00144 VII = V( I, I ) 00145 V( I, I ) = ONE 00146 IF( LSAME( STOREV, 'C' ) ) THEN 00147 ! Skip any trailing zeros. 00148 DO LASTV = N, I+1, -1 00149 IF( V( LASTV, I ).NE.ZERO ) EXIT 00150 END DO 00151 J = MIN( LASTV, PREVLASTV ) 00152 * 00153 * T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i) 00154 * 00155 CALL DGEMV( 'Transpose', J-I+1, I-1, -TAU( I ), 00156 $ V( I, 1 ), LDV, V( I, I ), 1, ZERO, 00157 $ T( 1, I ), 1 ) 00158 ELSE 00159 ! Skip any trailing zeros. 00160 DO LASTV = N, I+1, -1 00161 IF( V( I, LASTV ).NE.ZERO ) EXIT 00162 END DO 00163 J = MIN( LASTV, PREVLASTV ) 00164 * 00165 * T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T 00166 * 00167 CALL DGEMV( 'No transpose', I-1, J-I+1, -TAU( I ), 00168 $ V( 1, I ), LDV, V( I, I ), LDV, ZERO, 00169 $ T( 1, I ), 1 ) 00170 END IF 00171 V( I, I ) = VII 00172 * 00173 * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) 00174 * 00175 CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, 00176 $ LDT, T( 1, I ), 1 ) 00177 T( I, I ) = TAU( I ) 00178 IF( I.GT.1 ) THEN 00179 PREVLASTV = MAX( PREVLASTV, LASTV ) 00180 ELSE 00181 PREVLASTV = LASTV 00182 END IF 00183 END IF 00184 20 CONTINUE 00185 ELSE 00186 PREVLASTV = 1 00187 DO 40 I = K, 1, -1 00188 IF( TAU( I ).EQ.ZERO ) THEN 00189 * 00190 * H(i) = I 00191 * 00192 DO 30 J = I, K 00193 T( J, I ) = ZERO 00194 30 CONTINUE 00195 ELSE 00196 * 00197 * general case 00198 * 00199 IF( I.LT.K ) THEN 00200 IF( LSAME( STOREV, 'C' ) ) THEN 00201 VII = V( N-K+I, I ) 00202 V( N-K+I, I ) = ONE 00203 ! Skip any leading zeros. 00204 DO LASTV = 1, I-1 00205 IF( V( LASTV, I ).NE.ZERO ) EXIT 00206 END DO 00207 J = MAX( LASTV, PREVLASTV ) 00208 * 00209 * T(i+1:k,i) := 00210 * - tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i) 00211 * 00212 CALL DGEMV( 'Transpose', N-K+I-J+1, K-I, -TAU( I ), 00213 $ V( J, I+1 ), LDV, V( J, I ), 1, ZERO, 00214 $ T( I+1, I ), 1 ) 00215 V( N-K+I, I ) = VII 00216 ELSE 00217 VII = V( I, N-K+I ) 00218 V( I, N-K+I ) = ONE 00219 ! Skip any leading zeros. 00220 DO LASTV = 1, I-1 00221 IF( V( I, LASTV ).NE.ZERO ) EXIT 00222 END DO 00223 J = MAX( LASTV, PREVLASTV ) 00224 * 00225 * T(i+1:k,i) := 00226 * - tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T 00227 * 00228 CALL DGEMV( 'No transpose', K-I, N-K+I-J+1, 00229 $ -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV, 00230 $ ZERO, T( I+1, I ), 1 ) 00231 V( I, N-K+I ) = VII 00232 END IF 00233 * 00234 * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) 00235 * 00236 CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, 00237 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) 00238 IF( I.GT.1 ) THEN 00239 PREVLASTV = MIN( PREVLASTV, LASTV ) 00240 ELSE 00241 PREVLASTV = LASTV 00242 END IF 00243 END IF 00244 T( I, I ) = TAU( I ) 00245 END IF 00246 40 CONTINUE 00247 END IF 00248 RETURN 00249 * 00250 * End of DLARFT 00251 * 00252 END
1.7.3 
