|
2329
|
1 SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) |
|
|
2 * |
|
7034
|
3 * -- LAPACK auxiliary routine (version 3.1) -- |
|
|
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
5 * November 2006 |
|
2329
|
6 * |
|
|
7 * .. Scalar Arguments .. |
|
|
8 CHARACTER DIRECT, STOREV |
|
|
9 INTEGER K, LDT, LDV, N |
|
|
10 * .. |
|
|
11 * .. Array Arguments .. |
|
|
12 DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) |
|
|
13 * .. |
|
|
14 * |
|
|
15 * Purpose |
|
|
16 * ======= |
|
|
17 * |
|
|
18 * DLARFT forms the triangular factor T of a real block reflector H |
|
|
19 * of order n, which is defined as a product of k elementary reflectors. |
|
|
20 * |
|
|
21 * If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; |
|
|
22 * |
|
|
23 * If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. |
|
|
24 * |
|
|
25 * If STOREV = 'C', the vector which defines the elementary reflector |
|
|
26 * H(i) is stored in the i-th column of the array V, and |
|
|
27 * |
|
|
28 * H = I - V * T * V' |
|
|
29 * |
|
|
30 * If STOREV = 'R', the vector which defines the elementary reflector |
|
|
31 * H(i) is stored in the i-th row of the array V, and |
|
|
32 * |
|
|
33 * H = I - V' * T * V |
|
|
34 * |
|
|
35 * Arguments |
|
|
36 * ========= |
|
|
37 * |
|
|
38 * DIRECT (input) CHARACTER*1 |
|
|
39 * Specifies the order in which the elementary reflectors are |
|
|
40 * multiplied to form the block reflector: |
|
|
41 * = 'F': H = H(1) H(2) . . . H(k) (Forward) |
|
|
42 * = 'B': H = H(k) . . . H(2) H(1) (Backward) |
|
|
43 * |
|
|
44 * STOREV (input) CHARACTER*1 |
|
|
45 * Specifies how the vectors which define the elementary |
|
|
46 * reflectors are stored (see also Further Details): |
|
|
47 * = 'C': columnwise |
|
|
48 * = 'R': rowwise |
|
|
49 * |
|
|
50 * N (input) INTEGER |
|
|
51 * The order of the block reflector H. N >= 0. |
|
|
52 * |
|
|
53 * K (input) INTEGER |
|
|
54 * The order of the triangular factor T (= the number of |
|
|
55 * elementary reflectors). K >= 1. |
|
|
56 * |
|
|
57 * V (input/output) DOUBLE PRECISION array, dimension |
|
|
58 * (LDV,K) if STOREV = 'C' |
|
|
59 * (LDV,N) if STOREV = 'R' |
|
|
60 * The matrix V. See further details. |
|
|
61 * |
|
|
62 * LDV (input) INTEGER |
|
|
63 * The leading dimension of the array V. |
|
|
64 * If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. |
|
|
65 * |
|
|
66 * TAU (input) DOUBLE PRECISION array, dimension (K) |
|
|
67 * TAU(i) must contain the scalar factor of the elementary |
|
|
68 * reflector H(i). |
|
|
69 * |
|
|
70 * T (output) DOUBLE PRECISION array, dimension (LDT,K) |
|
|
71 * The k by k triangular factor T of the block reflector. |
|
|
72 * If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is |
|
|
73 * lower triangular. The rest of the array is not used. |
|
|
74 * |
|
|
75 * LDT (input) INTEGER |
|
|
76 * The leading dimension of the array T. LDT >= K. |
|
|
77 * |
|
|
78 * Further Details |
|
|
79 * =============== |
|
|
80 * |
|
|
81 * The shape of the matrix V and the storage of the vectors which define |
|
|
82 * the H(i) is best illustrated by the following example with n = 5 and |
|
|
83 * k = 3. The elements equal to 1 are not stored; the corresponding |
|
|
84 * array elements are modified but restored on exit. The rest of the |
|
|
85 * array is not used. |
|
|
86 * |
|
|
87 * DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': |
|
|
88 * |
|
|
89 * V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) |
|
|
90 * ( v1 1 ) ( 1 v2 v2 v2 ) |
|
|
91 * ( v1 v2 1 ) ( 1 v3 v3 ) |
|
|
92 * ( v1 v2 v3 ) |
|
|
93 * ( v1 v2 v3 ) |
|
|
94 * |
|
|
95 * DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': |
|
|
96 * |
|
|
97 * V = ( v1 v2 v3 ) V = ( v1 v1 1 ) |
|
|
98 * ( v1 v2 v3 ) ( v2 v2 v2 1 ) |
|
|
99 * ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) |
|
|
100 * ( 1 v3 ) |
|
|
101 * ( 1 ) |
|
|
102 * |
|
|
103 * ===================================================================== |
|
|
104 * |
|
|
105 * .. Parameters .. |
|
|
106 DOUBLE PRECISION ONE, ZERO |
|
|
107 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
|
|
108 * .. |
|
|
109 * .. Local Scalars .. |
|
|
110 INTEGER I, J |
|
|
111 DOUBLE PRECISION VII |
|
|
112 * .. |
|
|
113 * .. External Subroutines .. |
|
|
114 EXTERNAL DGEMV, DTRMV |
|
|
115 * .. |
|
|
116 * .. External Functions .. |
|
|
117 LOGICAL LSAME |
|
|
118 EXTERNAL LSAME |
|
|
119 * .. |
|
|
120 * .. Executable Statements .. |
|
|
121 * |
|
|
122 * Quick return if possible |
|
|
123 * |
|
|
124 IF( N.EQ.0 ) |
|
|
125 $ RETURN |
|
|
126 * |
|
|
127 IF( LSAME( DIRECT, 'F' ) ) THEN |
|
|
128 DO 20 I = 1, K |
|
|
129 IF( TAU( I ).EQ.ZERO ) THEN |
|
|
130 * |
|
|
131 * H(i) = I |
|
|
132 * |
|
|
133 DO 10 J = 1, I |
|
|
134 T( J, I ) = ZERO |
|
|
135 10 CONTINUE |
|
|
136 ELSE |
|
|
137 * |
|
|
138 * general case |
|
|
139 * |
|
|
140 VII = V( I, I ) |
|
|
141 V( I, I ) = ONE |
|
|
142 IF( LSAME( STOREV, 'C' ) ) THEN |
|
|
143 * |
|
|
144 * T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i) |
|
|
145 * |
|
|
146 CALL DGEMV( 'Transpose', N-I+1, I-1, -TAU( I ), |
|
|
147 $ V( I, 1 ), LDV, V( I, I ), 1, ZERO, |
|
|
148 $ T( 1, I ), 1 ) |
|
|
149 ELSE |
|
|
150 * |
|
|
151 * T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)' |
|
|
152 * |
|
|
153 CALL DGEMV( 'No transpose', I-1, N-I+1, -TAU( I ), |
|
|
154 $ V( 1, I ), LDV, V( I, I ), LDV, ZERO, |
|
|
155 $ T( 1, I ), 1 ) |
|
|
156 END IF |
|
|
157 V( I, I ) = VII |
|
|
158 * |
|
|
159 * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) |
|
|
160 * |
|
|
161 CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, |
|
|
162 $ LDT, T( 1, I ), 1 ) |
|
|
163 T( I, I ) = TAU( I ) |
|
|
164 END IF |
|
|
165 20 CONTINUE |
|
|
166 ELSE |
|
|
167 DO 40 I = K, 1, -1 |
|
|
168 IF( TAU( I ).EQ.ZERO ) THEN |
|
|
169 * |
|
|
170 * H(i) = I |
|
|
171 * |
|
|
172 DO 30 J = I, K |
|
|
173 T( J, I ) = ZERO |
|
|
174 30 CONTINUE |
|
|
175 ELSE |
|
|
176 * |
|
|
177 * general case |
|
|
178 * |
|
|
179 IF( I.LT.K ) THEN |
|
|
180 IF( LSAME( STOREV, 'C' ) ) THEN |
|
|
181 VII = V( N-K+I, I ) |
|
|
182 V( N-K+I, I ) = ONE |
|
|
183 * |
|
|
184 * T(i+1:k,i) := |
|
|
185 * - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i) |
|
|
186 * |
|
|
187 CALL DGEMV( 'Transpose', N-K+I, K-I, -TAU( I ), |
|
|
188 $ V( 1, I+1 ), LDV, V( 1, I ), 1, ZERO, |
|
|
189 $ T( I+1, I ), 1 ) |
|
|
190 V( N-K+I, I ) = VII |
|
|
191 ELSE |
|
|
192 VII = V( I, N-K+I ) |
|
|
193 V( I, N-K+I ) = ONE |
|
|
194 * |
|
|
195 * T(i+1:k,i) := |
|
|
196 * - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)' |
|
|
197 * |
|
|
198 CALL DGEMV( 'No transpose', K-I, N-K+I, -TAU( I ), |
|
|
199 $ V( I+1, 1 ), LDV, V( I, 1 ), LDV, ZERO, |
|
|
200 $ T( I+1, I ), 1 ) |
|
|
201 V( I, N-K+I ) = VII |
|
|
202 END IF |
|
|
203 * |
|
|
204 * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) |
|
|
205 * |
|
|
206 CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, |
|
|
207 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) |
|
|
208 END IF |
|
|
209 T( I, I ) = TAU( I ) |
|
|
210 END IF |
|
|
211 40 CONTINUE |
|
|
212 END IF |
|
|
213 RETURN |
|
|
214 * |
|
|
215 * End of DLARFT |
|
|
216 * |
|
|
217 END |
