LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlarfb_gett.f
Go to the documentation of this file.
1*> \brief \b DLARFB_GETT
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLARFB_GETT + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfb_gett.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfb_gett.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfb_gett.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLARFB_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* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
31* $ WORK( LDWORK, * )
32* ..
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DLARFB_GETT applies a real Householder block reflector H from the
40*> left to a real (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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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**T ) * ( A_in ) =
190*> ( B_out ) ( B_in ) ( B_in )
191*> = ( I - ( V1 ) * T * ( V1**T, V2**T ) ) * ( 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**T ) * ( A1_in )
255*> ( B1_out ) ( 0 ) ( 0 )
256*>
257*> ( A2_out ) := H * ( A2_in ) = ( I - V * T * V**T ) * ( 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**T)*A1_in
265*>
266*> B1_out: = - V2*T*(V1**T)*A1_in
267*>
268*> The computation for column block 2, which exists if N > K:
269*>
270*> A2_out: = A2_in - V1*T*( (V1**T)*A2_in + (V2**T)*B2_in )
271*>
272*> B2_out: = B2_in - V2*T*( (V1**T)*A2_in + (V2**T)*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**T)*B2_in )
285*>
286*> B2_out: = B2_in - V2*T*( A2_in + (V2**T)*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**T) * W2 = (unit_lower_tr_of_(A1)**T) * W2
312*> col2_(3) W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * 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**T) * W1 = (unit_lower_tr_of_(A1)**T) * 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**T) * B2 = W2 + (B1**T) * 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**T) * W2
358*> = (unit_lower_tr_of_(A1)**T) * W2
359*> end if
360*> col2_(3) W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * 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**T) * W1
375*> = (unit_lower_tr_of_(A1)**T) * 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 dlarfb_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 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
404 $ work( ldwork, * )
405* ..
406*
407* =====================================================================
408*
409* .. Parameters ..
410 DOUBLE PRECISION ONE, ZERO
411 parameter( one = 1.0d+0, zero = 0.0d+0 )
412* ..
413* .. Local Scalars ..
414 LOGICAL LNOTIDENT
415 INTEGER I, J
416* ..
417* .. EXTERNAL FUNCTIONS ..
418 LOGICAL LSAME
419 EXTERNAL lsame
420* ..
421* .. External Subroutines ..
422 EXTERNAL dcopy, dgemm, dtrmm
423* ..
424* .. Executable Statements ..
425*
426* Quick return if possible
427*
428 IF( m.LT.0 .OR. n.LE.0 .OR. k.EQ.0 .OR. k.GT.n )
429 $ RETURN
430*
431 lnotident = .NOT.lsame( ident, 'I' )
432*
433* ------------------------------------------------------------------
434*
435* First Step. Computation of the Column Block 2:
436*
437* ( A2 ) := H * ( A2 )
438* ( B2 ) ( B2 )
439*
440* ------------------------------------------------------------------
441*
442 IF( n.GT.k ) THEN
443*
444* col2_(1) Compute W2: = A2. Therefore, copy A2 = A(1:K, K+1:N)
445* into W2=WORK(1:K, 1:N-K) column-by-column.
446*
447 DO j = 1, n-k
448 CALL dcopy( k, a( 1, k+j ), 1, work( 1, j ), 1 )
449 END DO
450
451 IF( lnotident ) THEN
452*
453* col2_(2) Compute W2: = (V1**T) * W2 = (A1**T) * W2,
454* V1 is not an identity matrix, but unit lower-triangular
455* V1 stored in A1 (diagonal ones are not stored).
456*
457*
458 CALL dtrmm( 'L', 'L', 'T', 'U', k, n-k, one, a, lda,
459 $ work, ldwork )
460 END IF
461*
462* col2_(3) Compute W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * B2
463* V2 stored in B1.
464*
465 IF( m.GT.0 ) THEN
466 CALL dgemm( 'T', 'N', k, n-k, m, one, b, ldb,
467 $ b( 1, k+1 ), ldb, one, work, ldwork )
468 END IF
469*
470* col2_(4) Compute W2: = T * W2,
471* T is upper-triangular.
472*
473 CALL dtrmm( 'L', 'U', 'N', 'N', k, n-k, one, t, ldt,
474 $ work, ldwork )
475*
476* col2_(5) Compute B2: = B2 - V2 * W2 = B2 - B1 * W2,
477* V2 stored in B1.
478*
479 IF( m.GT.0 ) THEN
480 CALL dgemm( 'N', 'N', m, n-k, k, -one, b, ldb,
481 $ work, ldwork, one, b( 1, k+1 ), ldb )
482 END IF
483*
484 IF( lnotident ) THEN
485*
486* col2_(6) Compute W2: = V1 * W2 = A1 * W2,
487* V1 is not an identity matrix, but unit lower-triangular,
488* V1 stored in A1 (diagonal ones are not stored).
489*
490 CALL dtrmm( 'L', 'L', 'N', 'U', k, n-k, one, a, lda,
491 $ work, ldwork )
492 END IF
493*
494* col2_(7) Compute A2: = A2 - W2 =
495* = A(1:K, K+1:N-K) - WORK(1:K, 1:N-K),
496* column-by-column.
497*
498 DO j = 1, n-k
499 DO i = 1, k
500 a( i, k+j ) = a( i, k+j ) - work( i, j )
501 END DO
502 END DO
503*
504 END IF
505*
506* ------------------------------------------------------------------
507*
508* Second Step. Computation of the Column Block 1:
509*
510* ( A1 ) := H * ( A1 )
511* ( B1 ) ( 0 )
512*
513* ------------------------------------------------------------------
514*
515* col1_(1) Compute W1: = A1. Copy the upper-triangular
516* A1 = A(1:K, 1:K) into the upper-triangular
517* W1 = WORK(1:K, 1:K) column-by-column.
518*
519 DO j = 1, k
520 CALL dcopy( j, a( 1, j ), 1, work( 1, j ), 1 )
521 END DO
522*
523* Set the subdiagonal elements of W1 to zero column-by-column.
524*
525 DO j = 1, k - 1
526 DO i = j + 1, k
527 work( i, j ) = zero
528 END DO
529 END DO
530*
531 IF( lnotident ) THEN
532*
533* col1_(2) Compute W1: = (V1**T) * W1 = (A1**T) * W1,
534* V1 is not an identity matrix, but unit lower-triangular
535* V1 stored in A1 (diagonal ones are not stored),
536* W1 is upper-triangular with zeroes below the diagonal.
537*
538 CALL dtrmm( 'L', 'L', 'T', 'U', k, k, one, a, lda,
539 $ work, ldwork )
540 END IF
541*
542* col1_(3) Compute W1: = T * W1,
543* T is upper-triangular,
544* W1 is upper-triangular with zeroes below the diagonal.
545*
546 CALL dtrmm( 'L', 'U', 'N', 'N', k, k, one, t, ldt,
547 $ work, ldwork )
548*
549* col1_(4) Compute B1: = - V2 * W1 = - B1 * W1,
550* V2 = B1, W1 is upper-triangular with zeroes below the diagonal.
551*
552 IF( m.GT.0 ) THEN
553 CALL dtrmm( 'R', 'U', 'N', 'N', m, k, -one, work, ldwork,
554 $ b, ldb )
555 END IF
556*
557 IF( lnotident ) THEN
558*
559* col1_(5) Compute W1: = V1 * W1 = A1 * W1,
560* V1 is not an identity matrix, but unit lower-triangular
561* V1 stored in A1 (diagonal ones are not stored),
562* W1 is upper-triangular on input with zeroes below the diagonal,
563* and square on output.
564*
565 CALL dtrmm( 'L', 'L', 'N', 'U', k, k, one, a, lda,
566 $ work, ldwork )
567*
568* col1_(6) Compute A1: = A1 - W1 = A(1:K, 1:K) - WORK(1:K, 1:K)
569* column-by-column. A1 is upper-triangular on input.
570* If IDENT, A1 is square on output, and W1 is square,
571* if NOT IDENT, A1 is upper-triangular on output,
572* W1 is upper-triangular.
573*
574* col1_(6)_a Compute elements of A1 below the diagonal.
575*
576 DO j = 1, k - 1
577 DO i = j + 1, k
578 a( i, j ) = - work( i, j )
579 END DO
580 END DO
581*
582 END IF
583*
584* col1_(6)_b Compute elements of A1 on and above the diagonal.
585*
586 DO j = 1, k
587 DO i = 1, j
588 a( i, j ) = a( i, j ) - work( i, j )
589 END DO
590 END DO
591*
592 RETURN
593*
594* End of DLARFB_GETT
595*
596 END
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dlarfb_gett(ident, m, n, k, t, ldt, a, lda, b, ldb, work, ldwork)
DLARFB_GETT
subroutine dtrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRMM
Definition dtrmm.f:177