LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlarfb()

subroutine dlarfb ( character side,
character trans,
character direct,
character storev,
integer m,
integer n,
integer k,
double precision, dimension( ldv, * ) v,
integer ldv,
double precision, dimension( ldt, * ) t,
integer ldt,
double precision, dimension( ldc, * ) c,
integer ldc,
double precision, dimension( ldwork, * ) work,
integer ldwork )

DLARFB applies a block reflector or its transpose to a general rectangular matrix.

Download DLARFB + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DLARFB applies a real block reflector H or its transpose H**T to a
!> real m by n matrix C, from either the left or the right.
!>
Parameters
[in]SIDE
!>          SIDE is CHARACTER*1
!>          = 'L': apply H or H**T from the Left
!>          = 'R': apply H or H**T from the Right
!>
[in]TRANS
!>          TRANS is CHARACTER*1
!>          = 'N': apply H (No transpose)
!>          = 'T': apply H**T (Transpose)
!>
[in]DIRECT
!>          DIRECT is CHARACTER*1
!>          Indicates how H is formed from a product of elementary
!>          reflectors
!>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
!>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
!>
[in]STOREV
!>          STOREV is CHARACTER*1
!>          Indicates how the vectors which define the elementary
!>          reflectors are stored:
!>          = 'C': Columnwise
!>          = 'R': Rowwise
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix C.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix C.
!>
[in]K
!>          K is INTEGER
!>          The order of the matrix T (= the number of elementary
!>          reflectors whose product defines the block reflector).
!>          If SIDE = 'L', M >= K >= 0;
!>          if SIDE = 'R', N >= K >= 0.
!>
[in]V
!>          V is DOUBLE PRECISION array, dimension
!>                                (LDV,K) if STOREV = 'C'
!>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
!>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
!>          The matrix V. See Further Details.
!>
[in]LDV
!>          LDV is INTEGER
!>          The leading dimension of the array V.
!>          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
!>          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
!>          if STOREV = 'R', LDV >= K.
!>
[in]T
!>          T is DOUBLE PRECISION array, dimension (LDT,K)
!>          The triangular k by k matrix T in the representation of the
!>          block reflector.
!>
[in]LDT
!>          LDT is INTEGER
!>          The leading dimension of the array T. LDT >= K.
!>
[in,out]C
!>          C is DOUBLE PRECISION array, dimension (LDC,N)
!>          On entry, the m by n matrix C.
!>          On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.
!>
[in]LDC
!>          LDC is INTEGER
!>          The leading dimension of the array C. LDC >= max(1,M).
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (LDWORK,K)
!>
[in]LDWORK
!>          LDWORK is INTEGER
!>          The leading dimension of the array WORK.
!>          If SIDE = 'L', LDWORK >= max(1,N);
!>          if SIDE = 'R', LDWORK >= max(1,M).
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The shape of the matrix V and the storage of the vectors which define
!>  the H(i) is best illustrated by the following example with n = 5 and
!>  k = 3. The triangular part of V (including its diagonal) is not
!>  referenced.
!>
!>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
!>
!>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
!>                   ( v1  1    )                     (     1 v2 v2 v2 )
!>                   ( v1 v2  1 )                     (        1 v3 v3 )
!>                   ( v1 v2 v3 )
!>                   ( v1 v2 v3 )
!>
!>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
!>
!>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
!>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
!>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
!>                   (     1 v3 )
!>                   (        1 )
!>

Definition at line 192 of file dlarfb.f.

