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