LAPACK 3.12.0
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 elements equal to 1 are not stored; the corresponding
  array elements are modified but restored on exit. The rest of the
  array is not used.

  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 195 of file dlarfb.f.

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