LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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).
[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.
Date
June 2013
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 197 of file dlarfb.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: