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