LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine clarfb ( character  SIDE,
character  TRANS,
character  DIRECT,
character  STOREV,
integer  M,
integer  N,
integer  K,
complex, dimension( ldv, * )  V,
integer  LDV,
complex, dimension( ldt, * )  T,
integer  LDT,
complex, dimension( ldc, * )  C,
integer  LDC,
complex, dimension( ldwork, * )  WORK,
integer  LDWORK 
)

CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix.

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

Purpose: 
 CLARFB applies a complex block reflector H or its transpose H**H to a
 complex M-by-N matrix C, from either the left or the right.
Parameters
[in]SIDE
          SIDE is CHARACTER*1
          = 'L': apply H or H**H from the Left
          = 'R': apply H or H**H from the Right
[in]TRANS
          TRANS is CHARACTER*1
          = 'N': apply H (No transpose)
          = 'C': apply H**H (Conjugate 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 COMPLEX 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 COMPLEX 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 COMPLEX array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
[in]LDC
          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).
[out]WORK
          WORK is COMPLEX 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 clarfb.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  COMPLEX c( ldc, * ), t( ldt, * ), v( ldv, * ),
209  $ work( ldwork, * )
210 * ..
211 *
212 * =====================================================================
213 *
214 * .. Parameters ..
215  COMPLEX one
216  parameter ( one = ( 1.0e+0, 0.0e+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 ccopy, cgemm, clacgv, ctrmm
228 * ..
229 * .. Intrinsic Functions ..
230  INTRINSIC conjg
231 * ..
232 * .. Executable Statements ..
233 *
234 * Quick return if possible
235 *
236  IF( m.LE.0 .OR. n.LE.0 )
237  $ RETURN
238 *
239  IF( lsame( trans, 'N' ) ) THEN
240  transt = 'C'
241  ELSE
242  transt = 'N'
243  END IF
244 *
245  IF( lsame( storev, 'C' ) ) THEN
246 *
247  IF( lsame( direct, 'F' ) ) THEN
248 *
249 * Let V = ( V1 ) (first K rows)
250 * ( V2 )
251 * where V1 is unit lower triangular.
252 *
253  IF( lsame( side, 'L' ) ) THEN
254 *
255 * Form H * C or H**H * C where C = ( C1 )
256 * ( C2 )
257 *
258 * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
259 *
260 * W := C1**H
261 *
262  DO 10 j = 1, k
263  CALL ccopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
264  CALL clacgv( n, work( 1, j ), 1 )
265  10 CONTINUE
266 *
267 * W := W * V1
268 *
269  CALL ctrmm( 'Right', 'Lower', 'No transpose', 'Unit', n,
270  $ k, one, v, ldv, work, ldwork )
271  IF( m.GT.k ) THEN
272 *
273 * W := W + C2**H *V2
274 *
275  CALL cgemm( 'Conjugate transpose', 'No transpose', n,
276  $ k, m-k, one, c( k+1, 1 ), ldc,
277  $ v( k+1, 1 ), ldv, one, work, ldwork )
278  END IF
279 *
280 * W := W * T**H or W * T
281 *
282  CALL ctrmm( 'Right', 'Upper', transt, 'Non-unit', n, k,
283  $ one, t, ldt, work, ldwork )
284 *
285 * C := C - V * W**H
286 *
287  IF( m.GT.k ) THEN
288 *
289 * C2 := C2 - V2 * W**H
290 *
291  CALL cgemm( 'No transpose', 'Conjugate transpose',
292  $ m-k, n, k, -one, v( k+1, 1 ), ldv, work,
293  $ ldwork, one, c( k+1, 1 ), ldc )
294  END IF
295 *
296 * W := W * V1**H
297 *
298  CALL ctrmm( 'Right', 'Lower', 'Conjugate transpose',
299  $ 'Unit', n, k, one, v, ldv, work, ldwork )
300 *
301 * C1 := C1 - W**H
302 *
303  DO 30 j = 1, k
304  DO 20 i = 1, n
305  c( j, i ) = c( j, i ) - conjg( work( i, j ) )
306  20 CONTINUE
307  30 CONTINUE
308 *
309  ELSE IF( lsame( side, 'R' ) ) THEN
310 *
311 * Form C * H or C * H**H where C = ( C1 C2 )
312 *
313 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
314 *
315 * W := C1
316 *
317  DO 40 j = 1, k
318  CALL ccopy( m, c( 1, j ), 1, work( 1, j ), 1 )
319  40 CONTINUE
320 *
321 * W := W * V1
322 *
323  CALL ctrmm( 'Right', 'Lower', 'No transpose', 'Unit', m,
324  $ k, one, v, ldv, work, ldwork )
325  IF( n.GT.k ) THEN
326 *
327 * W := W + C2 * V2
328 *
329  CALL cgemm( 'No transpose', 'No transpose', m, k, n-k,
330  $ one, c( 1, k+1 ), ldc, v( k+1, 1 ), ldv,
331  $ one, work, ldwork )
332  END IF
333 *
334 * W := W * T or W * T**H
335 *
336  CALL ctrmm( 'Right', 'Upper', trans, 'Non-unit', m, k,
337  $ one, t, ldt, work, ldwork )
338 *
339 * C := C - W * V**H
340 *
341  IF( n.GT.k ) THEN
342 *
343 * C2 := C2 - W * V2**H
344 *
345  CALL cgemm( 'No transpose', 'Conjugate transpose', m,
346  $ n-k, k, -one, work, ldwork, v( k+1, 1 ),
347  $ ldv, one, c( 1, k+1 ), ldc )
348  END IF
349 *
350 * W := W * V1**H
351 *
352  CALL ctrmm( 'Right', 'Lower', 'Conjugate transpose',
353  $ 'Unit', m, k, one, v, ldv, work, ldwork )
354 *
355 * C1 := C1 - W
356 *
357  DO 60 j = 1, k
358  DO 50 i = 1, m
359  c( i, j ) = c( i, j ) - work( i, j )
360  50 CONTINUE
361  60 CONTINUE
362  END IF
363 *
364  ELSE
365 *
366 * Let V = ( V1 )
367 * ( V2 ) (last K rows)
368 * where V2 is unit upper triangular.
369 *
370  IF( lsame( side, 'L' ) ) THEN
371 *
372 * Form H * C or H**H * C where C = ( C1 )
373 * ( C2 )
374 *
375 * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
376 *
377 * W := C2**H
378 *
379  DO 70 j = 1, k
380  CALL ccopy( n, c( m-k+j, 1 ), ldc, work( 1, j ), 1 )
381  CALL clacgv( n, work( 1, j ), 1 )
382  70 CONTINUE
383 *
384 * W := W * V2
385 *
386  CALL ctrmm( 'Right', 'Upper', 'No transpose', 'Unit', n,
387  $ k, one, v( m-k+1, 1 ), ldv, work, ldwork )
388  IF( m.GT.k ) THEN
389 *
390 * W := W + C1**H * V1
391 *
392  CALL cgemm( 'Conjugate transpose', 'No transpose', n,
393  $ k, m-k, one, c, ldc, v, ldv, one, work,
394  $ ldwork )
395  END IF
396 *
397 * W := W * T**H or W * T
398 *
399  CALL ctrmm( 'Right', 'Lower', transt, 'Non-unit', n, k,
400  $ one, t, ldt, work, ldwork )
401 *
402 * C := C - V * W**H
403 *
404  IF( m.GT.k ) THEN
405 *
406 * C1 := C1 - V1 * W**H
407 *
408  CALL cgemm( 'No transpose', 'Conjugate transpose',
409  $ m-k, n, k, -one, v, ldv, work, ldwork,
410  $ one, c, ldc )
411  END IF
412 *
413 * W := W * V2**H
414 *
415  CALL ctrmm( 'Right', 'Upper', 'Conjugate transpose',
416  $ 'Unit', n, k, one, v( m-k+1, 1 ), ldv, work,
417  $ ldwork )
418 *
419 * C2 := C2 - W**H
420 *
421  DO 90 j = 1, k
422  DO 80 i = 1, n
423  c( m-k+j, i ) = c( m-k+j, i ) -
424  $ conjg( work( i, j ) )
425  80 CONTINUE
426  90 CONTINUE
427 *
428  ELSE IF( lsame( side, 'R' ) ) THEN
429 *
430 * Form C * H or C * H**H where C = ( C1 C2 )
431 *
432 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
433 *
434 * W := C2
435 *
436  DO 100 j = 1, k
437  CALL ccopy( m, c( 1, n-k+j ), 1, work( 1, j ), 1 )
438  100 CONTINUE
439 *
440 * W := W * V2
441 *
442  CALL ctrmm( 'Right', 'Upper', 'No transpose', 'Unit', m,
443  $ k, one, v( n-k+1, 1 ), ldv, work, ldwork )
444  IF( n.GT.k ) THEN
445 *
446 * W := W + C1 * V1
447 *
448  CALL cgemm( 'No transpose', 'No transpose', m, k, n-k,
449  $ one, c, ldc, v, ldv, one, work, ldwork )
450  END IF
451 *
452 * W := W * T or W * T**H
453 *
454  CALL ctrmm( 'Right', 'Lower', trans, 'Non-unit', m, k,
455  $ one, t, ldt, work, ldwork )
456 *
457 * C := C - W * V**H
458 *
459  IF( n.GT.k ) THEN
460 *
461 * C1 := C1 - W * V1**H
462 *
463  CALL cgemm( 'No transpose', 'Conjugate transpose', m,
464  $ n-k, k, -one, work, ldwork, v, ldv, one,
465  $ c, ldc )
466  END IF
467 *
468 * W := W * V2**H
469 *
470  CALL ctrmm( 'Right', 'Upper', 'Conjugate transpose',
471  $ 'Unit', m, k, one, v( n-k+1, 1 ), ldv, work,
472  $ ldwork )
473 *
474 * C2 := C2 - W
475 *
476  DO 120 j = 1, k
477  DO 110 i = 1, m
478  c( i, n-k+j ) = c( i, n-k+j ) - work( i, j )
479  110 CONTINUE
480  120 CONTINUE
481  END IF
482  END IF
483 *
484  ELSE IF( lsame( storev, 'R' ) ) THEN
485 *
486  IF( lsame( direct, 'F' ) ) THEN
487 *
488 * Let V = ( V1 V2 ) (V1: first K columns)
489 * where V1 is unit upper triangular.
490 *
491  IF( lsame( side, 'L' ) ) THEN
492 *
493 * Form H * C or H**H * C where C = ( C1 )
494 * ( C2 )
495 *
496 * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
497 *
498 * W := C1**H
499 *
500  DO 130 j = 1, k
501  CALL ccopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
502  CALL clacgv( n, work( 1, j ), 1 )
503  130 CONTINUE
504 *
505 * W := W * V1**H
506 *
507  CALL ctrmm( 'Right', 'Upper', 'Conjugate transpose',
508  $ 'Unit', n, k, one, v, ldv, work, ldwork )
509  IF( m.GT.k ) THEN
510 *
511 * W := W + C2**H * V2**H
512 *
513  CALL cgemm( 'Conjugate transpose',
514  $ 'Conjugate transpose', n, k, m-k, one,
515  $ c( k+1, 1 ), ldc, v( 1, k+1 ), ldv, one,
516  $ work, ldwork )
517  END IF
518 *
519 * W := W * T**H or W * T
520 *
521  CALL ctrmm( 'Right', 'Upper', transt, 'Non-unit', n, k,
522  $ one, t, ldt, work, ldwork )
523 *
524 * C := C - V**H * W**H
525 *
526  IF( m.GT.k ) THEN
527 *
528 * C2 := C2 - V2**H * W**H
529 *
530  CALL cgemm( 'Conjugate transpose',
531  $ 'Conjugate transpose', m-k, n, k, -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 ctrmm( 'Right', 'Upper', 'No transpose', 'Unit', n,
539  $ k, one, v, ldv, work, ldwork )
540 *
541 * C1 := C1 - W**H
542 *
543  DO 150 j = 1, k
544  DO 140 i = 1, n
545  c( j, i ) = c( j, i ) - conjg( work( i, j ) )
546  140 CONTINUE
547  150 CONTINUE
548 *
549  ELSE IF( lsame( side, 'R' ) ) THEN
550 *
551 * Form C * H or C * H**H where C = ( C1 C2 )
552 *
553 * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
554 *
555 * W := C1
556 *
557  DO 160 j = 1, k
558  CALL ccopy( m, c( 1, j ), 1, work( 1, j ), 1 )
559  160 CONTINUE
560 *
561 * W := W * V1**H
562 *
563  CALL ctrmm( 'Right', 'Upper', 'Conjugate transpose',
564  $ 'Unit', m, k, one, v, ldv, work, ldwork )
565  IF( n.GT.k ) THEN
566 *
567 * W := W + C2 * V2**H
568 *
569  CALL cgemm( 'No transpose', 'Conjugate transpose', m,
570  $ k, n-k, one, c( 1, k+1 ), ldc,
571  $ v( 1, k+1 ), ldv, one, work, ldwork )
572  END IF
573 *
574 * W := W * T or W * T**H
575 *
576  CALL ctrmm( 'Right', 'Upper', trans, 'Non-unit', m, k,
577  $ one, t, ldt, work, ldwork )
578 *
579 * C := C - W * V
580 *
581  IF( n.GT.k ) THEN
582 *
583 * C2 := C2 - W * V2
584 *
585  CALL cgemm( 'No transpose', 'No transpose', m, n-k, k,
586  $ -one, work, ldwork, v( 1, k+1 ), ldv, one,
587  $ c( 1, k+1 ), ldc )
588  END IF
589 *
590 * W := W * V1
591 *
592  CALL ctrmm( 'Right', 'Upper', 'No transpose', 'Unit', m,
593  $ k, one, v, ldv, work, ldwork )
594 *
595 * C1 := C1 - W
596 *
597  DO 180 j = 1, k
598  DO 170 i = 1, m
599  c( i, j ) = c( i, j ) - work( i, j )
600  170 CONTINUE
601  180 CONTINUE
602 *
603  END IF
604 *
605  ELSE
606 *
607 * Let V = ( V1 V2 ) (V2: last K columns)
608 * where V2 is unit lower triangular.
609 *
610  IF( lsame( side, 'L' ) ) THEN
611 *
612 * Form H * C or H**H * C where C = ( C1 )
613 * ( C2 )
614 *
615 * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
616 *
617 * W := C2**H
618 *
619  DO 190 j = 1, k
620  CALL ccopy( n, c( m-k+j, 1 ), ldc, work( 1, j ), 1 )
621  CALL clacgv( n, work( 1, j ), 1 )
622  190 CONTINUE
623 *
624 * W := W * V2**H
625 *
626  CALL ctrmm( 'Right', 'Lower', 'Conjugate transpose',
627  $ 'Unit', n, k, one, v( 1, m-k+1 ), ldv, work,
628  $ ldwork )
629  IF( m.GT.k ) THEN
630 *
631 * W := W + C1**H * V1**H
632 *
633  CALL cgemm( 'Conjugate transpose',
634  $ 'Conjugate transpose', n, k, m-k, one, c,
635  $ ldc, v, ldv, one, work, ldwork )
636  END IF
637 *
638 * W := W * T**H or W * T
639 *
640  CALL ctrmm( 'Right', 'Lower', transt, 'Non-unit', n, k,
641  $ one, t, ldt, work, ldwork )
642 *
643 * C := C - V**H * W**H
644 *
645  IF( m.GT.k ) THEN
646 *
647 * C1 := C1 - V1**H * W**H
648 *
649  CALL cgemm( 'Conjugate transpose',
650  $ 'Conjugate transpose', m-k, n, k, -one, v,
651  $ ldv, work, ldwork, one, c, ldc )
652  END IF
653 *
654 * W := W * V2
655 *
656  CALL ctrmm( 'Right', 'Lower', 'No transpose', 'Unit', n,
657  $ k, one, v( 1, m-k+1 ), ldv, work, ldwork )
658 *
659 * C2 := C2 - W**H
660 *
661  DO 210 j = 1, k
662  DO 200 i = 1, n
663  c( m-k+j, i ) = c( m-k+j, i ) -
664  $ conjg( work( i, j ) )
665  200 CONTINUE
666  210 CONTINUE
667 *
668  ELSE IF( lsame( side, 'R' ) ) THEN
669 *
670 * Form C * H or C * H**H where C = ( C1 C2 )
671 *
672 * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
673 *
674 * W := C2
675 *
676  DO 220 j = 1, k
677  CALL ccopy( m, c( 1, n-k+j ), 1, work( 1, j ), 1 )
678  220 CONTINUE
679 *
680 * W := W * V2**H
681 *
682  CALL ctrmm( 'Right', 'Lower', 'Conjugate transpose',
683  $ 'Unit', m, k, one, v( 1, n-k+1 ), ldv, work,
684  $ ldwork )
685  IF( n.GT.k ) THEN
686 *
687 * W := W + C1 * V1**H
688 *
689  CALL cgemm( 'No transpose', 'Conjugate transpose', m,
690  $ k, n-k, one, c, ldc, v, ldv, one, work,
691  $ ldwork )
692  END IF
693 *
694 * W := W * T or W * T**H
695 *
696  CALL ctrmm( 'Right', 'Lower', trans, 'Non-unit', m, k,
697  $ one, t, ldt, work, ldwork )
698 *
699 * C := C - W * V
700 *
701  IF( n.GT.k ) THEN
702 *
703 * C1 := C1 - W * V1
704 *
705  CALL cgemm( 'No transpose', 'No transpose', m, n-k, k,
706  $ -one, work, ldwork, v, ldv, one, c, ldc )
707  END IF
708 *
709 * W := W * V2
710 *
711  CALL ctrmm( 'Right', 'Lower', 'No transpose', 'Unit', m,
712  $ k, one, v( 1, n-k+1 ), ldv, work, ldwork )
713 *
714 * C1 := C1 - W
715 *
716  DO 240 j = 1, k
717  DO 230 i = 1, m
718  c( i, n-k+j ) = c( i, n-k+j ) - work( i, j )
719  230 CONTINUE
720  240 CONTINUE
721 *
722  END IF
723 *
724  END IF
725  END IF
726 *
727  RETURN
728 *
729 * End of CLARFB
730 *
ccopy
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:52
clacgv
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:76
ctrmm
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:179
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.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: