![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarft.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dlarft.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
2: IMPLICIT NONE
3: *
4: * -- LAPACK auxiliary routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
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'
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 * 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)' * 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)'
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)' * 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)'
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