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