LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ clarfb()

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