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

◆ zlarfb()

subroutine zlarfb ( character side,
character trans,
character direct,
character storev,
integer m,
integer n,
integer k,
complex*16, dimension( ldv, * ) v,
integer ldv,
complex*16, dimension( ldt, * ) t,
integer ldt,
complex*16, dimension( ldc, * ) c,
integer ldc,
complex*16, dimension( ldwork, * ) work,
integer ldwork )

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

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

Purpose:
!>
!> ZLARFB 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*16 array, dimension
!>                                (LDV,K) if STOREV = 'C'
!>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
!>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
!>          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*16 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*16 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*16 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 triangular part of V (including its diagonal) is not
!>  referenced.
!>
!>  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 192 of file zlarfb.f.

195*
196* -- LAPACK auxiliary routine --
197* -- LAPACK is a software package provided by Univ. of Tennessee, --
198* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
199*
200* .. Scalar Arguments ..
201 CHARACTER DIRECT, SIDE, STOREV, TRANS
202 INTEGER K, LDC, LDT, LDV, LDWORK, M, N
203* ..
204* .. Array Arguments ..
205 COMPLEX*16 C( LDC, * ), T( LDT, * ), V( LDV, * ),
206 $ WORK( LDWORK, * )
207* ..
208*
209* =====================================================================
210*
211* .. Parameters ..
212 COMPLEX*16 ONE
213 parameter( one = ( 1.0d+0, 0.0d+0 ) )
214* ..
215* .. Local Scalars ..
216 CHARACTER TRANST
217 INTEGER I, J
218* ..
219* .. External Functions ..
220 LOGICAL LSAME
221 EXTERNAL lsame
222* ..
223* .. External Subroutines ..
224 EXTERNAL zcopy, zgemm, zlacgv, ztrmm
225* ..
226* .. Intrinsic Functions ..
227 INTRINSIC dconjg
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 = 'C'
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**H * C where C = ( C1 )
253* ( C2 )
254*
255* W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
256*
257* W := C1**H
258*
259 DO 10 j = 1, k
260 CALL zcopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
261 CALL zlacgv( n, work( 1, j ), 1 )
262 10 CONTINUE
263*
264* W := W * V1
265*
266 CALL ztrmm( 'Right', 'Lower', 'No transpose', 'Unit',
267 $ n,
268 $ k, one, v, ldv, work, ldwork )
269 IF( m.GT.k ) THEN
270*
271* W := W + C2**H * V2
272*
273 CALL zgemm( 'Conjugate transpose', 'No transpose',
274 $ 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 ztrmm( 'Right', 'Upper', transt, 'Non-unit', n,
282 $ 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 zgemm( '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 ztrmm( '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 ) - dconjg( 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 zcopy( m, c( 1, j ), 1, work( 1, j ), 1 )
319 40 CONTINUE
320*
321* W := W * V1
322*
323 CALL ztrmm( 'Right', 'Lower', 'No transpose', 'Unit',
324 $ m,
325 $ k, one, v, ldv, work, ldwork )
326 IF( n.GT.k ) THEN
327*
328* W := W + C2 * V2
329*
330 CALL zgemm( 'No transpose', 'No transpose', m, k,
331 $ n-k,
332 $ one, c( 1, k+1 ), ldc, v( k+1, 1 ), ldv,
333 $ one, work, ldwork )
334 END IF
335*
336* W := W * T or W * T**H
337*
338 CALL ztrmm( 'Right', 'Upper', trans, 'Non-unit', m, k,
339 $ one, t, ldt, work, ldwork )
340*
341* C := C - W * V**H
342*
343 IF( n.GT.k ) THEN
344*
345* C2 := C2 - W * V2**H
346*
347 CALL zgemm( 'No transpose', 'Conjugate transpose',
348 $ m,
349 $ n-k, k, -one, work, ldwork, v( k+1, 1 ),
350 $ ldv, one, c( 1, k+1 ), ldc )
351 END IF
352*
353* W := W * V1**H
354*
355 CALL ztrmm( 'Right', 'Lower', 'Conjugate transpose',
356 $ 'Unit', m, k, one, v, ldv, work, ldwork )
357*
358* C1 := C1 - W
359*
360 DO 60 j = 1, k
361 DO 50 i = 1, m
362 c( i, j ) = c( i, j ) - work( i, j )
363 50 CONTINUE
364 60 CONTINUE
365 END IF
366*
367 ELSE
368*
369* Let V = ( V1 )
370* ( V2 ) (last K rows)
371* where V2 is unit upper triangular.
372*
373 IF( lsame( side, 'L' ) ) THEN
374*
375* Form H * C or H**H * C where C = ( C1 )
376* ( C2 )
377*
378* W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
379*
380* W := C2**H
381*
382 DO 70 j = 1, k
383 CALL zcopy( n, c( m-k+j, 1 ), ldc, work( 1, j ),
384 $ 1 )
385 CALL zlacgv( n, work( 1, j ), 1 )
386 70 CONTINUE
387*
388* W := W * V2
389*
390 CALL ztrmm( 'Right', 'Upper', 'No transpose', 'Unit',
391 $ n,
392 $ k, one, v( m-k+1, 1 ), ldv, work, ldwork )
393 IF( m.GT.k ) THEN
394*
395* W := W + C1**H * V1
396*
397 CALL zgemm( 'Conjugate transpose', 'No transpose',
398 $ n,
399 $ k, m-k, one, c, ldc, v, ldv, one, work,
400 $ ldwork )
401 END IF
402*
403* W := W * T**H or W * T
404*
405 CALL ztrmm( 'Right', 'Lower', transt, 'Non-unit', n,
406 $ k,
407 $ one, t, ldt, work, ldwork )
408*
409* C := C - V * W**H
410*
411 IF( m.GT.k ) THEN
412*
413* C1 := C1 - V1 * W**H
414*
415 CALL zgemm( 'No transpose', 'Conjugate transpose',
416 $ m-k, n, k, -one, v, ldv, work, ldwork,
417 $ one, c, ldc )
418 END IF
419*
420* W := W * V2**H
421*
422 CALL ztrmm( 'Right', 'Upper', 'Conjugate transpose',
423 $ 'Unit', n, k, one, v( m-k+1, 1 ), ldv, work,
424 $ ldwork )
425*
426* C2 := C2 - W**H
427*
428 DO 90 j = 1, k
429 DO 80 i = 1, n
430 c( m-k+j, i ) = c( m-k+j, i ) -
431 $ dconjg( work( i, j ) )
432 80 CONTINUE
433 90 CONTINUE
434*
435 ELSE IF( lsame( side, 'R' ) ) THEN
436*
437* Form C * H or C * H**H where C = ( C1 C2 )
438*
439* W := C * V = (C1*V1 + C2*V2) (stored in WORK)
440*
441* W := C2
442*
443 DO 100 j = 1, k
444 CALL zcopy( m, c( 1, n-k+j ), 1, work( 1, j ), 1 )
445 100 CONTINUE
446*
447* W := W * V2
448*
449 CALL ztrmm( 'Right', 'Upper', 'No transpose', 'Unit',
450 $ m,
451 $ k, one, v( n-k+1, 1 ), ldv, work, ldwork )
452 IF( n.GT.k ) THEN
453*
454* W := W + C1 * V1
455*
456 CALL zgemm( 'No transpose', 'No transpose', m, k,
457 $ n-k,
458 $ one, c, ldc, v, ldv, one, work, ldwork )
459 END IF
460*
461* W := W * T or W * T**H
462*
463 CALL ztrmm( 'Right', 'Lower', trans, 'Non-unit', m, k,
464 $ one, t, ldt, work, ldwork )
465*
466* C := C - W * V**H
467*
468 IF( n.GT.k ) THEN
469*
470* C1 := C1 - W * V1**H
471*
472 CALL zgemm( 'No transpose', 'Conjugate transpose',
473 $ m,
474 $ n-k, k, -one, work, ldwork, v, ldv, one,
475 $ c, ldc )
476 END IF
477*
478* W := W * V2**H
479*
480 CALL ztrmm( 'Right', 'Upper', 'Conjugate transpose',
481 $ 'Unit', m, k, one, v( n-k+1, 1 ), ldv, work,
482 $ ldwork )
483*
484* C2 := C2 - W
485*
486 DO 120 j = 1, k
487 DO 110 i = 1, m
488 c( i, n-k+j ) = c( i, n-k+j ) - work( i, j )
489 110 CONTINUE
490 120 CONTINUE
491 END IF
492 END IF
493*
494 ELSE IF( lsame( storev, 'R' ) ) THEN
495*
496 IF( lsame( direct, 'F' ) ) THEN
497*
498* Let V = ( V1 V2 ) (V1: first K columns)
499* where V1 is unit upper triangular.
500*
501 IF( lsame( side, 'L' ) ) THEN
502*
503* Form H * C or H**H * C where C = ( C1 )
504* ( C2 )
505*
506* W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
507*
508* W := C1**H
509*
510 DO 130 j = 1, k
511 CALL zcopy( n, c( j, 1 ), ldc, work( 1, j ), 1 )
512 CALL zlacgv( n, work( 1, j ), 1 )
513 130 CONTINUE
514*
515* W := W * V1**H
516*
517 CALL ztrmm( 'Right', 'Upper', 'Conjugate transpose',
518 $ 'Unit', n, k, one, v, ldv, work, ldwork )
519 IF( m.GT.k ) THEN
520*
521* W := W + C2**H * V2**H
522*
523 CALL zgemm( 'Conjugate transpose',
524 $ 'Conjugate transpose', n, k, m-k, one,
525 $ c( k+1, 1 ), ldc, v( 1, k+1 ), ldv, one,
526 $ work, ldwork )
527 END IF
528*
529* W := W * T**H or W * T
530*
531 CALL ztrmm( 'Right', 'Upper', transt, 'Non-unit', n,
532 $ k,
533 $ one, t, ldt, work, ldwork )
534*
535* C := C - V**H * W**H
536*
537 IF( m.GT.k ) THEN
538*
539* C2 := C2 - V2**H * W**H
540*
541 CALL zgemm( 'Conjugate transpose',
542 $ 'Conjugate transpose', m-k, n, k, -one,
543 $ v( 1, k+1 ), ldv, work, ldwork, one,
544 $ c( k+1, 1 ), ldc )
545 END IF
546*
547* W := W * V1
548*
549 CALL ztrmm( 'Right', 'Upper', 'No transpose', 'Unit',
550 $ n,
551 $ k, one, v, ldv, work, ldwork )
552*
553* C1 := C1 - W**H
554*
555 DO 150 j = 1, k
556 DO 140 i = 1, n
557 c( j, i ) = c( j, i ) - dconjg( work( i, j ) )
558 140 CONTINUE
559 150 CONTINUE
560*
561 ELSE IF( lsame( side, 'R' ) ) THEN
562*
563* Form C * H or C * H**H where C = ( C1 C2 )
564*
565* W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
566*
567* W := C1
568*
569 DO 160 j = 1, k
570 CALL zcopy( m, c( 1, j ), 1, work( 1, j ), 1 )
571 160 CONTINUE
572*
573* W := W * V1**H
574*
575 CALL ztrmm( 'Right', 'Upper', 'Conjugate transpose',
576 $ 'Unit', m, k, one, v, ldv, work, ldwork )
577 IF( n.GT.k ) THEN
578*
579* W := W + C2 * V2**H
580*
581 CALL zgemm( 'No transpose', 'Conjugate transpose',
582 $ m,
583 $ k, n-k, one, c( 1, k+1 ), ldc,
584 $ v( 1, k+1 ), ldv, one, work, ldwork )
585 END IF
586*
587* W := W * T or W * T**H
588*
589 CALL ztrmm( 'Right', 'Upper', trans, 'Non-unit', m, k,
590 $ one, t, ldt, work, ldwork )
591*
592* C := C - W * V
593*
594 IF( n.GT.k ) THEN
595*
596* C2 := C2 - W * V2
597*
598 CALL zgemm( 'No transpose', 'No transpose', m, n-k,
599 $ k,
600 $ -one, work, ldwork, v( 1, k+1 ), ldv, one,
601 $ c( 1, k+1 ), ldc )
602 END IF
603*
604* W := W * V1
605*
606 CALL ztrmm( 'Right', 'Upper', 'No transpose', 'Unit',
607 $ m,
608 $ k, one, v, ldv, work, ldwork )
609*
610* C1 := C1 - W
611*
612 DO 180 j = 1, k
613 DO 170 i = 1, m
614 c( i, j ) = c( i, j ) - work( i, j )
615 170 CONTINUE
616 180 CONTINUE
617*
618 END IF
619*
620 ELSE
621*
622* Let V = ( V1 V2 ) (V2: last K columns)
623* where V2 is unit lower triangular.
624*
625 IF( lsame( side, 'L' ) ) THEN
626*
627* Form H * C or H**H * C where C = ( C1 )
628* ( C2 )
629*
630* W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
631*
632* W := C2**H
633*
634 DO 190 j = 1, k
635 CALL zcopy( n, c( m-k+j, 1 ), ldc, work( 1, j ),
636 $ 1 )
637 CALL zlacgv( n, work( 1, j ), 1 )
638 190 CONTINUE
639*
640* W := W * V2**H
641*
642 CALL ztrmm( 'Right', 'Lower', 'Conjugate transpose',
643 $ 'Unit', n, k, one, v( 1, m-k+1 ), ldv, work,
644 $ ldwork )
645 IF( m.GT.k ) THEN
646*
647* W := W + C1**H * V1**H
648*
649 CALL zgemm( 'Conjugate transpose',
650 $ 'Conjugate transpose', n, k, m-k, one, c,
651 $ ldc, v, ldv, one, work, ldwork )
652 END IF
653*
654* W := W * T**H or W * T
655*
656 CALL ztrmm( 'Right', 'Lower', transt, 'Non-unit', n,
657 $ k,
658 $ one, t, ldt, work, ldwork )
659*
660* C := C - V**H * W**H
661*
662 IF( m.GT.k ) THEN
663*
664* C1 := C1 - V1**H * W**H
665*
666 CALL zgemm( 'Conjugate transpose',
667 $ 'Conjugate transpose', m-k, n, k, -one, v,
668 $ ldv, work, ldwork, one, c, ldc )
669 END IF
670*
671* W := W * V2
672*
673 CALL ztrmm( 'Right', 'Lower', 'No transpose', 'Unit',
674 $ n,
675 $ k, one, v( 1, m-k+1 ), ldv, work, ldwork )
676*
677* C2 := C2 - W**H
678*
679 DO 210 j = 1, k
680 DO 200 i = 1, n
681 c( m-k+j, i ) = c( m-k+j, i ) -
682 $ dconjg( work( i, j ) )
683 200 CONTINUE
684 210 CONTINUE
685*
686 ELSE IF( lsame( side, 'R' ) ) THEN
687*
688* Form C * H or C * H**H where C = ( C1 C2 )
689*
690* W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
691*
692* W := C2
693*
694 DO 220 j = 1, k
695 CALL zcopy( m, c( 1, n-k+j ), 1, work( 1, j ), 1 )
696 220 CONTINUE
697*
698* W := W * V2**H
699*
700 CALL ztrmm( 'Right', 'Lower', 'Conjugate transpose',
701 $ 'Unit', m, k, one, v( 1, n-k+1 ), ldv, work,
702 $ ldwork )
703 IF( n.GT.k ) THEN
704*
705* W := W + C1 * V1**H
706*
707 CALL zgemm( 'No transpose', 'Conjugate transpose',
708 $ m,
709 $ k, n-k, one, c, ldc, v, ldv, one, work,
710 $ ldwork )
711 END IF
712*
713* W := W * T or W * T**H
714*
715 CALL ztrmm( 'Right', 'Lower', trans, 'Non-unit', m, k,
716 $ one, t, ldt, work, ldwork )
717*
718* C := C - W * V
719*
720 IF( n.GT.k ) THEN
721*
722* C1 := C1 - W * V1
723*
724 CALL zgemm( 'No transpose', 'No transpose', m, n-k,
725 $ k,
726 $ -one, work, ldwork, v, ldv, one, c, ldc )
727 END IF
728*
729* W := W * V2
730*
731 CALL ztrmm( 'Right', 'Lower', 'No transpose', 'Unit',
732 $ m,
733 $ k, one, v( 1, n-k+1 ), ldv, work, ldwork )
734*
735* C1 := C1 - W
736*
737 DO 240 j = 1, k
738 DO 230 i = 1, m
739 c( i, n-k+j ) = c( i, n-k+j ) - work( i, j )
740 230 CONTINUE
741 240 CONTINUE
742*
743 END IF
744*
745 END IF
746 END IF
747*
748 RETURN
749*
750* End of ZLARFB
751*
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
zlacgv
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
Definition zlacgv.f:72
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
ztrmm
subroutine ztrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRMM
Definition ztrmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: