SRC\dlaref.f
#lines: 278   size: 10 Kb   creation: 18/01/2006 23:36:04   last modification: 08/05/2008 18:37:40   attribute: ARCH    Find   Reload   
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
203:
204:
205:
206:
207:
208:
209:
210:
211:
212:
213:
214:
215:
216:
217:
218:
219:
220:
221:
222:
223:
224:
225:
226:
227:
228:
229:
230:
231:
232:
233:
234:
235:
236:
237:
238:
239:
240:
241:
242:
243:
244:
245:
246:
247:
248:
249:
250:
251:
252:
253:
254:
255:
256:
257:
258:
259:
260:
261:
262:
263:
264:
265:
266:
267:
268:
269:
270:
271:
272:
273:
274:
275:
276:
277:
278:
 
      SUBROUTINE DLAREF( TYPE, A, LDA, WANTZ, Z, LDZ, BLOCK, IROW1,
     $                    ICOL1, ISTART, ISTOP, ITMP1, ITMP2, LILOZ,
     $                    LIHIZ, VECS, V2, V3, T1, T2, T3 )
*
*  -- ScaLAPACK routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     December 31, 1998
*
*     .. Scalar Arguments ..
      LOGICAL            BLOCK, WANTZ
      CHARACTER          TYPE
      INTEGER            ICOL1, IROW1, ISTART, ISTOP, ITMP1, ITMP2, LDA,
     $                   LDZ, LIHIZ, LILOZ
      DOUBLE PRECISION   T1, T2, T3, V2, V3
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), VECS( * ), Z( LDZ, * )
*     ..
*
*  Purpose
*  =======
*
*  DLAREF applies one or several Householder reflectors of size 3
*     to one or two matrices (if column is specified) on either their
*     rows or columns.
*
*  Arguments
*  =========
*
*  TYPE    (global input) CHARACTER*1
*          If 'R': Apply reflectors to the rows of the matrix
*              (apply from left)
*          Otherwise: Apply reflectors to the columns of the matrix
*          Unchanged on exit.
*
*  A       (global input/output) DOUBLE PRECISION array, (LDA,*)
*          On entry, the matrix to receive the reflections.
*          The updated matrix on exit.
*
*  LDA     (local input) INTEGER
*          On entry, the leading dimension of A.  Unchanged on exit.
*
*  WANTZ   (global input) LOGICAL
*          If .TRUE., then apply any column reflections to Z as well.
*          If .FALSE., then do no additional work on Z.
*
*  Z       (global input/output) DOUBLE PRECISION array, (LDZ,*)
*          On entry, the second matrix to receive column reflections.
*          This is changed only if WANTZ is set.
*
*  LDZ     (local input) INTEGER
*          On entry, the leading dimension of Z.  Unchanged on exit.
*
*  BLOCK   (global input) LOGICAL
*          If .TRUE., then apply several reflectors at once and read
*             their data from the VECS array.
*          If .FALSE., apply the single reflector given by V2, V3,
*             T1, T2, and T3.
*
*  IROW1   (local input/output) INTEGER
*          On entry, the local row element of A.
*          Undefined on output.
*
*
*  ICOL1   (local input/output) INTEGER
*          On entry, the local column element of A.
*          Undefined on output.
*
*  ISTART  (global input) INTEGER
*          Specifies the "number" of the first reflector.  This is
*              used as an index into VECS if BLOCK is set.
*              ISTART is ignored if BLOCK is .FALSE..
*
*  ISTOP   (global input) INTEGER
*          Specifies the "number" of the last reflector.  This is
*              used as an index into VECS if BLOCK is set.
*              ISTOP is ignored if BLOCK is .FALSE..
*
*  ITMP1   (local input) INTEGER
*          Starting range into A.  For rows, this is the local
*              first column.  For columns, this is the local first row.
*
*  ITMP2   (local input) INTEGER
*          Ending range into A.  For rows, this is the local last
*              column.  For columns, this is the local last row.
*
*  LILOZ
*  LIHIZ   (local input) INTEGER
*          These serve the same purpose as ITMP1,ITMP2 but for Z
*              when WANTZ is set.
*
*  VECS    (global input) DOUBLE PRECISION array of size 3*N (matrix
*                                                             size)
*          This holds the size 3 reflectors one after another and this
*              is only accessed when BLOCK is .TRUE.
*
*  V2
*  V3
*  T1
*  T2
*  T3      (global input/output) DOUBLE PRECISION
*          This holds information on a single size 3 Householder
*              reflector and is read when BLOCK is .FALSE., and
*              overwritten when BLOCK is .TRUE.
*
*  Implemented by:  G. Henry, November 17, 1996
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            J, K
      DOUBLE PRECISION   H11, H22, SUM, T12, T13, T22, T23, T32, T33,
     $                   V22, V23, V32, V33
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MOD
*     ..
*     .. Executable Statements ..
*
      IF( LSAME( TYPE, 'R' ) ) THEN
         IF( BLOCK ) THEN
            DO 20 K = ISTART, ISTOP - MOD( ISTOP-ISTART+1, 3 ), 3
               V2 = VECS( ( K-1 )*3+1 )
               V3 = VECS( ( K-1 )*3+2 )
               T1 = VECS( ( K-1 )*3+3 )
               V22 = VECS( ( K-1 )*3+4 )
               V32 = VECS( ( K-1 )*3+5 )
               T12 = VECS( ( K-1 )*3+6 )
               V23 = VECS( ( K-1 )*3+7 )
               V33 = VECS( ( K-1 )*3+8 )
               T13 = VECS( ( K-1 )*3+9 )
               T2 = T1*V2
               T3 = T1*V3
               T22 = T12*V22
               T32 = T12*V32
               T23 = T13*V23
               T33 = T13*V33
               DO 10 J = ITMP1, ITMP2
                  SUM = A( IROW1, J ) + V2*A( IROW1+1, J ) +
     $                  V3*A( IROW1+2, J )
                  A( IROW1, J ) = A( IROW1, J ) - SUM*T1
                  H11 = A( IROW1+1, J ) - SUM*T2
                  H22 = A( IROW1+2, J ) - SUM*T3
                  SUM = H11 + V22*H22 + V32*A( IROW1+3, J )
                  A( IROW1+1, J ) = H11 - SUM*T12
                  H11 = H22 - SUM*T22
                  H22 = A( IROW1+3, J ) - SUM*T32
                  SUM = H11 + V23*H22 + V33*A( IROW1+4, J )
                  A( IROW1+2, J ) = H11 - SUM*T13
                  A( IROW1+3, J ) = H22 - SUM*T23
                  A( IROW1+4, J ) = A( IROW1+4, J ) - SUM*T33
   10          CONTINUE
               IROW1 = IROW1 + 3
   20       CONTINUE
            DO 40 K = ISTOP - MOD( ISTOP-ISTART+1, 3 ) + 1, ISTOP
               V2 = VECS( ( K-1 )*3+1 )
               V3 = VECS( ( K-1 )*3+2 )
               T1 = VECS( ( K-1 )*3+3 )
               T2 = T1*V2
               T3 = T1*V3
               DO 30 J = ITMP1, ITMP2
                  SUM = A( IROW1, J ) + V2*A( IROW1+1, J ) +
     $                  V3*A( IROW1+2, J )
                  A( IROW1, J ) = A( IROW1, J ) - SUM*T1
                  A( IROW1+1, J ) = A( IROW1+1, J ) - SUM*T2
                  A( IROW1+2, J ) = A( IROW1+2, J ) - SUM*T3
   30          CONTINUE
               IROW1 = IROW1 + 1
   40       CONTINUE
         ELSE
            DO 50 J = ITMP1, ITMP2
               SUM = A( IROW1, J ) + V2*A( IROW1+1, J ) +
     $               V3*A( IROW1+2, J )
               A( IROW1, J ) = A( IROW1, J ) - SUM*T1
               A( IROW1+1, J ) = A( IROW1+1, J ) - SUM*T2
               A( IROW1+2, J ) = A( IROW1+2, J ) - SUM*T3
   50       CONTINUE
         END IF
      ELSE
