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