LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ slarfb()

subroutine slarfb ( character  SIDE,
character  TRANS,
character  DIRECT,
character  STOREV,
integer  M,
integer  N,
integer  K,
real, dimension( ldv, * )  V,
integer  LDV,
real, dimension( ldt, * )  T,
integer  LDT,
real, dimension( ldc, * )  C,
integer  LDC,
real, dimension( ldwork, * )  WORK,
integer  LDWORK 
)

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

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

Purpose:
 SLARFB 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 REAL 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 REAL 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 REAL 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 REAL 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 199 of file slarfb.f.

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