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