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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlarft.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack