Statistiques
| Révision :

root / src / lapack / double / dlarfb.f @ 10

Historique | Voir | Annoter | Télécharger (18,81 ko)

1 1 pfleura2
      SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
2 1 pfleura2
     $                   T, LDT, C, LDC, WORK, LDWORK )
3 1 pfleura2
      IMPLICIT NONE
4 1 pfleura2
*
5 1 pfleura2
*  -- LAPACK auxiliary routine (version 3.2) --
6 1 pfleura2
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
7 1 pfleura2
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 1 pfleura2
*     November 2006
9 1 pfleura2
*
10 1 pfleura2
*     .. Scalar Arguments ..
11 1 pfleura2
      CHARACTER          DIRECT, SIDE, STOREV, TRANS
12 1 pfleura2
      INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
13 1 pfleura2
*     ..
14 1 pfleura2
*     .. Array Arguments ..
15 1 pfleura2
      DOUBLE PRECISION   C( LDC, * ), T( LDT, * ), V( LDV, * ),
16 1 pfleura2
     $                   WORK( LDWORK, * )
17 1 pfleura2
*     ..
18 1 pfleura2
*
19 1 pfleura2
*  Purpose
20 1 pfleura2
*  =======
21 1 pfleura2
*
22 1 pfleura2
*  DLARFB applies a real block reflector H or its transpose H' to a
23 1 pfleura2
*  real m by n matrix C, from either the left or the right.
24 1 pfleura2
*
25 1 pfleura2
*  Arguments
26 1 pfleura2
*  =========
27 1 pfleura2
*
28 1 pfleura2
*  SIDE    (input) CHARACTER*1
29 1 pfleura2
*          = 'L': apply H or H' from the Left
30 1 pfleura2
*          = 'R': apply H or H' from the Right
31 1 pfleura2
*
32 1 pfleura2
*  TRANS   (input) CHARACTER*1
33 1 pfleura2
*          = 'N': apply H (No transpose)
34 1 pfleura2
*          = 'T': apply H' (Transpose)
35 1 pfleura2
*
36 1 pfleura2
*  DIRECT  (input) CHARACTER*1
37 1 pfleura2
*          Indicates how H is formed from a product of elementary
38 1 pfleura2
*          reflectors
39 1 pfleura2
*          = 'F': H = H(1) H(2) . . . H(k) (Forward)
40 1 pfleura2
*          = 'B': H = H(k) . . . H(2) H(1) (Backward)
41 1 pfleura2
*
42 1 pfleura2
*  STOREV  (input) CHARACTER*1
43 1 pfleura2
*          Indicates how the vectors which define the elementary
44 1 pfleura2
*          reflectors are stored:
45 1 pfleura2
*          = 'C': Columnwise
46 1 pfleura2
*          = 'R': Rowwise
47 1 pfleura2
*
48 1 pfleura2
*  M       (input) INTEGER
49 1 pfleura2
*          The number of rows of the matrix C.
50 1 pfleura2
*
51 1 pfleura2
*  N       (input) INTEGER
52 1 pfleura2
*          The number of columns of the matrix C.
53 1 pfleura2
*
54 1 pfleura2
*  K       (input) INTEGER
55 1 pfleura2
*          The order of the matrix T (= the number of elementary
56 1 pfleura2
*          reflectors whose product defines the block reflector).
57 1 pfleura2
*
58 1 pfleura2
*  V       (input) DOUBLE PRECISION array, dimension
59 1 pfleura2
*                                (LDV,K) if STOREV = 'C'
60 1 pfleura2
*                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
61 1 pfleura2
*                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
62 1 pfleura2
*          The matrix V. See further details.
63 1 pfleura2
*
64 1 pfleura2
*  LDV     (input) INTEGER
65 1 pfleura2
*          The leading dimension of the array V.
66 1 pfleura2
*          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
67 1 pfleura2
*          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
68 1 pfleura2
*          if STOREV = 'R', LDV >= K.
69 1 pfleura2
*
70 1 pfleura2
*  T       (input) DOUBLE PRECISION array, dimension (LDT,K)
71 1 pfleura2
*          The triangular k by k matrix T in the representation of the
72 1 pfleura2
*          block reflector.
73 1 pfleura2
*
74 1 pfleura2
*  LDT     (input) INTEGER
75 1 pfleura2
*          The leading dimension of the array T. LDT >= K.
76 1 pfleura2
*
77 1 pfleura2
*  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
78 1 pfleura2
*          On entry, the m by n matrix C.
79 1 pfleura2
*          On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
80 1 pfleura2
*
81 1 pfleura2
*  LDC     (input) INTEGER
82 1 pfleura2
*          The leading dimension of the array C. LDA >= max(1,M).
83 1 pfleura2
*
84 1 pfleura2
*  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
85 1 pfleura2
*
86 1 pfleura2
*  LDWORK  (input) INTEGER
87 1 pfleura2
*          The leading dimension of the array WORK.
88 1 pfleura2
*          If SIDE = 'L', LDWORK >= max(1,N);
89 1 pfleura2
*          if SIDE = 'R', LDWORK >= max(1,M).
90 1 pfleura2
*
91 1 pfleura2
*  =====================================================================
92 1 pfleura2
*
93 1 pfleura2
*     .. Parameters ..
94 1 pfleura2
      DOUBLE PRECISION   ONE
95 1 pfleura2
      PARAMETER          ( ONE = 1.0D+0 )
96 1 pfleura2
*     ..
97 1 pfleura2
*     .. Local Scalars ..
98 1 pfleura2
      CHARACTER          TRANST
99 1 pfleura2
      INTEGER            I, J, LASTV, LASTC
100 1 pfleura2
*     ..
101 1 pfleura2
*     .. External Functions ..
102 1 pfleura2
      LOGICAL            LSAME
103 1 pfleura2
      INTEGER            ILADLR, ILADLC
104 1 pfleura2
      EXTERNAL           LSAME, ILADLR, ILADLC
105 1 pfleura2
*     ..
106 1 pfleura2
*     .. External Subroutines ..
107 1 pfleura2
      EXTERNAL           DCOPY, DGEMM, DTRMM
108 1 pfleura2
*     ..
109 1 pfleura2
*     .. Executable Statements ..
110 1 pfleura2
*
111 1 pfleura2
*     Quick return if possible
112 1 pfleura2
*
113 1 pfleura2
      IF( M.LE.0 .OR. N.LE.0 )
114 1 pfleura2
     $   RETURN
115 1 pfleura2
*
116 1 pfleura2
      IF( LSAME( TRANS, 'N' ) ) THEN
117 1 pfleura2
         TRANST = 'T'
118 1 pfleura2
      ELSE
119 1 pfleura2
         TRANST = 'N'
120 1 pfleura2
      END IF
121 1 pfleura2
*
122 1 pfleura2
      IF( LSAME( STOREV, 'C' ) ) THEN
123 1 pfleura2
*
124 1 pfleura2
         IF( LSAME( DIRECT, 'F' ) ) THEN
125 1 pfleura2
*
126 1 pfleura2
*           Let  V =  ( V1 )    (first K rows)
127 1 pfleura2
*                     ( V2 )
128 1 pfleura2
*           where  V1  is unit lower triangular.
129 1 pfleura2
*
130 1 pfleura2
            IF( LSAME( SIDE, 'L' ) ) THEN
131 1 pfleura2
*
132 1 pfleura2
*              Form  H * C  or  H' * C  where  C = ( C1 )
133 1 pfleura2
*                                                  ( C2 )
134 1 pfleura2
*
135 1 pfleura2
               LASTV = MAX( K, ILADLR( M, K, V, LDV ) )
136 1 pfleura2
               LASTC = ILADLC( LASTV, N, C, LDC )
137 1 pfleura2
*
138 1 pfleura2
*              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK)
139 1 pfleura2
*
140 1 pfleura2
*              W := C1'
141 1 pfleura2
*
142 1 pfleura2
               DO 10 J = 1, K
143 1 pfleura2
                  CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
144 1 pfleura2
   10          CONTINUE
145 1 pfleura2
*
146 1 pfleura2
*              W := W * V1
147 1 pfleura2
*
148 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
149 1 pfleura2
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
150 1 pfleura2
               IF( LASTV.GT.K ) THEN
151 1 pfleura2
*
152 1 pfleura2
*                 W := W + C2'*V2
153 1 pfleura2
*
154 1 pfleura2
                  CALL DGEMM( 'Transpose', 'No transpose',
155 1 pfleura2
     $                 LASTC, K, LASTV-K,
156 1 pfleura2
     $                 ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV,
157 1 pfleura2
     $                 ONE, WORK, LDWORK )
158 1 pfleura2
               END IF
159 1 pfleura2
*
160 1 pfleura2
*              W := W * T'  or  W * T
161 1 pfleura2
*
162 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
163 1 pfleura2
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
164 1 pfleura2
*
165 1 pfleura2
*              C := C - V * W'
166 1 pfleura2
*
167 1 pfleura2
               IF( LASTV.GT.K ) THEN
168 1 pfleura2
*
169 1 pfleura2
*                 C2 := C2 - V2 * W'
170 1 pfleura2
*
171 1 pfleura2
                  CALL DGEMM( 'No transpose', 'Transpose',
172 1 pfleura2
     $                 LASTV-K, LASTC, K,
173 1 pfleura2
     $                 -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE,
174 1 pfleura2
     $                 C( K+1, 1 ), LDC )
175 1 pfleura2
               END IF
176 1 pfleura2
*
177 1 pfleura2
*              W := W * V1'
178 1 pfleura2
*
179 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
180 1 pfleura2
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
181 1 pfleura2
*
182 1 pfleura2
*              C1 := C1 - W'
183 1 pfleura2
*
184 1 pfleura2
               DO 30 J = 1, K
185 1 pfleura2
                  DO 20 I = 1, LASTC
186 1 pfleura2
                     C( J, I ) = C( J, I ) - WORK( I, J )
187 1 pfleura2
   20             CONTINUE
188 1 pfleura2
   30          CONTINUE
189 1 pfleura2
*
190 1 pfleura2
            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
191 1 pfleura2
*
192 1 pfleura2
*              Form  C * H  or  C * H'  where  C = ( C1  C2 )
193 1 pfleura2
*
194 1 pfleura2
               LASTV = MAX( K, ILADLR( N, K, V, LDV ) )
195 1 pfleura2
               LASTC = ILADLR( M, LASTV, C, LDC )
196 1 pfleura2
*
197 1 pfleura2
*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
198 1 pfleura2
*
199 1 pfleura2
*              W := C1
200 1 pfleura2
*
201 1 pfleura2
               DO 40 J = 1, K
202 1 pfleura2
                  CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
203 1 pfleura2
   40          CONTINUE
204 1 pfleura2
*
205 1 pfleura2
*              W := W * V1
206 1 pfleura2
*
207 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
208 1 pfleura2
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
209 1 pfleura2
               IF( LASTV.GT.K ) THEN
210 1 pfleura2
*
211 1 pfleura2
*                 W := W + C2 * V2
212 1 pfleura2
*
213 1 pfleura2
                  CALL DGEMM( 'No transpose', 'No transpose',
214 1 pfleura2
     $                 LASTC, K, LASTV-K,
215 1 pfleura2
     $                 ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
216 1 pfleura2
     $                 ONE, WORK, LDWORK )
217 1 pfleura2
               END IF
218 1 pfleura2
*
219 1 pfleura2
*              W := W * T  or  W * T'
220 1 pfleura2
*
221 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
222 1 pfleura2
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
223 1 pfleura2
*
224 1 pfleura2
*              C := C - W * V'
225 1 pfleura2
*
226 1 pfleura2
               IF( LASTV.GT.K ) THEN
227 1 pfleura2
*
228 1 pfleura2
*                 C2 := C2 - W * V2'
229 1 pfleura2
*
230 1 pfleura2
                  CALL DGEMM( 'No transpose', 'Transpose',
231 1 pfleura2
     $                 LASTC, LASTV-K, K,
232 1 pfleura2
     $                 -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE,
233 1 pfleura2
     $                 C( 1, K+1 ), LDC )
234 1 pfleura2
               END IF
235 1 pfleura2
*
236 1 pfleura2
*              W := W * V1'
237 1 pfleura2
*
238 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
239 1 pfleura2
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
240 1 pfleura2
*
241 1 pfleura2
*              C1 := C1 - W
242 1 pfleura2
*
243 1 pfleura2
               DO 60 J = 1, K
244 1 pfleura2
                  DO 50 I = 1, LASTC
245 1 pfleura2
                     C( I, J ) = C( I, J ) - WORK( I, J )
246 1 pfleura2
   50             CONTINUE
247 1 pfleura2
   60          CONTINUE
248 1 pfleura2
            END IF
249 1 pfleura2
*
250 1 pfleura2
         ELSE
251 1 pfleura2
*
252 1 pfleura2
*           Let  V =  ( V1 )
253 1 pfleura2
*                     ( V2 )    (last K rows)
254 1 pfleura2
*           where  V2  is unit upper triangular.
255 1 pfleura2
*
256 1 pfleura2
            IF( LSAME( SIDE, 'L' ) ) THEN
257 1 pfleura2
*
258 1 pfleura2
*              Form  H * C  or  H' * C  where  C = ( C1 )
259 1 pfleura2
*                                                  ( C2 )
260 1 pfleura2
*
261 1 pfleura2
               LASTV = MAX( K, ILADLR( M, K, V, LDV ) )
262 1 pfleura2
               LASTC = ILADLC( LASTV, N, C, LDC )
263 1 pfleura2
*
264 1 pfleura2
*              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK)
265 1 pfleura2
*
266 1 pfleura2
*              W := C2'
267 1 pfleura2
*
268 1 pfleura2
               DO 70 J = 1, K
269 1 pfleura2
                  CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
270 1 pfleura2
     $                 WORK( 1, J ), 1 )
271 1 pfleura2
   70          CONTINUE
272 1 pfleura2
*
273 1 pfleura2
*              W := W * V2
274 1 pfleura2
*
275 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
276 1 pfleura2
     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
277 1 pfleura2
     $              WORK, LDWORK )
278 1 pfleura2
               IF( LASTV.GT.K ) THEN
279 1 pfleura2
*
280 1 pfleura2
*                 W := W + C1'*V1
281 1 pfleura2
*
282 1 pfleura2
                  CALL DGEMM( 'Transpose', 'No transpose',
283 1 pfleura2
     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
284 1 pfleura2
     $                 ONE, WORK, LDWORK )
285 1 pfleura2
               END IF
286 1 pfleura2
*
287 1 pfleura2
*              W := W * T'  or  W * T
288 1 pfleura2
*
289 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
290 1 pfleura2
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
291 1 pfleura2
*
292 1 pfleura2
*              C := C - V * W'
293 1 pfleura2
*
294 1 pfleura2
               IF( LASTV.GT.K ) THEN
295 1 pfleura2
*
296 1 pfleura2
*                 C1 := C1 - V1 * W'
297 1 pfleura2
*
298 1 pfleura2
                  CALL DGEMM( 'No transpose', 'Transpose',
299 1 pfleura2
     $                 LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
300 1 pfleura2
     $                 ONE, C, LDC )
301 1 pfleura2
               END IF
302 1 pfleura2
*
303 1 pfleura2
*              W := W * V2'
304 1 pfleura2
*
305 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
306 1 pfleura2
     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
307 1 pfleura2
     $              WORK, LDWORK )
308 1 pfleura2
*
309 1 pfleura2
*              C2 := C2 - W'
310 1 pfleura2
*
311 1 pfleura2
               DO 90 J = 1, K
312 1 pfleura2
                  DO 80 I = 1, LASTC
313 1 pfleura2
                     C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J)
314 1 pfleura2
   80             CONTINUE
315 1 pfleura2
   90          CONTINUE
316 1 pfleura2
*
317 1 pfleura2
            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
318 1 pfleura2
*
319 1 pfleura2
*              Form  C * H  or  C * H'  where  C = ( C1  C2 )
320 1 pfleura2
*
321 1 pfleura2
               LASTV = MAX( K, ILADLR( N, K, V, LDV ) )
322 1 pfleura2
               LASTC = ILADLR( M, LASTV, C, LDC )
323 1 pfleura2
*
324 1 pfleura2
*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
325 1 pfleura2
*
326 1 pfleura2
*              W := C2
327 1 pfleura2
*
328 1 pfleura2
               DO 100 J = 1, K
329 1 pfleura2
                  CALL DCOPY( LASTC, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
330 1 pfleura2
  100          CONTINUE
331 1 pfleura2
*
332 1 pfleura2
*              W := W * V2
333 1 pfleura2
*
334 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
335 1 pfleura2
     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
336 1 pfleura2
     $              WORK, LDWORK )
337 1 pfleura2
               IF( LASTV.GT.K ) THEN
338 1 pfleura2
*
339 1 pfleura2
*                 W := W + C1 * V1
340 1 pfleura2
*
341 1 pfleura2
                  CALL DGEMM( 'No transpose', 'No transpose',
342 1 pfleura2
     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
343 1 pfleura2
     $                 ONE, WORK, LDWORK )
344 1 pfleura2
               END IF
345 1 pfleura2
*
346 1 pfleura2
*              W := W * T  or  W * T'
347 1 pfleura2
*
348 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
349 1 pfleura2
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
350 1 pfleura2
*
351 1 pfleura2
*              C := C - W * V'
352 1 pfleura2
*
353 1 pfleura2
               IF( LASTV.GT.K ) THEN
354 1 pfleura2
*
355 1 pfleura2
*                 C1 := C1 - W * V1'
356 1 pfleura2
*
357 1 pfleura2
                  CALL DGEMM( 'No transpose', 'Transpose',
358 1 pfleura2
     $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
359 1 pfleura2
     $                 ONE, C, LDC )
360 1 pfleura2
               END IF
361 1 pfleura2
*
362 1 pfleura2
*              W := W * V2'
363 1 pfleura2
*
364 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
365 1 pfleura2
     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
366 1 pfleura2
     $              WORK, LDWORK )
367 1 pfleura2
*
368 1 pfleura2
*              C2 := C2 - W
369 1 pfleura2
*
370 1 pfleura2
               DO 120 J = 1, K
371 1 pfleura2
                  DO 110 I = 1, LASTC
372 1 pfleura2
                     C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J)
373 1 pfleura2
  110             CONTINUE
374 1 pfleura2
  120          CONTINUE
375 1 pfleura2
            END IF
376 1 pfleura2
         END IF
377 1 pfleura2
*
378 1 pfleura2
      ELSE IF( LSAME( STOREV, 'R' ) ) THEN
379 1 pfleura2
*
380 1 pfleura2
         IF( LSAME( DIRECT, 'F' ) ) THEN
381 1 pfleura2
*
382 1 pfleura2
*           Let  V =  ( V1  V2 )    (V1: first K columns)
383 1 pfleura2
*           where  V1  is unit upper triangular.
384 1 pfleura2
*
385 1 pfleura2
            IF( LSAME( SIDE, 'L' ) ) THEN
386 1 pfleura2
*
387 1 pfleura2
*              Form  H * C  or  H' * C  where  C = ( C1 )
388 1 pfleura2
*                                                  ( C2 )
389 1 pfleura2
*
390 1 pfleura2
               LASTV = MAX( K, ILADLC( K, M, V, LDV ) )
391 1 pfleura2
               LASTC = ILADLC( LASTV, N, C, LDC )
392 1 pfleura2
*
393 1 pfleura2
*              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK)
394 1 pfleura2
*
395 1 pfleura2
*              W := C1'
396 1 pfleura2
*
397 1 pfleura2
               DO 130 J = 1, K
398 1 pfleura2
                  CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
399 1 pfleura2
  130          CONTINUE
400 1 pfleura2
*
401 1 pfleura2
*              W := W * V1'
402 1 pfleura2
*
403 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
404 1 pfleura2
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
405 1 pfleura2
               IF( LASTV.GT.K ) THEN
406 1 pfleura2
*
407 1 pfleura2
*                 W := W + C2'*V2'
408 1 pfleura2
*
409 1 pfleura2
                  CALL DGEMM( 'Transpose', 'Transpose',
410 1 pfleura2
     $                 LASTC, K, LASTV-K,
411 1 pfleura2
     $                 ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV,
412 1 pfleura2
     $                 ONE, WORK, LDWORK )
413 1 pfleura2
               END IF
414 1 pfleura2
*
415 1 pfleura2
*              W := W * T'  or  W * T
416 1 pfleura2
*
417 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
418 1 pfleura2
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
419 1 pfleura2
*
420 1 pfleura2
*              C := C - V' * W'
421 1 pfleura2
*
422 1 pfleura2
               IF( LASTV.GT.K ) THEN
423 1 pfleura2
*
424 1 pfleura2
*                 C2 := C2 - V2' * W'
425 1 pfleura2
*
426 1 pfleura2
                  CALL DGEMM( 'Transpose', 'Transpose',
427 1 pfleura2
     $                 LASTV-K, LASTC, K,
428 1 pfleura2
     $                 -ONE, V( 1, K+1 ), LDV, WORK, LDWORK,
429 1 pfleura2
     $                 ONE, C( K+1, 1 ), LDC )
430 1 pfleura2
               END IF
431 1 pfleura2
*
432 1 pfleura2
*              W := W * V1
433 1 pfleura2
*
434 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
435 1 pfleura2
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
436 1 pfleura2
*
437 1 pfleura2
*              C1 := C1 - W'
438 1 pfleura2
*
439 1 pfleura2
               DO 150 J = 1, K
440 1 pfleura2
                  DO 140 I = 1, LASTC
441 1 pfleura2
                     C( J, I ) = C( J, I ) - WORK( I, J )
442 1 pfleura2
  140             CONTINUE
443 1 pfleura2
  150          CONTINUE
444 1 pfleura2
*
445 1 pfleura2
            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
446 1 pfleura2
*
447 1 pfleura2
*              Form  C * H  or  C * H'  where  C = ( C1  C2 )
448 1 pfleura2
*
449 1 pfleura2
               LASTV = MAX( K, ILADLC( K, N, V, LDV ) )
450 1 pfleura2
               LASTC = ILADLR( M, LASTV, C, LDC )
451 1 pfleura2
*
452 1 pfleura2
*              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK)
453 1 pfleura2
*
454 1 pfleura2
*              W := C1
455 1 pfleura2
*
456 1 pfleura2
               DO 160 J = 1, K
457 1 pfleura2
                  CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
458 1 pfleura2
  160          CONTINUE
459 1 pfleura2
*
460 1 pfleura2
*              W := W * V1'
461 1 pfleura2
*
462 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
463 1 pfleura2
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
464 1 pfleura2
               IF( LASTV.GT.K ) THEN
465 1 pfleura2
*
466 1 pfleura2
*                 W := W + C2 * V2'
467 1 pfleura2
*
468 1 pfleura2
                  CALL DGEMM( 'No transpose', 'Transpose',
469 1 pfleura2
     $                 LASTC, K, LASTV-K,
470 1 pfleura2
     $                 ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV,
471 1 pfleura2
     $                 ONE, WORK, LDWORK )
472 1 pfleura2
               END IF
473 1 pfleura2
*
474 1 pfleura2
*              W := W * T  or  W * T'
475 1 pfleura2
*
476 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
477 1 pfleura2
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
478 1 pfleura2
*
479 1 pfleura2
*              C := C - W * V
480 1 pfleura2
*
481 1 pfleura2
               IF( LASTV.GT.K ) THEN
482 1 pfleura2
*
483 1 pfleura2
*                 C2 := C2 - W * V2
484 1 pfleura2
*
485 1 pfleura2
                  CALL DGEMM( 'No transpose', 'No transpose',
486 1 pfleura2
     $                 LASTC, LASTV-K, K,
487 1 pfleura2
     $                 -ONE, WORK, LDWORK, V( 1, K+1 ), LDV,
488 1 pfleura2
     $                 ONE, C( 1, K+1 ), LDC )
489 1 pfleura2
               END IF
490 1 pfleura2
*
491 1 pfleura2
*              W := W * V1
492 1 pfleura2
*
493 1 pfleura2
               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
494 1 pfleura2
     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
495 1 pfleura2
*
496 1 pfleura2
*              C1 := C1 - W
497 1 pfleura2
*
498 1 pfleura2
               DO 180 J = 1, K
499 1 pfleura2
                  DO 170 I = 1, LASTC
500 1 pfleura2
                     C( I, J ) = C( I, J ) - WORK( I, J )
501 1 pfleura2
  170             CONTINUE
502 1 pfleura2
  180          CONTINUE
503 1 pfleura2
*
504 1 pfleura2
            END IF
505 1 pfleura2
*
506 1 pfleura2
         ELSE
507 1 pfleura2
*
508 1 pfleura2
*           Let  V =  ( V1  V2 )    (V2: last K columns)
509 1 pfleura2
*           where  V2  is unit lower triangular.
510 1 pfleura2
*
511 1 pfleura2
            IF( LSAME( SIDE, 'L' ) ) THEN
512 1 pfleura2
*
513 1 pfleura2
*              Form  H * C  or  H' * C  where  C = ( C1 )
514 1 pfleura2
*                                                  ( C2 )
515 1 pfleura2
*
516 1 pfleura2
               LASTV = MAX( K, ILADLC( K, M, V, LDV ) )
517 1 pfleura2
               LASTC = ILADLC( LASTV, N, C, LDC )
518 1 pfleura2
*
519 1 pfleura2
*              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK)
520 1 pfleura2
*
521 1 pfleura2
*              W := C2'
522 1 pfleura2
*
523 1 pfleura2
               DO 190 J = 1, K
524 1 pfleura2
                  CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
525 1 pfleura2
     $                 WORK( 1, J ), 1 )
526 1 pfleura2
  190          CONTINUE
527 1 pfleura2
*
528 1 pfleura2
*              W := W * V2'
529 1 pfleura2
*
530 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
531 1 pfleura2
     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
532 1 pfleura2
     $              WORK, LDWORK )
533 1 pfleura2
               IF( LASTV.GT.K ) THEN
534 1 pfleura2
*
535 1 pfleura2
*                 W := W + C1'*V1'
536 1 pfleura2
*
537 1 pfleura2
                  CALL DGEMM( 'Transpose', 'Transpose',
538 1 pfleura2
     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
539 1 pfleura2
     $                 ONE, WORK, LDWORK )
540 1 pfleura2
               END IF
541 1 pfleura2
*
542 1 pfleura2
*              W := W * T'  or  W * T
543 1 pfleura2
*
544 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
545 1 pfleura2
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
546 1 pfleura2
*
547 1 pfleura2
*              C := C - V' * W'
548 1 pfleura2
*
549 1 pfleura2
               IF( LASTV.GT.K ) THEN
550 1 pfleura2
*
551 1 pfleura2
*                 C1 := C1 - V1' * W'
552 1 pfleura2
*
553 1 pfleura2
                  CALL DGEMM( 'Transpose', 'Transpose',
554 1 pfleura2
     $                 LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
555 1 pfleura2
     $                 ONE, C, LDC )
556 1 pfleura2
               END IF
557 1 pfleura2
*
558 1 pfleura2
*              W := W * V2
559 1 pfleura2
*
560 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
561 1 pfleura2
     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
562 1 pfleura2
     $              WORK, LDWORK )
563 1 pfleura2
*
564 1 pfleura2
*              C2 := C2 - W'
565 1 pfleura2
*
566 1 pfleura2
               DO 210 J = 1, K
567 1 pfleura2
                  DO 200 I = 1, LASTC
568 1 pfleura2
                     C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J)
569 1 pfleura2
  200             CONTINUE
570 1 pfleura2
  210          CONTINUE
571 1 pfleura2
*
572 1 pfleura2
            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
573 1 pfleura2
*
574 1 pfleura2
*              Form  C * H  or  C * H'  where  C = ( C1  C2 )
575 1 pfleura2
*
576 1 pfleura2
               LASTV = MAX( K, ILADLC( K, N, V, LDV ) )
577 1 pfleura2
               LASTC = ILADLR( M, LASTV, C, LDC )
578 1 pfleura2
*
579 1 pfleura2
*              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK)
580 1 pfleura2
*
581 1 pfleura2
*              W := C2
582 1 pfleura2
*
583 1 pfleura2
               DO 220 J = 1, K
584 1 pfleura2
                  CALL DCOPY( LASTC, C( 1, LASTV-K+J ), 1,
585 1 pfleura2
     $                 WORK( 1, J ), 1 )
586 1 pfleura2
  220          CONTINUE
587 1 pfleura2
*
588 1 pfleura2
*              W := W * V2'
589 1 pfleura2
*
590 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
591 1 pfleura2
     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
592 1 pfleura2
     $              WORK, LDWORK )
593 1 pfleura2
               IF( LASTV.GT.K ) THEN
594 1 pfleura2
*
595 1 pfleura2
*                 W := W + C1 * V1'
596 1 pfleura2
*
597 1 pfleura2
                  CALL DGEMM( 'No transpose', 'Transpose',
598 1 pfleura2
     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
599 1 pfleura2
     $                 ONE, WORK, LDWORK )
600 1 pfleura2
               END IF
601 1 pfleura2
*
602 1 pfleura2
*              W := W * T  or  W * T'
603 1 pfleura2
*
604 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
605 1 pfleura2
     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
606 1 pfleura2
*
607 1 pfleura2
*              C := C - W * V
608 1 pfleura2
*
609 1 pfleura2
               IF( LASTV.GT.K ) THEN
610 1 pfleura2
*
611 1 pfleura2
*                 C1 := C1 - W * V1
612 1 pfleura2
*
613 1 pfleura2
                  CALL DGEMM( 'No transpose', 'No transpose',
614 1 pfleura2
     $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
615 1 pfleura2
     $                 ONE, C, LDC )
616 1 pfleura2
               END IF
617 1 pfleura2
*
618 1 pfleura2
*              W := W * V2
619 1 pfleura2
*
620 1 pfleura2
               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
621 1 pfleura2
     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
622 1 pfleura2
     $              WORK, LDWORK )
623 1 pfleura2
*
624 1 pfleura2
*              C1 := C1 - W
625 1 pfleura2
*
626 1 pfleura2
               DO 240 J = 1, K
627 1 pfleura2
                  DO 230 I = 1, LASTC
628 1 pfleura2
                     C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J)
629 1 pfleura2
  230             CONTINUE
630 1 pfleura2
  240          CONTINUE
631 1 pfleura2
*
632 1 pfleura2
            END IF
633 1 pfleura2
*
634 1 pfleura2
         END IF
635 1 pfleura2
      END IF
636 1 pfleura2
*
637 1 pfleura2
      RETURN
638 1 pfleura2
*
639 1 pfleura2
*     End of DLARFB
640 1 pfleura2
*
641 1 pfleura2
      END