source: pacpussensors/trunk/Vislab/lib3dv/eigen/lapack/clarft.f@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

Doc
File size: 10.2 KB
Line 
1*> \brief \b CLARFT
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CLARFT + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarft.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarft.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarft.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
22*
23* .. Scalar Arguments ..
24* CHARACTER DIRECT, STOREV
25* INTEGER K, LDT, LDV, N
26* ..
27* .. Array Arguments ..
28* COMPLEX T( LDT, * ), TAU( * ), V( LDV, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> CLARFT forms the triangular factor T of a complex block reflector H
38*> of order n, which is defined as a product of k elementary reflectors.
39*>
40*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
41*>
42*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
43*>
44*> If STOREV = 'C', the vector which defines the elementary reflector
45*> H(i) is stored in the i-th column of the array V, and
46*>
47*> H = I - V * T * V**H
48*>
49*> If STOREV = 'R', the vector which defines the elementary reflector
50*> H(i) is stored in the i-th row of the array V, and
51*>
52*> H = I - V**H * T * V
53*> \endverbatim
54*
55* Arguments:
56* ==========
57*
58*> \param[in] DIRECT
59*> \verbatim
60*> DIRECT is CHARACTER*1
61*> Specifies the order in which the elementary reflectors are
62*> multiplied to form the block reflector:
63*> = 'F': H = H(1) H(2) . . . H(k) (Forward)
64*> = 'B': H = H(k) . . . H(2) H(1) (Backward)
65*> \endverbatim
66*>
67*> \param[in] STOREV
68*> \verbatim
69*> STOREV is CHARACTER*1
70*> Specifies how the vectors which define the elementary
71*> reflectors are stored (see also Further Details):
72*> = 'C': columnwise
73*> = 'R': rowwise
74*> \endverbatim
75*>
76*> \param[in] N
77*> \verbatim
78*> N is INTEGER
79*> The order of the block reflector H. N >= 0.
80*> \endverbatim
81*>
82*> \param[in] K
83*> \verbatim
84*> K is INTEGER
85*> The order of the triangular factor T (= the number of
86*> elementary reflectors). K >= 1.
87*> \endverbatim
88*>
89*> \param[in] V
90*> \verbatim
91*> V is COMPLEX array, dimension
92*> (LDV,K) if STOREV = 'C'
93*> (LDV,N) if STOREV = 'R'
94*> The matrix V. See further details.
95*> \endverbatim
96*>
97*> \param[in] LDV
98*> \verbatim
99*> LDV is INTEGER
100*> The leading dimension of the array V.
101*> If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
102*> \endverbatim
103*>
104*> \param[in] TAU
105*> \verbatim
106*> TAU is COMPLEX array, dimension (K)
107*> TAU(i) must contain the scalar factor of the elementary
108*> reflector H(i).
109*> \endverbatim
110*>
111*> \param[out] T
112*> \verbatim
113*> T is COMPLEX array, dimension (LDT,K)
114*> The k by k triangular factor T of the block reflector.
115*> If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
116*> lower triangular. The rest of the array is not used.
117*> \endverbatim
118*>
119*> \param[in] LDT
120*> \verbatim
121*> LDT is INTEGER
122*> The leading dimension of the array T. LDT >= K.
123*> \endverbatim
124*
125* Authors:
126* ========
127*
128*> \author Univ. of Tennessee
129*> \author Univ. of California Berkeley
130*> \author Univ. of Colorado Denver
131*> \author NAG Ltd.
132*
133*> \date April 2012
134*
135*> \ingroup complexOTHERauxiliary
136*
137*> \par Further Details:
138* =====================
139*>
140*> \verbatim
141*>
142*> The shape of the matrix V and the storage of the vectors which define
143*> the H(i) is best illustrated by the following example with n = 5 and
144*> k = 3. The elements equal to 1 are not stored.
145*>
146*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R':
147*>
148*> V = ( 1 ) V = ( 1 v1 v1 v1 v1 )
149*> ( v1 1 ) ( 1 v2 v2 v2 )
150*> ( v1 v2 1 ) ( 1 v3 v3 )
151*> ( v1 v2 v3 )
152*> ( v1 v2 v3 )
153*>
154*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R':
155*>
156*> V = ( v1 v2 v3 ) V = ( v1 v1 1 )
157*> ( v1 v2 v3 ) ( v2 v2 v2 1 )
158*> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 )
159*> ( 1 v3 )
160*> ( 1 )
161*> \endverbatim
162*>
163* =====================================================================
164 SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
165*
166* -- LAPACK auxiliary routine (version 3.4.1) --
167* -- LAPACK is a software package provided by Univ. of Tennessee, --
168* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169* April 2012
170*
171* .. Scalar Arguments ..
172 CHARACTER DIRECT, STOREV
173 INTEGER K, LDT, LDV, N
174* ..
175* .. Array Arguments ..
176 COMPLEX T( LDT, * ), TAU( * ), V( LDV, * )
177* ..
178*
179* =====================================================================
180*
181* .. Parameters ..
182 COMPLEX ONE, ZERO
183 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ),
184 $ ZERO = ( 0.0E+0, 0.0E+0 ) )
185* ..
186* .. Local Scalars ..
187 INTEGER I, J, PREVLASTV, LASTV
188* ..
189* .. External Subroutines ..
190 EXTERNAL CGEMV, CLACGV, CTRMV
191* ..
192* .. External Functions ..
193 LOGICAL LSAME
194 EXTERNAL LSAME
195* ..
196* .. Executable Statements ..
197*
198* Quick return if possible
199*
200 IF( N.EQ.0 )
201 $ RETURN
202*
203 IF( LSAME( DIRECT, 'F' ) ) THEN
204 PREVLASTV = N
205 DO I = 1, K
206 PREVLASTV = MAX( PREVLASTV, I )
207 IF( TAU( I ).EQ.ZERO ) THEN
208*
209* H(i) = I
210*
211 DO J = 1, I
212 T( J, I ) = ZERO
213 END DO
214 ELSE
215*
216* general case
217*
218 IF( LSAME( STOREV, 'C' ) ) THEN
219* Skip any trailing zeros.
220 DO LASTV = N, I+1, -1
221 IF( V( LASTV, I ).NE.ZERO ) EXIT
222 END DO
223 DO J = 1, I-1
224 T( J, I ) = -TAU( I ) * CONJG( V( I , J ) )
225 END DO
226 J = MIN( LASTV, PREVLASTV )
227*
228* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i)
229*
230 CALL CGEMV( 'Conjugate transpose', J-I, I-1,
231 $ -TAU( I ), V( I+1, 1 ), LDV,
232 $ V( I+1, I ), 1,
233 $ ONE, T( 1, I ), 1 )
234 ELSE
235* Skip any trailing zeros.
236 DO LASTV = N, I+1, -1
237 IF( V( I, LASTV ).NE.ZERO ) EXIT
238 END DO
239 DO J = 1, I-1
240 T( J, I ) = -TAU( I ) * V( J , I )
241 END DO
242 J = MIN( LASTV, PREVLASTV )
243*
244* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H
245*
246 CALL CGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ),
247 $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV,
248 $ ONE, T( 1, I ), LDT )
249 END IF
250*
251* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
252*
253 CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
254 $ LDT, T( 1, I ), 1 )
255 T( I, I ) = TAU( I )
256 IF( I.GT.1 ) THEN
257 PREVLASTV = MAX( PREVLASTV, LASTV )
258 ELSE
259 PREVLASTV = LASTV
260 END IF
261 END IF
262 END DO
263 ELSE
264 PREVLASTV = 1
265 DO I = K, 1, -1
266 IF( TAU( I ).EQ.ZERO ) THEN
267*
268* H(i) = I
269*
270 DO J = I, K
271 T( J, I ) = ZERO
272 END DO
273 ELSE
274*
275* general case
276*
277 IF( I.LT.K ) THEN
278 IF( LSAME( STOREV, 'C' ) ) THEN
279* Skip any leading zeros.
280 DO LASTV = 1, I-1
281 IF( V( LASTV, I ).NE.ZERO ) EXIT
282 END DO
283 DO J = I+1, K
284 T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) )
285 END DO
286 J = MAX( LASTV, PREVLASTV )
287*
288* T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i)
289*
290 CALL CGEMV( 'Conjugate transpose', N-K+I-J, K-I,
291 $ -TAU( I ), V( J, I+1 ), LDV, V( J, I ),
292 $ 1, ONE, T( I+1, I ), 1 )
293 ELSE
294* Skip any leading zeros.
295 DO LASTV = 1, I-1
296 IF( V( I, LASTV ).NE.ZERO ) EXIT
297 END DO
298 DO J = I+1, K
299 T( J, I ) = -TAU( I ) * V( J, N-K+I )
300 END DO
301 J = MAX( LASTV, PREVLASTV )
302*
303* T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H
304*
305 CALL CGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ),
306 $ V( I+1, J ), LDV, V( I, J ), LDV,
307 $ ONE, T( I+1, I ), LDT )
308 END IF
309*
310* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
311*
312 CALL CTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
313 $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
314 IF( I.GT.1 ) THEN
315 PREVLASTV = MIN( PREVLASTV, LASTV )
316 ELSE
317 PREVLASTV = LASTV
318 END IF
319 END IF
320 T( I, I ) = TAU( I )
321 END IF
322 END DO
323 END IF
324 RETURN
325*
326* End of CLARFT
327*
328 END
Note: See TracBrowser for help on using the repository browser.