195*
196* -- LAPACK auxiliary routine --
197* -- LAPACK is a software package provided by Univ. of Tennessee, --
198* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
199*
200* .. Scalar Arguments ..
201 CHARACTER DIRECT, SIDE, STOREV, TRANS
202 INTEGER K, LDC, LDT, LDV, LDWORK, M, N
203* ..
204* .. Array Arguments ..
205 DOUBLE PRECISION C( LDC, * ), T( LDT, * ), V( LDV, * ),
206 $ WORK( LDWORK, * )
207* ..
208*
209* =====================================================================
210*
211* .. Parameters ..
212 DOUBLE PRECISION ONE
213 parameter( one = 1.0d+0 )
214* ..
215* .. Local Scalars ..
216 CHARACTER TRANST
217 INTEGER I, J
218* ..
219* .. External Functions ..
220 LOGICAL LSAME
221 EXTERNAL lsame
222* ..
223* .. External Subroutines ..
224 EXTERNAL dcopy, dgemm, dtrmm
225* ..
226* .. Executable Statements ..
227*
228* Quick return if possible
229*
230 IF( m.LE.0 .OR. n.LE.0 )
231 $ RETURN
232*
233 IF( lsame( trans, 'N' ) ) THEN
234 transt = 'T'
235 ELSE
236 transt = 'N'
237 END IF
238*
239 IF( lsame( storev, 'C' ) ) THEN
240*
241 IF( lsame( direct, 'F' ) ) THEN
242*
243* Let V = ( V1 ) (first K rows)
244* ( V2 )
245* where V1 is unit lower triangular.
246*
247 IF( lsame( side, 'L' ) ) THEN
248*
249* Form H * C or H**T * C where C = ( C1 )
250* ( C2 )
251*
252* W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK)
253*
254* W := C1**T
255*
256 DO 10 j = 1, k
257 CALL dcopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
258 10 CONTINUE
259*
260* W := W * V1
261*
262 CALL dtrmm( 'Right', 'Lower', 'No transpose', 'Unit',
263 $ n,
264 $ k, one, v, ldv, work, ldwork )
265 IF( m.GT.k ) THEN
266*
267* W := W + C2**T * V2
268*
269 CALL dgemm( 'Transpose', 'No transpose', n, k, m-k,
270 $ one, c( k+1, 1 ), ldc, v( k+1, 1 ), ldv,
271 $ one, work, ldwork )
272 END IF
273*
274* W := W * T**T or W * T
275*
276 CALL dtrmm( 'Right', 'Upper', transt, 'Non-unit', n,
277 $ k,
278 $ one, t, ldt, work, ldwork )
279*
280* C := C - V * W**T
281*
282 IF( m.GT.k ) THEN
283*
284* C2 := C2 - V2 * W**T
285*
286 CALL dgemm( 'No transpose', 'Transpose', m-k, n, k,
287 $ -one, v( k+1, 1 ), ldv, work, ldwork, one,
288 $ c( k+1, 1 ), ldc )
289 END IF
290*
291* W := W * V1**T
292*
293 CALL dtrmm( 'Right', 'Lower', 'Transpose', 'Unit', n,
294 $ k,
295 $ one, v, ldv, work, ldwork )
296*
297* C1 := C1 - W**T
298*
299 DO 30 j = 1, k
300 DO 20 i = 1, n
301 c( j, i ) = c( j, i ) - work( i, j )
302 20 CONTINUE
303 30 CONTINUE
304*
305 ELSE IF( lsame( side, 'R' ) ) THEN
306*
307* Form C * H or C * H**T where C = ( C1 C2 )
308*
309* W := C * V = (C1*V1 + C2*V2) (stored in WORK)
310*
311* W := C1
312*
313 DO 40 j = 1, k
314 CALL dcopy( m, c( 1, j ), 1, work( 1, j ), 1 )
315 40 CONTINUE
316*
317* W := W * V1
318*
319 CALL dtrmm( 'Right', 'Lower', 'No transpose', 'Unit',
320 $ m,
321 $ k, one, v, ldv, work, ldwork )
322 IF( n.GT.k ) THEN
323*
324* W := W + C2 * V2
325*
326 CALL dgemm( 'No transpose', 'No transpose', m, k,
327 $ n-k,
328 $ one, c( 1, k+1 ), ldc, v( k+1, 1 ), ldv,
329 $ one, work, ldwork )
330 END IF
331*
332* W := W * T or W * T**T
333*
334 CALL dtrmm( 'Right', 'Upper', trans, 'Non-unit', m, k,
335 $ one, t, ldt, work, ldwork )
336*
337* C := C - W * V**T
338*
339 IF( n.GT.k ) THEN
340*
341* C2 := C2 - W * V2**T
342*
343 CALL dgemm( 'No transpose', 'Transpose', m, n-k, k,
344 $ -one, work, ldwork, v( k+1, 1 ), ldv, one,
345 $ c( 1, k+1 ), ldc )
346 END IF
347*
348* W := W * V1**T
349*
350 CALL dtrmm( 'Right', 'Lower', 'Transpose', 'Unit', m,
351 $ k,
352 $ one, v, ldv, work, ldwork )
353*
354* C1 := C1 - W
355*
356 DO 60 j = 1, k
357 DO 50 i = 1, m
358 c( i, j ) = c( i, j ) - work( i, j )
359 50 CONTINUE
360 60 CONTINUE
361 END IF
362*
363 ELSE
364*
365* Let V = ( V1 )
366* ( V2 ) (last K rows)
367* where V2 is unit upper triangular.
368*
369 IF( lsame( side, 'L' ) ) THEN
370*
371* Form H * C or H**T * C where C = ( C1 )
372* ( C2 )
373*
374* W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK)
375*
376* W := C2**T
377*
378 DO 70 j = 1, k
379 CALL dcopy( n, c( m-k+j, 1 ), ldc, work( 1, j ),
380 $ 1 )
381 70 CONTINUE
382*
383* W := W * V2
384*
385 CALL dtrmm( 'Right', 'Upper', 'No transpose', 'Unit',
386 $ n,
387 $ k, one, v( m-k+1, 1 ), ldv, work, ldwork )
388 IF( m.GT.k ) THEN
389*
390* W := W + C1**T * V1
391*
392 CALL dgemm( 'Transpose', 'No transpose', n, k, m-k,
393 $ one, c, ldc, v, ldv, one, work, ldwork )
394 END IF
395*
396* W := W * T**T or W * T
397*
398 CALL dtrmm( 'Right', 'Lower', transt, 'Non-unit', n,
399 $ k,
400 $ one, t, ldt, work, ldwork )
401*
402* C := C - V * W**T
403*
404 IF( m.GT.k ) THEN
405*
406* C1 := C1 - V1 * W**T
407*
408 CALL dgemm( 'No transpose', 'Transpose', m-k, n, k,
409 $ -one, v, ldv, work, ldwork, one, c, ldc )
410 END IF
411*
412* W := W * V2**T
413*
414 CALL dtrmm( 'Right', 'Upper', 'Transpose', 'Unit', n,
415 $ k,
416 $ one, v( m-k+1, 1 ), ldv, work, ldwork )
417*
418* C2 := C2 - W**T
419*
420 DO 90 j = 1, k
421 DO 80 i = 1, n
422 c( m-k+j, i ) = c( m-k+j, i ) - work( i, j )
423 80 CONTINUE
424 90 CONTINUE
425*
426 ELSE IF( lsame( side, 'R' ) ) THEN
427*
428* Form C * H or C * H**T where C = ( C1 C2 )
429*
430* W := C * V = (C1*V1 + C2*V2) (stored in WORK)
431*
432* W := C2
433*
434 DO 100 j = 1, k
435 CALL dcopy( m, c( 1, n-k+j ), 1, work( 1, j ), 1 )
436 100 CONTINUE
437*
438* W := W * V2
439*
440 CALL dtrmm( 'Right', 'Upper', 'No transpose', 'Unit',
441 $ m,
442 $ k, one, v( n-k+1, 1 ), ldv, work, ldwork )
443 IF( n.GT.k ) THEN
444*
445* W := W + C1 * V1
446*
447 CALL dgemm( 'No transpose', 'No transpose', m, k,
448 $ n-k,
449 $ one, c, ldc, v, ldv, one, work, ldwork )
450 END IF
451*
452* W := W * T or W * T**T
453*
454 CALL dtrmm( 'Right', 'Lower', trans, 'Non-unit', m, k,
455 $ one, t, ldt, work, ldwork )
456*
457* C := C - W * V**T
458*
459 IF( n.GT.k ) THEN
460*
461* C1 := C1 - W * V1**T
462*
463 CALL dgemm( 'No transpose', 'Transpose', m, n-k, k,
464 $ -one, work, ldwork, v, ldv, one, c, ldc )
465 END IF
466*
467* W := W * V2**T
468*
469 CALL dtrmm( 'Right', 'Upper', 'Transpose', 'Unit', m,
470 $ k,
471 $ one, v( n-k+1, 1 ), ldv, work, ldwork )
472*
473* C2 := C2 - W
474*
475 DO 120 j = 1, k
476 DO 110 i = 1, m
477 c( i, n-k+j ) = c( i, n-k+j ) - work( i, j )
478 110 CONTINUE
479 120 CONTINUE
480 END IF
481 END IF
482*
483 ELSE IF( lsame( storev, 'R' ) ) THEN
484*
485 IF( lsame( direct, 'F' ) ) THEN
486*
487* Let V = ( V1 V2 ) (V1: first K columns)
488* where V1 is unit upper triangular.
489*
490 IF( lsame( side, 'L' ) ) THEN
491*
492* Form H * C or H**T * C where C = ( C1 )
493* ( C2 )
494*
495* W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
496*
497* W := C1**T
498*
499 DO 130 j = 1, k
500 CALL dcopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
501 130 CONTINUE
502*
503* W := W * V1**T
504*
505 CALL dtrmm( 'Right', 'Upper', 'Transpose', 'Unit', n,
506 $ k,
507 $ one, v, ldv, work, ldwork )
508 IF( m.GT.k ) THEN
509*
510* W := W + C2**T * V2**T
511*
512 CALL dgemm( 'Transpose', 'Transpose', n, k, m-k,
513 $ one,
514 $ c( k+1, 1 ), ldc, v( 1, k+1 ), ldv, one,
515 $ work, ldwork )
516 END IF
517*
518* W := W * T**T or W * T
519*
520 CALL dtrmm( 'Right', 'Upper', transt, 'Non-unit', n,
521 $ k,
522 $ one, t, ldt, work, ldwork )
523*
524* C := C - V**T * W**T
525*
526 IF( m.GT.k ) THEN
527*
528* C2 := C2 - V2**T * W**T
529*
530 CALL dgemm( 'Transpose', 'Transpose', m-k, n, k,
531 $ -one,
532 $ v( 1, k+1 ), ldv, work, ldwork, one,
533 $ c( k+1, 1 ), ldc )
534 END IF
535*
536* W := W * V1
537*
538 CALL dtrmm( 'Right', 'Upper', 'No transpose', 'Unit',
539 $ n,
540 $ k, one, v, ldv, work, ldwork )
541*
542* C1 := C1 - W**T
543*
544 DO 150 j = 1, k
545 DO 140 i = 1, n
546 c( j, i ) = c( j, i ) - work( i, j )
547 140 CONTINUE
548 150 CONTINUE
549*
550 ELSE IF( lsame( side, 'R' ) ) THEN
551*
552* Form C * H or C * H**T where C = ( C1 C2 )
553*
554* W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK)
555*
556* W := C1
557*
558 DO 160 j = 1, k
559 CALL dcopy( m, c( 1, j ), 1, work( 1, j ), 1 )
560 160 CONTINUE
561*
562* W := W * V1**T
563*
564 CALL dtrmm( 'Right', 'Upper', 'Transpose', 'Unit', m,
565 $ k,
566 $ one, v, ldv, work, ldwork )
567 IF( n.GT.k ) THEN
568*
569* W := W + C2 * V2**T
570*
571 CALL dgemm( 'No transpose', 'Transpose', m, k, n-k,
572 $ one, c( 1, k+1 ), ldc, v( 1, k+1 ), ldv,
573 $ one, work, ldwork )
574 END IF
575*
576* W := W * T or W * T**T
577*
578 CALL dtrmm( 'Right', 'Upper', trans, 'Non-unit', m, k,
579 $ one, t, ldt, work, ldwork )
580*
581* C := C - W * V
582*
583 IF( n.GT.k ) THEN
584*
585* C2 := C2 - W * V2
586*
587 CALL dgemm( 'No transpose', 'No transpose', m, n-k,
588 $ k,
589 $ -one, work, ldwork, v( 1, k+1 ), ldv, one,
590 $ c( 1, k+1 ), ldc )
591 END IF
592*
593* W := W * V1
594*
595 CALL dtrmm( 'Right', 'Upper', 'No transpose', 'Unit',
596 $ m,
597 $ k, one, v, ldv, work, ldwork )
598*
599* C1 := C1 - W
600*
601 DO 180 j = 1, k
602 DO 170 i = 1, m
603 c( i, j ) = c( i, j ) - work( i, j )
604 170 CONTINUE
605 180 CONTINUE
606*
607 END IF
608*
609 ELSE
610*
611* Let V = ( V1 V2 ) (V2: last K columns)
612* where V2 is unit lower triangular.
613*
614 IF( lsame( side, 'L' ) ) THEN
615*
616* Form H * C or H**T * C where C = ( C1 )
617* ( C2 )
618*
619* W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
620*
621* W := C2**T
622*
623 DO 190 j = 1, k
624 CALL dcopy( n, c( m-k+j, 1 ), ldc, work( 1, j ),
625 $ 1 )
626 190 CONTINUE
627*
628* W := W * V2**T
629*
630 CALL dtrmm( 'Right', 'Lower', 'Transpose', 'Unit', n,
631 $ k,
632 $ one, v( 1, m-k+1 ), ldv, work, ldwork )
633 IF( m.GT.k ) THEN
634*
635* W := W + C1**T * V1**T
636*
637 CALL dgemm( 'Transpose', 'Transpose', n, k, m-k,
638 $ one,
639 $ c, ldc, v, ldv, one, work, ldwork )
640 END IF
641*
642* W := W * T**T or W * T
643*
644 CALL dtrmm( 'Right', 'Lower', transt, 'Non-unit', n,
645 $ k,
646 $ one, t, ldt, work, ldwork )
647*
648* C := C - V**T * W**T
649*
650 IF( m.GT.k ) THEN
651*
652* C1 := C1 - V1**T * W**T
653*
654 CALL dgemm( 'Transpose', 'Transpose', m-k, n, k,
655 $ -one,
656 $ v, ldv, work, ldwork, one, c, ldc )
657 END IF
658*
659* W := W * V2
660*
661 CALL dtrmm( 'Right', 'Lower', 'No transpose', 'Unit',
662 $ n,
663 $ k, one, v( 1, m-k+1 ), ldv, work, ldwork )
664*
665* C2 := C2 - W**T
666*
667 DO 210 j = 1, k
668 DO 200 i = 1, n
669 c( m-k+j, i ) = c( m-k+j, i ) - work( i, j )
670 200 CONTINUE
671 210 CONTINUE
672*
673 ELSE IF( lsame( side, 'R' ) ) THEN
674*
675* Form C * H or C * H' where C = ( C1 C2 )
676*
677* W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK)
678*
679* W := C2
680*
681 DO 220 j = 1, k
682 CALL dcopy( m, c( 1, n-k+j ), 1, work( 1, j ), 1 )
683 220 CONTINUE
684*
685* W := W * V2**T
686*
687 CALL dtrmm( 'Right', 'Lower', 'Transpose', 'Unit', m,
688 $ k,
689 $ one, v( 1, n-k+1 ), ldv, work, ldwork )
690 IF( n.GT.k ) THEN
691*
692* W := W + C1 * V1**T
693*
694 CALL dgemm( 'No transpose', 'Transpose', m, k, n-k,
695 $ one, c, ldc, v, ldv, one, work, ldwork )
696 END IF
697*
698* W := W * T or W * T**T
699*
700 CALL dtrmm( 'Right', 'Lower', trans, 'Non-unit', m, k,
701 $ one, t, ldt, work, ldwork )
702*
703* C := C - W * V
704*
705 IF( n.GT.k ) THEN
706*
707* C1 := C1 - W * V1
708*
709 CALL dgemm( 'No transpose', 'No transpose', m, n-k,
710 $ k,
711 $ -one, work, ldwork, v, ldv, one, c, ldc )
712 END IF
713*
714* W := W * V2
715*
716 CALL dtrmm( 'Right', 'Lower', 'No transpose', 'Unit',
717 $ m,
718 $ k, one, v( 1, n-k+1 ), ldv, work, ldwork )
719*
720* C1 := C1 - W
721*
722 DO 240 j = 1, k
723 DO 230 i = 1, m
724 c( i, n-k+j ) = c( i, n-k+j ) - work( i, j )
725 230 CONTINUE
726 240 CONTINUE
727*
728 END IF
729*
730 END IF
731 END IF
732*
733 RETURN
734*
735* End of DLARFB
736*
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dtrmm
subroutine dtrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRMM
Definition dtrmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: