LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
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: