KTH framework for Nek5000 toolboxes; testing version  0.0.1
dlarft.f
Go to the documentation of this file.
1  SUBROUTINE dlarft( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
2 *
3 * -- LAPACK auxiliary routine (version 3.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * February 29, 1992
7 *
8 * .. Scalar Arguments ..
9  CHARACTER DIRECT, STOREV
10  INTEGER K, LDT, LDV, N
11 * ..
12 * .. Array Arguments ..
13  DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DLARFT forms the triangular factor T of a real 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'
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' * 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (K)
68 * TAU(i) must contain the scalar factor of the elementary
69 * reflector H(i).
70 *
71 * T (output) DOUBLE PRECISION 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  DOUBLE PRECISION ONE, ZERO
108  parameter( one = 1.0d+0, zero = 0.0d+0 )
109 * ..
110 * .. Local Scalars ..
111  INTEGER I, J
112  DOUBLE PRECISION VII
113 * ..
114 * .. External Subroutines ..
115  EXTERNAL dgemv, dtrmv
116 * ..
117 * .. External Functions ..
118  LOGICAL LSAME
119  EXTERNAL lsame
120 * ..
121 * .. Executable Statements ..
122 *
123 * Quick return if possible
124 *
125  IF( n.EQ.0 )
126  $ RETURN
127 *
128  IF( lsame( direct, 'F' ) ) THEN
129  DO 20 i = 1, k
130  IF( tau( i ).EQ.zero ) THEN
131 *
132 * H(i) = I
133 *
134  DO 10 j = 1, i
135  t( j, i ) = zero
136  10 CONTINUE
137  ELSE
138 *
139 * general case
140 *
141  vii = v( i, i )
142  v( i, i ) = one
143  IF( lsame( storev, 'C' ) ) THEN
144 *
145 * T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i)
146 *
147  CALL dgemv( 'Transpose', n-i+1, i-1, -tau( i ),
148  $ v( i, 1 ), ldv, v( i, i ), 1, zero,
149  $ t( 1, i ), 1 )
150  ELSE
151 *
152 * T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)'
153 *
154  CALL dgemv( 'No transpose', i-1, n-i+1, -tau( i ),
155  $ v( 1, i ), ldv, v( i, i ), ldv, zero,
156  $ t( 1, i ), 1 )
157  END IF
158  v( i, i ) = vii
159 *
160 * T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
161 *
162  CALL dtrmv( 'Upper', 'No transpose', 'Non-unit', i-1, t,
163  $ ldt, t( 1, i ), 1 )
164  t( i, i ) = tau( i )
165  END IF
166  20 CONTINUE
167  ELSE
168  DO 40 i = k, 1, -1
169  IF( tau( i ).EQ.zero ) THEN
170 *
171 * H(i) = I
172 *
173  DO 30 j = i, k
174  t( j, i ) = zero
175  30 CONTINUE
176  ELSE
177 *
178 * general case
179 *
180  IF( i.LT.k ) THEN
181  IF( lsame( storev, 'C' ) ) THEN
182  vii = v( n-k+i, i )
183  v( n-k+i, i ) = one
184 *
185 * T(i+1:k,i) :=
186 * - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i)
187 *
188  CALL dgemv( 'Transpose', n-k+i, k-i, -tau( i ),
189  $ v( 1, i+1 ), ldv, v( 1, i ), 1, zero,
190  $ t( i+1, i ), 1 )
191  v( n-k+i, i ) = vii
192  ELSE
193  vii = v( i, n-k+i )
194  v( i, n-k+i ) = one
195 *
196 * T(i+1:k,i) :=
197 * - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)'
198 *
199  CALL dgemv( 'No transpose', k-i, n-k+i, -tau( i ),
200  $ v( i+1, 1 ), ldv, v( i, 1 ), ldv, zero,
201  $ t( i+1, i ), 1 )
202  v( i, n-k+i ) = vii
203  END IF
204 *
205 * T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
206 *
207  CALL dtrmv( 'Lower', 'No transpose', 'Non-unit', k-i,
208  $ t( i+1, i+1 ), ldt, t( i+1, i ), 1 )
209  END IF
210  t( i, i ) = tau( i )
211  END IF
212  40 CONTINUE
213  END IF
214  RETURN
215 *
216 * End of DLARFT
217 *
218  END
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
Definition: dgemv.f:3
subroutine dlarft(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
Definition: dlarft.f:2
subroutine dtrmv(UPLO, TRANS, DIAG, N, A, LDA, X, INCX)
Definition: dtrmv.f:2