2329
|
1 SUBROUTINE ZLARFT( 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 COMPLEX*16 T( LDT, * ), TAU( * ), V( LDV, * ) |
|
13 * .. |
|
14 * |
|
15 * Purpose |
|
16 * ======= |
|
17 * |
|
18 * ZLARFT forms the triangular factor T of a complex 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) COMPLEX*16 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) COMPLEX*16 array, dimension (K) |
|
67 * TAU(i) must contain the scalar factor of the elementary |
|
68 * reflector H(i). |
|
69 * |
|
70 * T (output) COMPLEX*16 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 COMPLEX*16 ONE, ZERO |
|
107 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), |
|
108 $ ZERO = ( 0.0D+0, 0.0D+0 ) ) |
|
109 * .. |
|
110 * .. Local Scalars .. |
|
111 INTEGER I, J |
|
112 COMPLEX*16 VII |
|
113 * .. |
|
114 * .. External Subroutines .. |
|
115 EXTERNAL ZGEMV, ZLACGV, ZTRMV |
|
116 * .. |
|
117 * .. External Functions .. |
|
118 LOGICAL LSAME |
|
119 EXTERNAL LSAME |
|
120 * .. |
|
121 * .. Executable Statements .. |
|
122 * |
|
123 * Quick return if possible |
|
124 * |
|
125 IF( N.EQ.0 ) |
|
126 $ RETURN |
|
127 * |
|
128 IF( LSAME( DIRECT, 'F' ) ) THEN |
|
129 DO 20 I = 1, K |
|
130 IF( TAU( I ).EQ.ZERO ) THEN |
|
131 * |
|
132 * H(i) = I |
|
133 * |
|
134 DO 10 J = 1, I |
|
135 T( J, I ) = ZERO |
|
136 10 CONTINUE |
|
137 ELSE |
|
138 * |
|
139 * general case |
|
140 * |
|
141 VII = V( I, I ) |
|
142 V( I, I ) = ONE |
|
143 IF( LSAME( STOREV, 'C' ) ) THEN |
|
144 * |
|
145 * T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i) |
|
146 * |
|
147 CALL ZGEMV( 'Conjugate transpose', N-I+1, I-1, |
|
148 $ -TAU( I ), V( I, 1 ), LDV, V( I, I ), 1, |
|
149 $ ZERO, T( 1, I ), 1 ) |
|
150 ELSE |
|
151 * |
|
152 * T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)' |
|
153 * |
|
154 IF( I.LT.N ) |
|
155 $ CALL ZLACGV( N-I, V( I, I+1 ), LDV ) |
|
156 CALL ZGEMV( 'No transpose', I-1, N-I+1, -TAU( I ), |
|
157 $ V( 1, I ), LDV, V( I, I ), LDV, ZERO, |
|
158 $ T( 1, I ), 1 ) |
|
159 IF( I.LT.N ) |
|
160 $ CALL ZLACGV( N-I, V( I, I+1 ), LDV ) |
|
161 END IF |
|
162 V( I, I ) = VII |
|
163 * |
|
164 * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) |
|
165 * |
|
166 CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, |
|
167 $ LDT, T( 1, I ), 1 ) |
|
168 T( I, I ) = TAU( I ) |
|
169 END IF |
|
170 20 CONTINUE |
|
171 ELSE |
|
172 DO 40 I = K, 1, -1 |
|
173 IF( TAU( I ).EQ.ZERO ) THEN |
|
174 * |
|
175 * H(i) = I |
|
176 * |
|
177 DO 30 J = I, K |
|
178 T( J, I ) = ZERO |
|
179 30 CONTINUE |
|
180 ELSE |
|
181 * |
|
182 * general case |
|
183 * |
|
184 IF( I.LT.K ) THEN |
|
185 IF( LSAME( STOREV, 'C' ) ) THEN |
|
186 VII = V( N-K+I, I ) |
|
187 V( N-K+I, I ) = ONE |
|
188 * |
|
189 * T(i+1:k,i) := |
|
190 * - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i) |
|
191 * |
|
192 CALL ZGEMV( 'Conjugate transpose', N-K+I, K-I, |
|
193 $ -TAU( I ), V( 1, I+1 ), LDV, V( 1, I ), |
|
194 $ 1, ZERO, T( I+1, I ), 1 ) |
|
195 V( N-K+I, I ) = VII |
|
196 ELSE |
|
197 VII = V( I, N-K+I ) |
|
198 V( I, N-K+I ) = ONE |
|
199 * |
|
200 * T(i+1:k,i) := |
|
201 * - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)' |
|
202 * |
|
203 CALL ZLACGV( N-K+I-1, V( I, 1 ), LDV ) |
|
204 CALL ZGEMV( 'No transpose', K-I, N-K+I, -TAU( I ), |
|
205 $ V( I+1, 1 ), LDV, V( I, 1 ), LDV, ZERO, |
|
206 $ T( I+1, I ), 1 ) |
|
207 CALL ZLACGV( N-K+I-1, V( I, 1 ), LDV ) |
|
208 V( I, N-K+I ) = VII |
|
209 END IF |
|
210 * |
|
211 * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) |
|
212 * |
|
213 CALL ZTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, |
|
214 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) |
|
215 END IF |
|
216 T( I, I ) = TAU( I ) |
|
217 END IF |
|
218 40 CONTINUE |
|
219 END IF |
|
220 RETURN |
|
221 * |
|
222 * End of ZLARFT |
|
223 * |
|
224 END |