*
*        Do column transforms
*
         IF( BLOCK ) THEN
            DO 80 K = ISTART, ISTOP - MOD( ISTOP-ISTART+1, 3 ), 3
               V2 = VECS( ( K-1 )*3+1 )
               V3 = VECS( ( K-1 )*3+2 )
               T1 = VECS( ( K-1 )*3+3 )
               V22 = VECS( ( K-1 )*3+4 )
               V32 = VECS( ( K-1 )*3+5 )
               T12 = VECS( ( K-1 )*3+6 )
               V23 = VECS( ( K-1 )*3+7 )
               V33 = VECS( ( K-1 )*3+8 )
               T13 = VECS( ( K-1 )*3+9 )
               T2 = T1*V2
               T3 = T1*V3
               T22 = T12*V22
               T32 = T12*V32
               T23 = T13*V23
               T33 = T13*V33
               DO 60 J = ITMP1, ITMP2
                  SUM = A( J, ICOL1 ) + V2*A( J, ICOL1+1 ) +
     $                  V3*A( J, ICOL1+2 )
                  A( J, ICOL1 ) = A( J, ICOL1 ) - SUM*T1
                  H11 = A( J, ICOL1+1 ) - SUM*T2
                  H22 = A( J, ICOL1+2 ) - SUM*T3
                  SUM = H11 + V22*H22 + V32*A( J, ICOL1+3 )
                  A( J, ICOL1+1 ) = H11 - SUM*T12
                  H11 = H22 - SUM*T22
                  H22 = A( J, ICOL1+3 ) - SUM*T32
                  SUM = H11 + V23*H22 + V33*A( J, ICOL1+4 )
                  A( J, ICOL1+2 ) = H11 - SUM*T13
                  A( J, ICOL1+3 ) = H22 - SUM*T23
                  A( J, ICOL1+4 ) = A( J, ICOL1+4 ) - SUM*T33
   60          CONTINUE
               IF( WANTZ ) THEN
                  DO 70 J = LILOZ, LIHIZ
                     SUM = Z( J, ICOL1 ) + V2*Z( J, ICOL1+1 ) +
     $                     V3*Z( J, ICOL1+2 )
                     Z( J, ICOL1 ) = Z( J, ICOL1 ) - SUM*T1
                     H11 = Z( J, ICOL1+1 ) - SUM*T2
                     H22 = Z( J, ICOL1+2 ) - SUM*T3
                     SUM = H11 + V22*H22 + V32*Z( J, ICOL1+3 )
                     Z( J, ICOL1+1 ) = H11 - SUM*T12
                     H11 = H22 - SUM*T22
                     H22 = Z( J, ICOL1+3 ) - SUM*T32
                     SUM = H11 + V23*H22 + V33*Z( J, ICOL1+4 )
                     Z( J, ICOL1+2 ) = H11 - SUM*T13
                     Z( J, ICOL1+3 ) = H22 - SUM*T23
                     Z( J, ICOL1+4 ) = Z( J, ICOL1+4 ) - SUM*T33
   70             CONTINUE
               END IF
               ICOL1 = ICOL1 + 3
   80       CONTINUE
            DO 110 K = ISTOP - MOD( ISTOP-ISTART+1, 3 ) + 1, ISTOP
               V2 = VECS( ( K-1 )*3+1 )
               V3 = VECS( ( K-1 )*3+2 )
               T1 = VECS( ( K-1 )*3+3 )
               T2 = T1*V2
               T3 = T1*V3
               DO 90 J = ITMP1, ITMP2
                  SUM = A( J, ICOL1 ) + V2*A( J, ICOL1+1 ) +
     $                  V3*A( J, ICOL1+2 )
                  A( J, ICOL1 ) = A( J, ICOL1 ) - SUM*T1
                  A( J, ICOL1+1 ) = A( J, ICOL1+1 ) - SUM*T2
                  A( J, ICOL1+2 ) = A( J, ICOL1+2 ) - SUM*T3
   90          CONTINUE
               IF( WANTZ ) THEN
                  DO 100 J = LILOZ, LIHIZ
                     SUM = Z( J, ICOL1 ) + V2*Z( J, ICOL1+1 ) +
     $                     V3*Z( J, ICOL1+2 )
                     Z( J, ICOL1 ) = Z( J, ICOL1 ) - SUM*T1
                     Z( J, ICOL1+1 ) = Z( J, ICOL1+1 ) - SUM*T2
                     Z( J, ICOL1+2 ) = Z( J, ICOL1+2 ) - SUM*T3
  100             CONTINUE
               END IF
               ICOL1 = ICOL1 + 1
  110       CONTINUE
         ELSE
            DO 120 J = ITMP1, ITMP2
               SUM = A( J, ICOL1 ) + V2*A( J, ICOL1+1 ) +
     $               V3*A( J, ICOL1+2 )
               A( J, ICOL1 ) = A( J, ICOL1 ) - SUM*T1
               A( J, ICOL1+1 ) = A( J, ICOL1+1 ) - SUM*T2
               A( J, ICOL1+2 ) = A( J, ICOL1+2 ) - SUM*T3
  120       CONTINUE
         END IF
      END IF
      RETURN
*
*     End of DLAREF
*
      END
*