ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
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