LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
clarfb_gett.f
Go to the documentation of this file.
1*> \brief \b CLARFB_GETT
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CLARFB_GETT + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarfb_gett.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarfb_gett.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarfb_gett.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*>
18* Definition:
19* ===========
20*
21* SUBROUTINE CLARFB_GETT( IDENT, M, N, K, T, LDT, A, LDA, B, LDB,
22* $ WORK, LDWORK )
23* IMPLICIT NONE
24*
25* .. Scalar Arguments ..
26* CHARACTER IDENT
27* INTEGER K, LDA, LDB, LDT, LDWORK, M, N
28* ..
29* .. Array Arguments ..
30* COMPLEX A( LDA, * ), B( LDB, * ), T( LDT, * ),
31* $ WORK( LDWORK, * )
32* ..
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> CLARFB_GETT applies a complex Householder block reflector H from the
40*> left to a complex (K+M)-by-N "triangular-pentagonal" matrix
41*> composed of two block matrices: an upper trapezoidal K-by-N matrix A
42*> stored in the array A, and a rectangular M-by-(N-K) matrix B, stored
43*> in the array B. The block reflector H is stored in a compact
44*> WY-representation, where the elementary reflectors are in the
45*> arrays A, B and T. See Further Details section.
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] IDENT
52*> \verbatim
53*> IDENT is CHARACTER*1
54*> If IDENT = not 'I', or not 'i', then V1 is unit
55*> lower-triangular and stored in the left K-by-K block of
56*> the input matrix A,
57*> If IDENT = 'I' or 'i', then V1 is an identity matrix and
58*> not stored.
59*> See Further Details section.
60*> \endverbatim
61*>
62*> \param[in] M
63*> \verbatim
64*> M is INTEGER
65*> The number of rows of the matrix B.
66*> M >= 0.
67*> \endverbatim
68*>
69*> \param[in] N
70*> \verbatim
71*> N is INTEGER
72*> The number of columns of the matrices A and B.
73*> N >= 0.
74*> \endverbatim
75*>
76*> \param[in] K
77*> \verbatim
78*> K is INTEGER
79*> The number or rows of the matrix A.
80*> K is also order of the matrix T, i.e. the number of
81*> elementary reflectors whose product defines the block
82*> reflector. 0 <= K <= N.
83*> \endverbatim
84*>
85*> \param[in] T
86*> \verbatim
87*> T is COMPLEX array, dimension (LDT,K)
88*> The upper-triangular K-by-K matrix T in the representation
89*> of the block reflector.
90*> \endverbatim
91*>
92*> \param[in] LDT
93*> \verbatim
94*> LDT is INTEGER
95*> The leading dimension of the array T. LDT >= K.
96*> \endverbatim
97*>
98*> \param[in,out] A
99*> \verbatim
100*> A is COMPLEX array, dimension (LDA,N)
101*>
102*> On entry:
103*> a) In the K-by-N upper-trapezoidal part A: input matrix A.
104*> b) In the columns below the diagonal: columns of V1
105*> (ones are not stored on the diagonal).
106*>
107*> On exit:
108*> A is overwritten by rectangular K-by-N product H*A.
109*>
110*> See Further Details section.
111*> \endverbatim
112*>
113*> \param[in] LDA
114*> \verbatim
115*> LDB is INTEGER
116*> The leading dimension of the array A. LDA >= max(1,K).
117*> \endverbatim
118*>
119*> \param[in,out] B
120*> \verbatim
121*> B is COMPLEX array, dimension (LDB,N)
122*>
123*> On entry:
124*> a) In the M-by-(N-K) right block: input matrix B.
125*> b) In the M-by-N left block: columns of V2.
126*>
127*> On exit:
128*> B is overwritten by rectangular M-by-N product H*B.
129*>
130*> See Further Details section.
131*> \endverbatim
132*>
133*> \param[in] LDB
134*> \verbatim
135*> LDB is INTEGER
136*> The leading dimension of the array B. LDB >= max(1,M).
137*> \endverbatim
138*>
139*> \param[out] WORK
140*> \verbatim
141*> WORK is COMPLEX array,
142*> dimension (LDWORK,max(K,N-K))
143*> \endverbatim
144*>
145*> \param[in] LDWORK
146*> \verbatim
147*> LDWORK is INTEGER
148*> The leading dimension of the array WORK. LDWORK>=max(1,K).
149*>
150*> \endverbatim
151*
152* Authors:
153* ========
154*
155*> \author Univ. of Tennessee
156*> \author Univ. of California Berkeley
157*> \author Univ. of Colorado Denver
158*> \author NAG Ltd.
159*
160*> \ingroup larfb_gett
161*
162*> \par Contributors:
163* ==================
164*>
165*> \verbatim
166*>
167*> November 2020, Igor Kozachenko,
168*> Computer Science Division,
169*> University of California, Berkeley
170*>
171*> \endverbatim
172*
173*> \par Further Details:
174* =====================
175*>
176*> \verbatim
177*>
178*> (1) Description of the Algebraic Operation.
179*>
180*> The matrix A is a K-by-N matrix composed of two column block
181*> matrices, A1, which is K-by-K, and A2, which is K-by-(N-K):
182*> A = ( A1, A2 ).
183*> The matrix B is an M-by-N matrix composed of two column block
184*> matrices, B1, which is M-by-K, and B2, which is M-by-(N-K):
185*> B = ( B1, B2 ).
186*>
187*> Perform the operation:
188*>
189*> ( A_out ) := H * ( A_in ) = ( I - V * T * V**H ) * ( A_in ) =
190*> ( B_out ) ( B_in ) ( B_in )
191*> = ( I - ( V1 ) * T * ( V1**H, V2**H ) ) * ( A_in )
192*> ( V2 ) ( B_in )
193*> On input:
194*>
195*> a) ( A_in ) consists of two block columns:
196*> ( B_in )
197*>
198*> ( A_in ) = (( A1_in ) ( A2_in )) = (( A1_in ) ( A2_in ))
199*> ( B_in ) (( B1_in ) ( B2_in )) (( 0 ) ( B2_in )),
200*>
201*> where the column blocks are:
202*>
203*> ( A1_in ) is a K-by-K upper-triangular matrix stored in the
204*> upper triangular part of the array A(1:K,1:K).
205*> ( B1_in ) is an M-by-K rectangular ZERO matrix and not stored.
206*>
207*> ( A2_in ) is a K-by-(N-K) rectangular matrix stored
208*> in the array A(1:K,K+1:N).
209*> ( B2_in ) is an M-by-(N-K) rectangular matrix stored
210*> in the array B(1:M,K+1:N).
211*>
212*> b) V = ( V1 )
213*> ( V2 )
214*>
215*> where:
216*> 1) if IDENT == 'I',V1 is a K-by-K identity matrix, not stored;
217*> 2) if IDENT != 'I',V1 is a K-by-K unit lower-triangular matrix,
218*> stored in the lower-triangular part of the array
219*> A(1:K,1:K) (ones are not stored),
220*> and V2 is an M-by-K rectangular stored the array B(1:M,1:K),
221*> (because on input B1_in is a rectangular zero
222*> matrix that is not stored and the space is
223*> used to store V2).
224*>
225*> c) T is a K-by-K upper-triangular matrix stored
226*> in the array T(1:K,1:K).
227*>
228*> On output:
229*>
230*> a) ( A_out ) consists of two block columns:
231*> ( B_out )
232*>
233*> ( A_out ) = (( A1_out ) ( A2_out ))
234*> ( B_out ) (( B1_out ) ( B2_out )),
235*>
236*> where the column blocks are:
237*>
238*> ( A1_out ) is a K-by-K square matrix, or a K-by-K
239*> upper-triangular matrix, if V1 is an
240*> identity matrix. AiOut is stored in
241*> the array A(1:K,1:K).
242*> ( B1_out ) is an M-by-K rectangular matrix stored
243*> in the array B(1:M,K:N).
244*>
245*> ( A2_out ) is a K-by-(N-K) rectangular matrix stored
246*> in the array A(1:K,K+1:N).
247*> ( B2_out ) is an M-by-(N-K) rectangular matrix stored
248*> in the array B(1:M,K+1:N).
249*>
250*>
251*> The operation above can be represented as the same operation
252*> on each block column:
253*>
254*> ( A1_out ) := H * ( A1_in ) = ( I - V * T * V**H ) * ( A1_in )
255*> ( B1_out ) ( 0 ) ( 0 )
256*>
257*> ( A2_out ) := H * ( A2_in ) = ( I - V * T * V**H ) * ( A2_in )
258*> ( B2_out ) ( B2_in ) ( B2_in )
259*>
260*> If IDENT != 'I':
261*>
262*> The computation for column block 1:
263*>
264*> A1_out: = A1_in - V1*T*(V1**H)*A1_in
265*>
266*> B1_out: = - V2*T*(V1**H)*A1_in
267*>
268*> The computation for column block 2, which exists if N > K:
269*>
270*> A2_out: = A2_in - V1*T*( (V1**H)*A2_in + (V2**H)*B2_in )
271*>
272*> B2_out: = B2_in - V2*T*( (V1**H)*A2_in + (V2**H)*B2_in )
273*>
274*> If IDENT == 'I':
275*>
276*> The operation for column block 1:
277*>
278*> A1_out: = A1_in - V1*T*A1_in
279*>
280*> B1_out: = - V2*T*A1_in
281*>
282*> The computation for column block 2, which exists if N > K:
283*>
284*> A2_out: = A2_in - T*( A2_in + (V2**H)*B2_in )
285*>
286*> B2_out: = B2_in - V2*T*( A2_in + (V2**H)*B2_in )
287*>
288*> (2) Description of the Algorithmic Computation.
289*>
290*> In the first step, we compute column block 2, i.e. A2 and B2.
291*> Here, we need to use the K-by-(N-K) rectangular workspace
292*> matrix W2 that is of the same size as the matrix A2.
293*> W2 is stored in the array WORK(1:K,1:(N-K)).
294*>
295*> In the second step, we compute column block 1, i.e. A1 and B1.
296*> Here, we need to use the K-by-K square workspace matrix W1
297*> that is of the same size as the as the matrix A1.
298*> W1 is stored in the array WORK(1:K,1:K).
299*>
300*> NOTE: Hence, in this routine, we need the workspace array WORK
301*> only of size WORK(1:K,1:max(K,N-K)) so it can hold both W2 from
302*> the first step and W1 from the second step.
303*>
304*> Case (A), when V1 is unit lower-triangular, i.e. IDENT != 'I',
305*> more computations than in the Case (B).
306*>
307*> if( IDENT != 'I' ) then
308*> if ( N > K ) then
309*> (First Step - column block 2)
310*> col2_(1) W2: = A2
311*> col2_(2) W2: = (V1**H) * W2 = (unit_lower_tr_of_(A1)**H) * W2
312*> col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2
313*> col2_(4) W2: = T * W2
314*> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
315*> col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
316*> col2_(7) A2: = A2 - W2
317*> else
318*> (Second Step - column block 1)
319*> col1_(1) W1: = A1
320*> col1_(2) W1: = (V1**H) * W1 = (unit_lower_tr_of_(A1)**H) * W1
321*> col1_(3) W1: = T * W1
322*> col1_(4) B1: = - V2 * W1 = - B1 * W1
323*> col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
324*> col1_(6) square A1: = A1 - W1
325*> end if
326*> end if
327*>
328*> Case (B), when V1 is an identity matrix, i.e. IDENT == 'I',
329*> less computations than in the Case (A)
330*>
331*> if( IDENT == 'I' ) then
332*> if ( N > K ) then
333*> (First Step - column block 2)
334*> col2_(1) W2: = A2
335*> col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2
336*> col2_(4) W2: = T * W2
337*> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
338*> col2_(7) A2: = A2 - W2
339*> else
340*> (Second Step - column block 1)
341*> col1_(1) W1: = A1
342*> col1_(3) W1: = T * W1
343*> col1_(4) B1: = - V2 * W1 = - B1 * W1
344*> col1_(6) upper-triangular_of_(A1): = A1 - W1
345*> end if
346*> end if
347*>
348*> Combine these cases (A) and (B) together, this is the resulting
349*> algorithm:
350*>
351*> if ( N > K ) then
352*>
353*> (First Step - column block 2)
354*>
355*> col2_(1) W2: = A2
356*> if( IDENT != 'I' ) then
357*> col2_(2) W2: = (V1**H) * W2
358*> = (unit_lower_tr_of_(A1)**H) * W2
359*> end if
360*> col2_(3) W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2]
361*> col2_(4) W2: = T * W2
362*> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
363*> if( IDENT != 'I' ) then
364*> col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
365*> end if
366*> col2_(7) A2: = A2 - W2
367*>
368*> else
369*>
370*> (Second Step - column block 1)
371*>
372*> col1_(1) W1: = A1
373*> if( IDENT != 'I' ) then
374*> col1_(2) W1: = (V1**H) * W1
375*> = (unit_lower_tr_of_(A1)**H) * W1
376*> end if
377*> col1_(3) W1: = T * W1
378*> col1_(4) B1: = - V2 * W1 = - B1 * W1
379*> if( IDENT != 'I' ) then
380*> col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
381*> col1_(6_a) below_diag_of_(A1): = - below_diag_of_(W1)
382*> end if
383*> col1_(6_b) up_tr_of_(A1): = up_tr_of_(A1) - up_tr_of_(W1)
384*>
385*> end if
386*>
387*> \endverbatim
388*>
389* =====================================================================
390 SUBROUTINE clarfb_gett( IDENT, M, N, K, T, LDT, A, LDA, B, LDB,
391 $ WORK, LDWORK )
392 IMPLICIT NONE
393*
394* -- LAPACK auxiliary routine --
395* -- LAPACK is a software package provided by Univ. of Tennessee, --
396* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
397*
398* .. Scalar Arguments ..
399 CHARACTER IDENT
400 INTEGER K, LDA, LDB, LDT, LDWORK, M, N
401* ..
402* .. Array Arguments ..
403 COMPLEX A( LDA, * ), B( LDB, * ), T( LDT, * ),
404 $ work( ldwork, * )
405* ..
406*
407* =====================================================================
408*
409* .. Parameters ..
410 COMPLEX CONE, CZERO
411 parameter( cone = ( 1.0e+0, 0.0e+0 ),
412 $ czero = ( 0.0e+0, 0.0e+0 ) )
413* ..
414* .. Local Scalars ..
415 LOGICAL LNOTIDENT
416 INTEGER I, J
417* ..
418* .. EXTERNAL FUNCTIONS ..
419 LOGICAL LSAME
420 EXTERNAL lsame
421* ..
422* .. External Subroutines ..
423 EXTERNAL ccopy, cgemm, ctrmm
424* ..
425* .. Executable Statements ..
426*
427* Quick return if possible
428*
429 IF( m.LT.0 .OR. n.LE.0 .OR. k.EQ.0 .OR. k.GT.n )
430 $ RETURN
431*
432 lnotident = .NOT.lsame( ident, 'I' )
433*
434* ------------------------------------------------------------------
435*
436* First Step. Computation of the Column Block 2:
437*
438* ( A2 ) := H * ( A2 )
439* ( B2 ) ( B2 )
440*
441* ------------------------------------------------------------------
442*
443 IF( n.GT.k ) THEN
444*
445* col2_(1) Compute W2: = A2. Therefore, copy A2 = A(1:K, K+1:N)
446* into W2=WORK(1:K, 1:N-K) column-by-column.
447*
448 DO j = 1, n-k
449 CALL ccopy( k, a( 1, k+j ), 1, work( 1, j ), 1 )
450 END DO
451
452 IF( lnotident ) THEN
453*
454* col2_(2) Compute W2: = (V1**H) * W2 = (A1**H) * W2,
455* V1 is not an identity matrix, but unit lower-triangular
456* V1 stored in A1 (diagonal ones are not stored).
457*
458*
459 CALL ctrmm( 'L', 'L', 'C', 'U', k, n-k, cone, a, lda,
460 $ work, ldwork )
461 END IF
462*
463* col2_(3) Compute W2: = W2 + (V2**H) * B2 = W2 + (B1**H) * B2
464* V2 stored in B1.
465*
466 IF( m.GT.0 ) THEN
467 CALL cgemm( 'C', 'N', k, n-k, m, cone, b, ldb,
468 $ b( 1, k+1 ), ldb, cone, work, ldwork )
469 END IF
470*
471* col2_(4) Compute W2: = T * W2,
472* T is upper-triangular.
473*
474 CALL ctrmm( 'L', 'U', 'N', 'N', k, n-k, cone, t, ldt,
475 $ work, ldwork )
476*
477* col2_(5) Compute B2: = B2 - V2 * W2 = B2 - B1 * W2,
478* V2 stored in B1.
479*
480 IF( m.GT.0 ) THEN
481 CALL cgemm( 'N', 'N', m, n-k, k, -cone, b, ldb,
482 $ work, ldwork, cone, b( 1, k+1 ), ldb )
483 END IF
484*
485 IF( lnotident ) THEN
486*
487* col2_(6) Compute W2: = V1 * W2 = A1 * W2,
488* V1 is not an identity matrix, but unit lower-triangular,
489* V1 stored in A1 (diagonal ones are not stored).
490*
491 CALL ctrmm( 'L', 'L', 'N', 'U', k, n-k, cone, a, lda,
492 $ work, ldwork )
493 END IF
494*
495* col2_(7) Compute A2: = A2 - W2 =
496* = A(1:K, K+1:N-K) - WORK(1:K, 1:N-K),
497* column-by-column.
498*
499 DO j = 1, n-k
500 DO i = 1, k
501 a( i, k+j ) = a( i, k+j ) - work( i, j )
502 END DO
503 END DO
504*
505 END IF
506*
507* ------------------------------------------------------------------
508*
509* Second Step. Computation of the Column Block 1:
510*
511* ( A1 ) := H * ( A1 )
512* ( B1 ) ( 0 )
513*
514* ------------------------------------------------------------------
515*
516* col1_(1) Compute W1: = A1. Copy the upper-triangular
517* A1 = A(1:K, 1:K) into the upper-triangular
518* W1 = WORK(1:K, 1:K) column-by-column.
519*
520 DO j = 1, k
521 CALL ccopy( j, a( 1, j ), 1, work( 1, j ), 1 )
522 END DO
523*
524* Set the subdiagonal elements of W1 to zero column-by-column.
525*
526 DO j = 1, k - 1
527 DO i = j + 1, k
528 work( i, j ) = czero
529 END DO
530 END DO
531*
532 IF( lnotident ) THEN
533*
534* col1_(2) Compute W1: = (V1**H) * W1 = (A1**H) * W1,
535* V1 is not an identity matrix, but unit lower-triangular
536* V1 stored in A1 (diagonal ones are not stored),
537* W1 is upper-triangular with zeroes below the diagonal.
538*
539 CALL ctrmm( 'L', 'L', 'C', 'U', k, k, cone, a, lda,
540 $ work, ldwork )
541 END IF
542*
543* col1_(3) Compute W1: = T * W1,
544* T is upper-triangular,
545* W1 is upper-triangular with zeroes below the diagonal.
546*
547 CALL ctrmm( 'L', 'U', 'N', 'N', k, k, cone, t, ldt,
548 $ work, ldwork )
549*
550* col1_(4) Compute B1: = - V2 * W1 = - B1 * W1,
551* V2 = B1, W1 is upper-triangular with zeroes below the diagonal.
552*
553 IF( m.GT.0 ) THEN
554 CALL ctrmm( 'R', 'U', 'N', 'N', m, k, -cone, work, ldwork,
555 $ b, ldb )
556 END IF
557*
558 IF( lnotident ) THEN
559*
560* col1_(5) Compute W1: = V1 * W1 = A1 * W1,
561* V1 is not an identity matrix, but unit lower-triangular
562* V1 stored in A1 (diagonal ones are not stored),
563* W1 is upper-triangular on input with zeroes below the diagonal,
564* and square on output.
565*
566 CALL ctrmm( 'L', 'L', 'N', 'U', k, k, cone, a, lda,
567 $ work, ldwork )
568*
569* col1_(6) Compute A1: = A1 - W1 = A(1:K, 1:K) - WORK(1:K, 1:K)
570* column-by-column. A1 is upper-triangular on input.
571* If IDENT, A1 is square on output, and W1 is square,
572* if NOT IDENT, A1 is upper-triangular on output,
573* W1 is upper-triangular.
574*
575* col1_(6)_a Compute elements of A1 below the diagonal.
576*
577 DO j = 1, k - 1
578 DO i = j + 1, k
579 a( i, j ) = - work( i, j )
580 END DO
581 END DO
582*
583 END IF
584*
585* col1_(6)_b Compute elements of A1 on and above the diagonal.
586*
587 DO j = 1, k
588 DO i = 1, j
589 a( i, j ) = a( i, j ) - work( i, j )
590 END DO
591 END DO
592*
593 RETURN
594*
595* End of CLARFB_GETT
596*
597 END
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
subroutine clarfb_gett(ident, m, n, k, t, ldt, a, lda, b, ldb, work, ldwork)
CLARFB_GETT
subroutine ctrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRMM
Definition ctrmm.f:177