1 *> \brief \b CLARFB applies a block reflector or its conjugate-transpose to a general rectangular matrix.
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download CLARFB + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarfb.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarfb.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarfb.f">
21 * SUBROUTINE CLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
22 * T, LDT, C, LDC, WORK, LDWORK )
24 * .. Scalar Arguments ..
25 * CHARACTER DIRECT, SIDE, STOREV, TRANS
26 * INTEGER K, LDC, LDT, LDV, LDWORK, M, N
28 * .. Array Arguments ..
29 * COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ),
39 *> CLARFB applies a complex block reflector H or its transpose H**H to a
40 *> complex M-by-N matrix C, from either the left or the right.
48 *> SIDE is CHARACTER*1
49 *> = 'L': apply H or H**H from the Left
50 *> = 'R': apply H or H**H from the Right
55 *> TRANS is CHARACTER*1
56 *> = 'N': apply H (No transpose)
57 *> = 'C': apply H**H (Conjugate transpose)
62 *> DIRECT is CHARACTER*1
63 *> Indicates how H is formed from a product of elementary
65 *> = 'F': H = H(1) H(2) . . . H(k) (Forward)
66 *> = 'B': H = H(k) . . . H(2) H(1) (Backward)
71 *> STOREV is CHARACTER*1
72 *> Indicates how the vectors which define the elementary
73 *> reflectors are stored:
81 *> The number of rows of the matrix C.
87 *> The number of columns of the matrix C.
93 *> The order of the matrix T (= the number of elementary
94 *> reflectors whose product defines the block reflector).
99 *> V is COMPLEX array, dimension
100 *> (LDV,K) if STOREV = 'C'
101 *> (LDV,M) if STOREV = 'R' and SIDE = 'L'
102 *> (LDV,N) if STOREV = 'R' and SIDE = 'R'
103 *> The matrix V. See Further Details.
109 *> The leading dimension of the array V.
110 *> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
111 *> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
112 *> if STOREV = 'R', LDV >= K.
117 *> T is COMPLEX array, dimension (LDT,K)
118 *> The triangular K-by-K matrix T in the representation of the
125 *> The leading dimension of the array T. LDT >= K.
130 *> C is COMPLEX array, dimension (LDC,N)
131 *> On entry, the M-by-N matrix C.
132 *> On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
138 *> The leading dimension of the array C. LDC >= max(1,M).
143 *> WORK is COMPLEX array, dimension (LDWORK,K)
149 *> The leading dimension of the array WORK.
150 *> If SIDE = 'L', LDWORK >= max(1,N);
151 *> if SIDE = 'R', LDWORK >= max(1,M).
157 *> \author Univ. of Tennessee
158 *> \author Univ. of California Berkeley
159 *> \author Univ. of Colorado Denver
164 *> \ingroup complexOTHERauxiliary
166 *> \par Further Details:
167 * =====================
171 *> The shape of the matrix V and the storage of the vectors which define
172 *> the H(i) is best illustrated by the following example with n = 5 and
173 *> k = 3. The elements equal to 1 are not stored; the corresponding
174 *> array elements are modified but restored on exit. The rest of the
175 *> array is not used.
177 *> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R':
179 *> V = ( 1 ) V = ( 1 v1 v1 v1 v1 )
180 *> ( v1 1 ) ( 1 v2 v2 v2 )
181 *> ( v1 v2 1 ) ( 1 v3 v3 )
185 *> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R':
187 *> V = ( v1 v2 v3 ) V = ( v1 v1 1 )
188 *> ( v1 v2 v3 ) ( v2 v2 v2 1 )
189 *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 )
194 * =====================================================================
195 SUBROUTINE CLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
196 $ T, LDT, C, LDC, WORK, LDWORK )
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..--
203 * .. Scalar Arguments ..
204 CHARACTER DIRECT, SIDE, STOREV, TRANS
205 INTEGER K, LDC, LDT, LDV, LDWORK, M, N
207 * .. Array Arguments ..
208 COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ),
212 * =====================================================================
216 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
218 * .. Local Scalars ..
222 * .. External Functions ..
226 * .. External Subroutines ..
227 EXTERNAL CCOPY, CGEMM, CLACGV, CTRMM
229 * .. Intrinsic Functions ..
232 * .. Executable Statements ..
234 * Quick return if possible
236 IF( M.LE.0 .OR. N.LE.0 )
239 IF( LSAME( TRANS, 'N' ) ) THEN
245 IF( LSAME( STOREV, 'C' ) ) THEN
247 IF( LSAME( DIRECT, 'F' ) ) THEN
249 * Let V = ( V1 ) (first K rows)
251 * where V1 is unit lower triangular.
253 IF( LSAME( SIDE, 'L' ) ) THEN
255 * Form H * C or H**H * C where C = ( C1 )
258 * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
263 CALL CCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
264 CALL CLACGV( N, WORK( 1, J ), 1 )
269 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
270 $ K, ONE, V, LDV, WORK, LDWORK )
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 )
280 * W := W * T**H or W * T
282 CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
283 $ ONE, T, LDT, WORK, LDWORK )
289 * C2 := C2 - V2 * W**H
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 )
298 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
299 $ 'Unit', N, K, ONE, V, LDV, WORK, LDWORK )
305 C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) )
309 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
311 * Form C * H or C * H**H where C = ( C1 C2 )
313 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
318 CALL CCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
323 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
324 $ K, ONE, V, LDV, WORK, LDWORK )
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 )
334 * W := W * T or W * T**H
336 CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
337 $ ONE, T, LDT, WORK, LDWORK )
343 * C2 := C2 - W * V2**H
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 )
352 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
353 $ 'Unit', M, K, ONE, V, LDV, WORK, LDWORK )
359 C( I, J ) = C( I, J ) - WORK( I, J )
367 * ( V2 ) (last K rows)
368 * where V2 is unit upper triangular.
370 IF( LSAME( SIDE, 'L' ) ) THEN
372 * Form H * C or H**H * C where C = ( C1 )
375 * W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK)
380 CALL CCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
381 CALL CLACGV( N, WORK( 1, J ), 1 )
386 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
387 $ K, ONE, V( M-K+1, 1 ), LDV, WORK, LDWORK )
390 * W := W + C1**H * V1
392 CALL CGEMM( 'Conjugate transpose', 'No transpose', N,
393 $ K, M-K, ONE, C, LDC, V, LDV, ONE, WORK,
397 * W := W * T**H or W * T
399 CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
400 $ ONE, T, LDT, WORK, LDWORK )
406 * C1 := C1 - V1 * W**H
408 CALL CGEMM( 'No transpose', 'Conjugate transpose',
409 $ M-K, N, K, -ONE, V, LDV, WORK, LDWORK,
415 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
416 $ 'Unit', N, K, ONE, V( M-K+1, 1 ), LDV, WORK,
423 C( M-K+J, I ) = C( M-K+J, I ) -
424 $ CONJG( WORK( I, J ) )
428 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
430 * Form C * H or C * H**H where C = ( C1 C2 )
432 * W := C * V = (C1*V1 + C2*V2) (stored in WORK)
437 CALL CCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
442 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
443 $ K, ONE, V( N-K+1, 1 ), LDV, WORK, LDWORK )
448 CALL CGEMM( 'No transpose', 'No transpose', M, K, N-K,
449 $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
452 * W := W * T or W * T**H
454 CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
455 $ ONE, T, LDT, WORK, LDWORK )
461 * C1 := C1 - W * V1**H
463 CALL CGEMM( 'No transpose', 'Conjugate transpose', M,
464 $ N-K, K, -ONE, WORK, LDWORK, V, LDV, ONE,
470 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
471 $ 'Unit', M, K, ONE, V( N-K+1, 1 ), LDV, WORK,
478 C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )
484 ELSE IF( LSAME( STOREV, 'R' ) ) THEN
486 IF( LSAME( DIRECT, 'F' ) ) THEN
488 * Let V = ( V1 V2 ) (V1: first K columns)
489 * where V1 is unit upper triangular.
491 IF( LSAME( SIDE, 'L' ) ) THEN
493 * Form H * C or H**H * C where C = ( C1 )
496 * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
501 CALL CCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )
502 CALL CLACGV( N, WORK( 1, J ), 1 )
507 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
508 $ 'Unit', N, K, ONE, V, LDV, WORK, LDWORK )
511 * W := W + C2**H * V2**H
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,
519 * W := W * T**H or W * T
521 CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit', N, K,
522 $ ONE, T, LDT, WORK, LDWORK )
524 * C := C - V**H * W**H
528 * C2 := C2 - V2**H * W**H
530 CALL CGEMM( 'Conjugate transpose',
531 $ 'Conjugate transpose', M-K, N, K, -ONE,
532 $ V( 1, K+1 ), LDV, WORK, LDWORK, ONE,
538 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', N,
539 $ K, ONE, V, LDV, WORK, LDWORK )
545 C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) )
549 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
551 * Form C * H or C * H**H where C = ( C1 C2 )
553 * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
558 CALL CCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )
563 CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
564 $ 'Unit', M, K, ONE, V, LDV, WORK, LDWORK )
567 * W := W + C2 * V2**H
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 )
574 * W := W * T or W * T**H
576 CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit', M, K,
577 $ ONE, T, LDT, WORK, LDWORK )
585 CALL CGEMM( 'No transpose', 'No transpose', M, N-K, K,
586 $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, ONE,
592 CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', M,
593 $ K, ONE, V, LDV, WORK, LDWORK )
599 C( I, J ) = C( I, J ) - WORK( I, J )
607 * Let V = ( V1 V2 ) (V2: last K columns)
608 * where V2 is unit lower triangular.
610 IF( LSAME( SIDE, 'L' ) ) THEN
612 * Form H * C or H**H * C where C = ( C1 )
615 * W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
620 CALL CCOPY( N, C( M-K+J, 1 ), LDC, WORK( 1, J ), 1 )
621 CALL CLACGV( N, WORK( 1, J ), 1 )
626 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
627 $ 'Unit', N, K, ONE, V( 1, M-K+1 ), LDV, WORK,
631 * W := W + C1**H * V1**H
633 CALL CGEMM( 'Conjugate transpose',
634 $ 'Conjugate transpose', N, K, M-K, ONE, C,
635 $ LDC, V, LDV, ONE, WORK, LDWORK )
638 * W := W * T**H or W * T
640 CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K,
641 $ ONE, T, LDT, WORK, LDWORK )
643 * C := C - V**H * W**H
647 * C1 := C1 - V1**H * W**H
649 CALL CGEMM( 'Conjugate transpose',
650 $ 'Conjugate transpose', M-K, N, K, -ONE, V,
651 $ LDV, WORK, LDWORK, ONE, C, LDC )
656 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', N,
657 $ K, ONE, V( 1, M-K+1 ), LDV, WORK, LDWORK )
663 C( M-K+J, I ) = C( M-K+J, I ) -
664 $ CONJG( WORK( I, J ) )
668 ELSE IF( LSAME( SIDE, 'R' ) ) THEN
670 * Form C * H or C * H**H where C = ( C1 C2 )
672 * W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK)
677 CALL CCOPY( M, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
682 CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
683 $ 'Unit', M, K, ONE, V( 1, N-K+1 ), LDV, WORK,
687 * W := W + C1 * V1**H
689 CALL CGEMM( 'No transpose', 'Conjugate transpose', M,
690 $ K, N-K, ONE, C, LDC, V, LDV, ONE, WORK,
694 * W := W * T or W * T**H
696 CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K,
697 $ ONE, T, LDT, WORK, LDWORK )
705 CALL CGEMM( 'No transpose', 'No transpose', M, N-K, K,
706 $ -ONE, WORK, LDWORK, V, LDV, ONE, C, LDC )
711 CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', M,
712 $ K, ONE, V( 1, N-K+1 ), LDV, WORK, LDWORK )
718 C( I, N-K+J ) = C( I, N-K+J ) - WORK( I, J )