SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlaref.f
Go to the documentation of this file.
1 SUBROUTINE dlaref( TYPE, A, LDA, WANTZ, Z, LDZ, BLOCK, IROW1,
2 $ ICOL1, ISTART, ISTOP, ITMP1, ITMP2, LILOZ,
3 $ LIHIZ, VECS, V2, V3, T1, T2, T3 )
4 IMPLICIT NONE
5*
6* -- ScaLAPACK auxiliary routine (version 1.5) --
7* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
8* and University of California, Berkeley.
9* May 1, 1997
10*
11* .. Scalar Arguments ..
12 LOGICAL BLOCK, WANTZ
13 CHARACTER TYPE
14 INTEGER ICOL1, IROW1, ISTART, ISTOP, ITMP1, ITMP2, LDA,
15 $ ldz, lihiz, liloz
16 DOUBLE PRECISION T1, T2, T3, V2, V3
17* ..
18* .. Array Arguments ..
19 DOUBLE PRECISION A( LDA, * ), VECS( * ), Z( LDZ, * )
20* ..
21*
22* Purpose
23* =======
24*
25* DLAREF applies one or several Householder reflectors of size 3
26* to one or two matrices (if column is specified) on either their
27* rows or columns.
28*
29* Arguments
30* =========
31*
32* TYPE (global input) CHARACTER*1
33* If 'R': Apply reflectors to the rows of the matrix
34* (apply from left)
35* Otherwise: Apply reflectors to the columns of the matrix
36* Unchanged on exit.
37*
38* A (global input/output) DOUBLE PRECISION array, (LDA,*)
39* On entry, the matrix to receive the reflections.
40* The updated matrix on exit.
41*
42* LDA (local input) INTEGER
43* On entry, the leading dimension of A. Unchanged on exit.
44*
45* WANTZ (global input) LOGICAL
46* If .TRUE., then apply any column reflections to Z as well.
47* If .FALSE., then do no additional work on Z.
48*
49* Z (global input/output) DOUBLE PRECISION array, (LDZ,*)
50* On entry, the second matrix to receive column reflections.
51* This is changed only if WANTZ is set.
52*
53* LDZ (local input) INTEGER
54* On entry, the leading dimension of Z. Unchanged on exit.
55*
56* BLOCK (global input) LOGICAL
57* If .TRUE., then apply several reflectors at once and read
58* their data from the VECS array.
59* If .FALSE., apply the single reflector given by V2, V3,
60* T1, T2, and T3.
61*
62* IROW1 (local input/output) INTEGER
63* On entry, the local row element of A.
64* Undefined on output.
65*
66*
67* ICOL1 (local input/output) INTEGER
68* On entry, the local column element of A.
69* Undefined on output.
70*
71* ISTART (global input) INTEGER
72* Specifies the "number" of the first reflector. This is
73* used as an index into VECS if BLOCK is set.
74* ISTART is ignored if BLOCK is .FALSE..
75*
76* ISTOP (global input) INTEGER
77* Specifies the "number" of the last reflector. This is
78* used as an index into VECS if BLOCK is set.
79* ISTOP is ignored if BLOCK is .FALSE..
80*
81* ITMP1 (local input) INTEGER
82* Starting range into A. For rows, this is the local
83* first column. For columns, this is the local first row.
84*
85* ITMP2 (local input) INTEGER
86* Ending range into A. For rows, this is the local last
87* column. For columns, this is the local last row.
88*
89* LILOZ
90* LIHIZ (local input) INTEGER
91* These serve the same purpose as ITMP1,ITMP2 but for Z
92* when WANTZ is set.
93*
94* VECS (global input) DOUBLE PRECISION array of size 3*N (matrix
95* size)
96* This holds the size 3 reflectors one after another and this
97* is only accessed when BLOCK is .TRUE.
98*
99* V2
100* V3
101* T1
102* T2
103* T3 (global input/output) DOUBLE PRECISION
104* This holds information on a single size 3 Householder
105* reflector and is read when BLOCK is .FALSE., and
106* overwritten when BLOCK is .TRUE.
107*
108* Implemented by: G. Henry, May 1, 1997
109*
110* =====================================================================
111*
112* .. Local Scalars ..
113 INTEGER J, K
114 DOUBLE PRECISION H11, H22, SUM, T12, T13, T22, T23, T32, T33,
115 $ v22, v23, v32, v33, a1, a2, a3, a4, a5, b1,
116 $ b2, b3, b4, b5, tmp1, tmp2, tmp3, sum1, sum2,
117 $ sum3, a11, a22
118* ..
119* .. External Functions ..
120 LOGICAL LSAME
121 EXTERNAL LSAME
122* ..
123* .. Intrinsic Functions ..
124 INTRINSIC mod
125* ..
126* .. Executable Statements ..
127*
128 IF( lsame( TYPE, 'R' ) ) then
129 if( block ) THEN
130 DO 30 k = istart, istop - mod( istop-istart+1, 3 ), 3
131 v2 = vecs( ( k-1 )*3+1 )
132 v3 = vecs( ( k-1 )*3+2 )
133 t1 = vecs( ( k-1 )*3+3 )
134 v22 = vecs( ( k-1 )*3+4 )
135 v32 = vecs( ( k-1 )*3+5 )
136 t12 = vecs( ( k-1 )*3+6 )
137 v23 = vecs( ( k-1 )*3+7 )
138 v33 = vecs( ( k-1 )*3+8 )
139 t13 = vecs( ( k-1 )*3+9 )
140 t2 = t1*v2
141 t3 = t1*v3
142 t22 = t12*v22
143 t32 = t12*v32
144 t23 = t13*v23
145 t33 = t13*v33
146 DO 10 j = itmp1, itmp2-mod(itmp2-itmp1+1,2), 2
147 a1 = a( irow1 , j )
148 a2 = a( irow1+1, j )
149 a3 = a( irow1+2, j )
150 a4 = a( irow1+3, j )
151 a5 = a( irow1+4, j )
152 b1 = a( irow1 , j+1 )
153 b2 = a( irow1+1, j+1 )
154 b3 = a( irow1+2, j+1 )
155 b4 = a( irow1+3, j+1 )
156 b5 = a( irow1+4, j+1 )
157 sum1 = a1 + v2*a2 + v3*a3
158 a( irow1 , j ) = a1 - sum1 * t1
159 h11 = a2 - sum1 * t2
160 h22 = a3 - sum1 * t3
161 tmp1 = b1 + v2*b2 + v3*b3
162 a( irow1 , j+1 ) = b1 - tmp1 * t1
163 a11 = b2 - tmp1 * t2
164 a22 = b3 - tmp1 * t3
165 sum2 = h11 + v22*h22 + v32*a4
166 a( irow1+1, j ) = h11 - sum2 * t12
167 h11 = h22 - sum2 * t22
168 h22 = a4 - sum2 * t32
169 tmp2 = a11 + v22*a22 + v32*b4
170 a( irow1+1, j+1 ) = a11 - tmp2 * t12
171 a11 = a22 - tmp2 * t22
172 a22 = b4 - tmp2 * t32
173 sum3 = h11 + v23*h22 + v33*a5
174 a( irow1+2, j ) = h11 - sum3 * t13
175 a( irow1+3, j ) = h22 - sum3 * t23
176 a( irow1+4, j ) = a5 - sum3 * t33
177 tmp3 = a11 + v23*a22 + v33*b5
178 a( irow1+2, j+1 ) = a11 - tmp3 * t13
179 a( irow1+3, j+1 ) = a22 - tmp3 * t23
180 a( irow1+4, j+1 ) = b5 - tmp3 * t33
181 10 CONTINUE
182 DO 20 j = itmp2-mod(itmp2-itmp1+1,2)+1, itmp2
183 sum = a( irow1, j ) + v2*a( irow1+1, j ) +
184 $ v3*a( irow1+2, j )
185 a( irow1, j ) = a( irow1, j ) - sum*t1
186 h11 = a( irow1+1, j ) - sum*t2
187 h22 = a( irow1+2, j ) - sum*t3
188 sum = h11 + v22*h22 + v32*a( irow1+3, j )
189 a( irow1+1, j ) = h11 - sum*t12
190 h11 = h22 - sum*t22
191 h22 = a( irow1+3, j ) - sum*t32
192 sum = h11 + v23*h22 + v33*a( irow1+4, j )
193 a( irow1+2, j ) = h11 - sum*t13
194 a( irow1+3, j ) = h22 - sum*t23
195 a( irow1+4, j ) = a( irow1+4, j ) - sum*t33
196 20 CONTINUE
197 irow1 = irow1 + 3
198 30 CONTINUE
199 DO 50 k = istop - mod( istop-istart+1, 3 ) + 1, istop
200 v2 = vecs( ( k-1 )*3+1 )
201 v3 = vecs( ( k-1 )*3+2 )
202 t1 = vecs( ( k-1 )*3+3 )
203 t2 = t1*v2
204 t3 = t1*v3
205 DO 40 j = itmp1, itmp2
206 sum = a( irow1, j ) + v2*a( irow1+1, j ) +
207 $ v3*a( irow1+2, j )
208 a( irow1, j ) = a( irow1, j ) - sum*t1
209 a( irow1+1, j ) = a( irow1+1, j ) - sum*t2
210 a( irow1+2, j ) = a( irow1+2, j ) - sum*t3
211 40 CONTINUE
212 irow1 = irow1 + 1
213 50 CONTINUE
214 ELSE
215 DO 60 j = itmp1, itmp2
216 sum = a( irow1, j ) + v2*a( irow1+1, j ) +
217 $ v3*a( irow1+2, j )
218 a( irow1, j ) = a( irow1, j ) - sum*t1
219 a( irow1+1, j ) = a( irow1+1, j ) - sum*t2
220 a( irow1+2, j ) = a( irow1+2, j ) - sum*t3
221 60 CONTINUE
222 END IF
223 ELSE
224*
225* Do column transforms
226*
227 IF( block ) THEN
228 DO 90 k = istart, istop - mod( istop-istart+1, 3 ), 3
229 v2 = vecs( ( k-1 )*3+1 )
230 v3 = vecs( ( k-1 )*3+2 )
231 t1 = vecs( ( k-1 )*3+3 )
232 v22 = vecs( ( k-1 )*3+4 )
233 v32 = vecs( ( k-1 )*3+5 )
234 t12 = vecs( ( k-1 )*3+6 )
235 v23 = vecs( ( k-1 )*3+7 )
236 v33 = vecs( ( k-1 )*3+8 )
237 t13 = vecs( ( k-1 )*3+9 )
238 t2 = t1*v2
239 t3 = t1*v3
240 t22 = t12*v22
241 t32 = t12*v32
242 t23 = t13*v23
243 t33 = t13*v33
244 DO 70 j = itmp1, itmp2
245 sum = a( j, icol1 ) + v2*a( j, icol1+1 ) +
246 $ v3*a( j, icol1+2 )
247 a( j, icol1 ) = a( j, icol1 ) - sum*t1
248 h11 = a( j, icol1+1 ) - sum*t2
249 h22 = a( j, icol1+2 ) - sum*t3
250 sum = h11 + v22*h22 + v32*a( j, icol1+3 )
251 a( j, icol1+1 ) = h11 - sum*t12
252 h11 = h22 - sum*t22
253 h22 = a( j, icol1+3 ) - sum*t32
254 sum = h11 + v23*h22 + v33*a( j, icol1+4 )
255 a( j, icol1+2 ) = h11 - sum*t13
256 a( j, icol1+3 ) = h22 - sum*t23
257 a( j, icol1+4 ) = a( j, icol1+4 ) - sum*t33
258 70 CONTINUE
259 IF( wantz ) THEN
260 DO 80 j = liloz, lihiz
261 sum = z( j, icol1 ) + v2*z( j, icol1+1 ) +
262 $ v3*z( j, icol1+2 )
263 z( j, icol1 ) = z( j, icol1 ) - sum*t1
264 h11 = z( j, icol1+1 ) - sum*t2
265 h22 = z( j, icol1+2 ) - sum*t3
266 sum = h11 + v22*h22 + v32*z( j, icol1+3 )
267 z( j, icol1+1 ) = h11 - sum*t12
268 h11 = h22 - sum*t22
269 h22 = z( j, icol1+3 ) - sum*t32
270 sum = h11 + v23*h22 + v33*z( j, icol1+4 )
271 z( j, icol1+2 ) = h11 - sum*t13
272 z( j, icol1+3 ) = h22 - sum*t23
273 z( j, icol1+4 ) = z( j, icol1+4 ) - sum*t33
274 80 CONTINUE
275 END IF
276 icol1 = icol1 + 3
277 90 CONTINUE
278 DO 120 k = istop - mod( istop-istart+1, 3 ) + 1, istop
279 v2 = vecs( ( k-1 )*3+1 )
280 v3 = vecs( ( k-1 )*3+2 )
281 t1 = vecs( ( k-1 )*3+3 )
282 t2 = t1*v2
283 t3 = t1*v3
284 DO 100 j = itmp1, itmp2
285 sum = a( j, icol1 ) + v2*a( j, icol1+1 ) +
286 $ v3*a( j, icol1+2 )
287 a( j, icol1 ) = a( j, icol1 ) - sum*t1
288 a( j, icol1+1 ) = a( j, icol1+1 ) - sum*t2
289 a( j, icol1+2 ) = a( j, icol1+2 ) - sum*t3
290 100 CONTINUE
291 IF( wantz ) THEN
292 DO 110 j = liloz, lihiz
293 sum = z( j, icol1 ) + v2*z( j, icol1+1 ) +
294 $ v3*z( j, icol1+2 )
295 z( j, icol1 ) = z( j, icol1 ) - sum*t1
296 z( j, icol1+1 ) = z( j, icol1+1 ) - sum*t2
297 z( j, icol1+2 ) = z( j, icol1+2 ) - sum*t3
298 110 CONTINUE
299 END IF
300 icol1 = icol1 + 1
301 120 CONTINUE
302 ELSE
303 DO 130 j = itmp1, itmp2
304 sum = a( j, icol1 ) + v2*a( j, icol1+1 ) +
305 $ v3*a( j, icol1+2 )
306 a( j, icol1 ) = a( j, icol1 ) - sum*t1
307 a( j, icol1+1 ) = a( j, icol1+1 ) - sum*t2
308 a( j, icol1+2 ) = a( j, icol1+2 ) - sum*t3
309 130 CONTINUE
310 END IF
311 END IF
312 RETURN
313*
314* End of DLAREF
315*
316 END
dlaref
subroutine dlaref(type, a, lda, wantz, z, ldz, block, irow1, icol1, istart, istop, itmp1, itmp2, liloz, lihiz, vecs, v2, v3, t1, t2, t3)
Definition dlaref.f:4