1: *> \brief \b ZLARFB_GETT
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZLARFB_GETT + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarfb_gett.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarfb_gett.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarfb_gett.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLARFB_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*16 A( LDA, * ), B( LDB, * ), T( LDT, * ),
31: * $ WORK( LDWORK, * )
32: * ..
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> ZLARFB_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*16 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*16 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*16 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*16 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 complex16OTHERauxiliary
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 ZLARFB_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*16 A( LDA, * ), B( LDB, * ), T( LDT, * ),
404: $ WORK( LDWORK, * )
405: * ..
406: *
407: * =====================================================================
408: *
409: * .. Parameters ..
410: COMPLEX*16 CONE, CZERO
411: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
412: $ CZERO = ( 0.0D+0, 0.0D+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 ZCOPY, ZGEMM, ZTRMM
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 ZCOPY( 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 identy matrix, but unit lower-triangular
456: * V1 stored in A1 (diagonal ones are not stored).
457: *
458: *
459: CALL ZTRMM( '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 ZGEMM( '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 ZTRMM( '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 ZGEMM( '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 ZTRMM( '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 ZCOPY( 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 ZTRMM( '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 ZTRMM( '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 ZTRMM( '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 ZTRMM( '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 ZLARFB_GETT
596: *
597: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlarfb_gett.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack