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