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