Annotation of rpl/lapack/lapack/ztprfb.f, revision 1.9
1.3 bertrand 1: *> \brief \b ZTPRFB applies a real or complex "triangular-pentagonal" blocked reflector to a real or complex matrix, which is composed of two blocks.
1.1 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.7 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.1 bertrand 7: *
8: *> \htmlonly
1.7 bertrand 9: *> Download ZTPRFB + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztprfb.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztprfb.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztprfb.f">
1.1 bertrand 15: *> [TXT]</a>
1.7 bertrand 16: *> \endhtmlonly
1.1 bertrand 17: *
18: * Definition:
19: * ===========
20: *
1.7 bertrand 21: * SUBROUTINE ZTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
1.1 bertrand 22: * V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
1.7 bertrand 23: *
1.1 bertrand 24: * .. Scalar Arguments ..
25: * CHARACTER DIRECT, SIDE, STOREV, TRANS
26: * INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
27: * ..
28: * .. Array Arguments ..
1.7 bertrand 29: * COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ),
1.1 bertrand 30: * $ V( LDV, * ), WORK( LDWORK, * )
31: * ..
1.7 bertrand 32: *
1.1 bertrand 33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
1.7 bertrand 39: *> ZTPRFB applies a complex "triangular-pentagonal" block reflector H or its
40: *> conjugate transpose H**H to a complex matrix C, which is composed of two
1.1 bertrand 41: *> blocks A and B, either from the left or right.
1.7 bertrand 42: *>
1.1 bertrand 43: *> \endverbatim
44: *
45: * Arguments:
46: * ==========
47: *
48: *> \param[in] SIDE
49: *> \verbatim
50: *> SIDE is CHARACTER*1
51: *> = 'L': apply H or H**H from the Left
52: *> = 'R': apply H or H**H from the Right
53: *> \endverbatim
54: *>
55: *> \param[in] TRANS
56: *> \verbatim
57: *> TRANS is CHARACTER*1
58: *> = 'N': apply H (No transpose)
59: *> = 'C': apply H**H (Conjugate transpose)
60: *> \endverbatim
61: *>
62: *> \param[in] DIRECT
63: *> \verbatim
64: *> DIRECT is CHARACTER*1
65: *> Indicates how H is formed from a product of elementary
66: *> reflectors
67: *> = 'F': H = H(1) H(2) . . . H(k) (Forward)
68: *> = 'B': H = H(k) . . . H(2) H(1) (Backward)
69: *> \endverbatim
70: *>
71: *> \param[in] STOREV
72: *> \verbatim
73: *> STOREV is CHARACTER*1
74: *> Indicates how the vectors which define the elementary
75: *> reflectors are stored:
76: *> = 'C': Columns
77: *> = 'R': Rows
78: *> \endverbatim
79: *>
80: *> \param[in] M
81: *> \verbatim
82: *> M is INTEGER
1.7 bertrand 83: *> The number of rows of the matrix B.
1.1 bertrand 84: *> M >= 0.
85: *> \endverbatim
86: *>
87: *> \param[in] N
88: *> \verbatim
89: *> N is INTEGER
1.7 bertrand 90: *> The number of columns of the matrix B.
1.1 bertrand 91: *> N >= 0.
92: *> \endverbatim
93: *>
94: *> \param[in] K
95: *> \verbatim
96: *> K is INTEGER
97: *> The order of the matrix T, i.e. the number of elementary
1.7 bertrand 98: *> reflectors whose product defines the block reflector.
1.1 bertrand 99: *> K >= 0.
100: *> \endverbatim
101: *>
102: *> \param[in] L
103: *> \verbatim
104: *> L is INTEGER
1.7 bertrand 105: *> The order of the trapezoidal part of V.
1.1 bertrand 106: *> K >= L >= 0. See Further Details.
107: *> \endverbatim
108: *>
109: *> \param[in] V
110: *> \verbatim
111: *> V is COMPLEX*16 array, dimension
112: *> (LDV,K) if STOREV = 'C'
113: *> (LDV,M) if STOREV = 'R' and SIDE = 'L'
114: *> (LDV,N) if STOREV = 'R' and SIDE = 'R'
115: *> The pentagonal matrix V, which contains the elementary reflectors
116: *> H(1), H(2), ..., H(K). See Further Details.
117: *> \endverbatim
118: *>
119: *> \param[in] LDV
120: *> \verbatim
121: *> LDV is INTEGER
122: *> The leading dimension of the array V.
123: *> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
124: *> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
125: *> if STOREV = 'R', LDV >= K.
126: *> \endverbatim
127: *>
128: *> \param[in] T
129: *> \verbatim
130: *> T is COMPLEX*16 array, dimension (LDT,K)
131: *> The triangular K-by-K matrix T in the representation of the
1.7 bertrand 132: *> block reflector.
1.1 bertrand 133: *> \endverbatim
134: *>
135: *> \param[in] LDT
136: *> \verbatim
137: *> LDT is INTEGER
1.7 bertrand 138: *> The leading dimension of the array T.
1.1 bertrand 139: *> LDT >= K.
140: *> \endverbatim
141: *>
142: *> \param[in,out] A
143: *> \verbatim
144: *> A is COMPLEX*16 array, dimension
145: *> (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'
146: *> On entry, the K-by-N or M-by-K matrix A.
1.7 bertrand 147: *> On exit, A is overwritten by the corresponding block of
148: *> H*C or H**H*C or C*H or C*H**H. See Further Details.
1.1 bertrand 149: *> \endverbatim
150: *>
151: *> \param[in] LDA
152: *> \verbatim
153: *> LDA is INTEGER
1.7 bertrand 154: *> The leading dimension of the array A.
1.1 bertrand 155: *> If SIDE = 'L', LDC >= max(1,K);
1.7 bertrand 156: *> If SIDE = 'R', LDC >= max(1,M).
1.1 bertrand 157: *> \endverbatim
158: *>
159: *> \param[in,out] B
160: *> \verbatim
161: *> B is COMPLEX*16 array, dimension (LDB,N)
162: *> On entry, the M-by-N matrix B.
163: *> On exit, B is overwritten by the corresponding block of
164: *> H*C or H**H*C or C*H or C*H**H. See Further Details.
165: *> \endverbatim
166: *>
167: *> \param[in] LDB
168: *> \verbatim
169: *> LDB is INTEGER
1.7 bertrand 170: *> The leading dimension of the array B.
1.1 bertrand 171: *> LDB >= max(1,M).
172: *> \endverbatim
173: *>
174: *> \param[out] WORK
175: *> \verbatim
176: *> WORK is COMPLEX*16 array, dimension
177: *> (LDWORK,N) if SIDE = 'L',
178: *> (LDWORK,K) if SIDE = 'R'.
179: *> \endverbatim
180: *>
181: *> \param[in] LDWORK
182: *> \verbatim
183: *> LDWORK is INTEGER
184: *> The leading dimension of the array WORK.
1.7 bertrand 185: *> If SIDE = 'L', LDWORK >= K;
1.1 bertrand 186: *> if SIDE = 'R', LDWORK >= M.
187: *> \endverbatim
188: *
189: * Authors:
190: * ========
191: *
1.7 bertrand 192: *> \author Univ. of Tennessee
193: *> \author Univ. of California Berkeley
194: *> \author Univ. of Colorado Denver
195: *> \author NAG Ltd.
1.1 bertrand 196: *
1.7 bertrand 197: *> \date December 2016
1.1 bertrand 198: *
199: *> \ingroup complex16OTHERauxiliary
200: *
201: *> \par Further Details:
202: * =====================
203: *>
204: *> \verbatim
205: *>
206: *> The matrix C is a composite matrix formed from blocks A and B.
1.7 bertrand 207: *> The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K,
1.1 bertrand 208: *> and if SIDE = 'L', A is of size K-by-N.
209: *>
210: *> If SIDE = 'R' and DIRECT = 'F', C = [A B].
211: *>
1.7 bertrand 212: *> If SIDE = 'L' and DIRECT = 'F', C = [A]
1.1 bertrand 213: *> [B].
214: *>
215: *> If SIDE = 'R' and DIRECT = 'B', C = [B A].
216: *>
217: *> If SIDE = 'L' and DIRECT = 'B', C = [B]
1.7 bertrand 218: *> [A].
1.1 bertrand 219: *>
1.7 bertrand 220: *> The pentagonal matrix V is composed of a rectangular block V1 and a
221: *> trapezoidal block V2. The size of the trapezoidal block is determined by
1.1 bertrand 222: *> the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular;
223: *> if L=0, there is no trapezoidal block, thus V = V1 is rectangular.
224: *>
225: *> If DIRECT = 'F' and STOREV = 'C': V = [V1]
226: *> [V2]
227: *> - V2 is upper trapezoidal (first L rows of K-by-K upper triangular)
228: *>
229: *> If DIRECT = 'F' and STOREV = 'R': V = [V1 V2]
230: *>
231: *> - V2 is lower trapezoidal (first L columns of K-by-K lower triangular)
232: *>
233: *> If DIRECT = 'B' and STOREV = 'C': V = [V2]
234: *> [V1]
235: *> - V2 is lower trapezoidal (last L rows of K-by-K lower triangular)
236: *>
237: *> If DIRECT = 'B' and STOREV = 'R': V = [V2 V1]
1.7 bertrand 238: *>
1.1 bertrand 239: *> - V2 is upper trapezoidal (last L columns of K-by-K upper triangular)
240: *>
241: *> If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K.
242: *>
243: *> If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K.
244: *>
245: *> If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L.
246: *>
247: *> If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.
248: *> \endverbatim
249: *>
250: * =====================================================================
1.7 bertrand 251: SUBROUTINE ZTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
1.1 bertrand 252: $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
253: *
1.7 bertrand 254: * -- LAPACK auxiliary routine (version 3.7.0) --
1.1 bertrand 255: * -- LAPACK is a software package provided by Univ. of Tennessee, --
256: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.7 bertrand 257: * December 2016
1.1 bertrand 258: *
259: * .. Scalar Arguments ..
260: CHARACTER DIRECT, SIDE, STOREV, TRANS
261: INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
262: * ..
263: * .. Array Arguments ..
1.7 bertrand 264: COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ),
1.1 bertrand 265: $ V( LDV, * ), WORK( LDWORK, * )
266: * ..
267: *
268: * ==========================================================================
269: *
270: * .. Parameters ..
271: COMPLEX*16 ONE, ZERO
272: PARAMETER ( ONE = (1.0,0.0), ZERO = (0.0,0.0) )
273: * ..
274: * .. Local Scalars ..
275: INTEGER I, J, MP, NP, KP
276: LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW
277: * ..
278: * .. External Functions ..
279: LOGICAL LSAME
280: EXTERNAL LSAME
281: * ..
282: * .. External Subroutines ..
283: EXTERNAL ZGEMM, ZTRMM
284: * ..
285: * .. Intrinsic Functions ..
286: INTRINSIC CONJG
287: * ..
288: * .. Executable Statements ..
289: *
290: * Quick return if possible
291: *
292: IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN
293: *
294: IF( LSAME( STOREV, 'C' ) ) THEN
295: COLUMN = .TRUE.
296: ROW = .FALSE.
297: ELSE IF ( LSAME( STOREV, 'R' ) ) THEN
298: COLUMN = .FALSE.
299: ROW = .TRUE.
300: ELSE
301: COLUMN = .FALSE.
302: ROW = .FALSE.
303: END IF
304: *
305: IF( LSAME( SIDE, 'L' ) ) THEN
306: LEFT = .TRUE.
307: RIGHT = .FALSE.
308: ELSE IF( LSAME( SIDE, 'R' ) ) THEN
309: LEFT = .FALSE.
310: RIGHT = .TRUE.
311: ELSE
312: LEFT = .FALSE.
313: RIGHT = .FALSE.
314: END IF
315: *
316: IF( LSAME( DIRECT, 'F' ) ) THEN
317: FORWARD = .TRUE.
318: BACKWARD = .FALSE.
319: ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
320: FORWARD = .FALSE.
321: BACKWARD = .TRUE.
322: ELSE
323: FORWARD = .FALSE.
324: BACKWARD = .FALSE.
325: END IF
326: *
327: * ---------------------------------------------------------------------------
1.7 bertrand 328: *
1.1 bertrand 329: IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN
330: *
331: * ---------------------------------------------------------------------------
332: *
333: * Let W = [ I ] (K-by-K)
334: * [ V ] (M-by-K)
335: *
336: * Form H C or H**H C where C = [ A ] (K-by-N)
337: * [ B ] (M-by-N)
338: *
339: * H = I - W T W**H or H**H = I - W T**H W**H
340: *
341: * A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
1.7 bertrand 342: * B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
1.1 bertrand 343: *
344: * ---------------------------------------------------------------------------
345: *
346: MP = MIN( M-L+1, M )
347: KP = MIN( L+1, K )
1.7 bertrand 348: *
1.1 bertrand 349: DO J = 1, N
350: DO I = 1, L
351: WORK( I, J ) = B( M-L+I, J )
352: END DO
353: END DO
354: CALL ZTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( MP, 1 ), LDV,
1.7 bertrand 355: $ WORK, LDWORK )
356: CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V, LDV, B, LDB,
1.1 bertrand 357: $ ONE, WORK, LDWORK )
1.7 bertrand 358: CALL ZGEMM( 'C', 'N', K-L, N, M, ONE, V( 1, KP ), LDV,
1.1 bertrand 359: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
1.7 bertrand 360: *
1.1 bertrand 361: DO J = 1, N
362: DO I = 1, K
363: WORK( I, J ) = WORK( I, J ) + A( I, J )
364: END DO
365: END DO
366: *
1.7 bertrand 367: CALL ZTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
1.1 bertrand 368: $ WORK, LDWORK )
1.7 bertrand 369: *
1.1 bertrand 370: DO J = 1, N
371: DO I = 1, K
372: A( I, J ) = A( I, J ) - WORK( I, J )
373: END DO
374: END DO
375: *
376: CALL ZGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
377: $ ONE, B, LDB )
378: CALL ZGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV,
1.7 bertrand 379: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
1.1 bertrand 380: CALL ZTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV,
381: $ WORK, LDWORK )
382: DO J = 1, N
383: DO I = 1, L
384: B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
385: END DO
386: END DO
387: *
388: * ---------------------------------------------------------------------------
1.7 bertrand 389: *
1.1 bertrand 390: ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN
391: *
392: * ---------------------------------------------------------------------------
393: *
394: * Let W = [ I ] (K-by-K)
395: * [ V ] (N-by-K)
396: *
397: * Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N)
398: *
399: * H = I - W T W**H or H**H = I - W T**H W**H
400: *
401: * A = A - (A + B V) T or A = A - (A + B V) T**H
402: * B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
403: *
404: * ---------------------------------------------------------------------------
405: *
406: NP = MIN( N-L+1, N )
407: KP = MIN( L+1, K )
1.7 bertrand 408: *
1.1 bertrand 409: DO J = 1, L
410: DO I = 1, M
411: WORK( I, J ) = B( I, N-L+J )
412: END DO
413: END DO
414: CALL ZTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV,
415: $ WORK, LDWORK )
1.7 bertrand 416: CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB,
1.1 bertrand 417: $ V, LDV, ONE, WORK, LDWORK )
1.7 bertrand 418: CALL ZGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
1.1 bertrand 419: $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK )
1.7 bertrand 420: *
1.1 bertrand 421: DO J = 1, K
422: DO I = 1, M
423: WORK( I, J ) = WORK( I, J ) + A( I, J )
424: END DO
425: END DO
426: *
1.7 bertrand 427: CALL ZTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
1.1 bertrand 428: $ WORK, LDWORK )
1.7 bertrand 429: *
1.1 bertrand 430: DO J = 1, K
431: DO I = 1, M
432: A( I, J ) = A( I, J ) - WORK( I, J )
433: END DO
434: END DO
435: *
436: CALL ZGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK,
437: $ V, LDV, ONE, B, LDB )
438: CALL ZGEMM( 'N', 'C', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
439: $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
440: CALL ZTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( NP, 1 ), LDV,
441: $ WORK, LDWORK )
442: DO J = 1, L
443: DO I = 1, M
444: B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
445: END DO
446: END DO
447: *
448: * ---------------------------------------------------------------------------
1.7 bertrand 449: *
1.1 bertrand 450: ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN
451: *
452: * ---------------------------------------------------------------------------
453: *
454: * Let W = [ V ] (M-by-K)
455: * [ I ] (K-by-K)
456: *
457: * Form H C or H**H C where C = [ B ] (M-by-N)
458: * [ A ] (K-by-N)
459: *
460: * H = I - W T W**H or H**H = I - W T**H W**H
461: *
462: * A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
1.7 bertrand 463: * B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
1.1 bertrand 464: *
465: * ---------------------------------------------------------------------------
466: *
467: MP = MIN( L+1, M )
468: KP = MIN( K-L+1, K )
469: *
470: DO J = 1, N
471: DO I = 1, L
472: WORK( K-L+I, J ) = B( I, J )
473: END DO
474: END DO
475: *
476: CALL ZTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, KP ), LDV,
477: $ WORK( KP, 1 ), LDWORK )
1.7 bertrand 478: CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V( MP, KP ), LDV,
1.1 bertrand 479: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
480: CALL ZGEMM( 'C', 'N', K-L, N, M, ONE, V, LDV,
1.7 bertrand 481: $ B, LDB, ZERO, WORK, LDWORK )
1.1 bertrand 482: *
483: DO J = 1, N
484: DO I = 1, K
485: WORK( I, J ) = WORK( I, J ) + A( I, J )
486: END DO
487: END DO
488: *
1.7 bertrand 489: CALL ZTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT,
1.1 bertrand 490: $ WORK, LDWORK )
1.7 bertrand 491: *
1.1 bertrand 492: DO J = 1, N
493: DO I = 1, K
494: A( I, J ) = A( I, J ) - WORK( I, J )
495: END DO
496: END DO
497: *
1.7 bertrand 498: CALL ZGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV,
1.1 bertrand 499: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
500: CALL ZGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV,
501: $ WORK, LDWORK, ONE, B, LDB )
502: CALL ZTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV,
503: $ WORK( KP, 1 ), LDWORK )
504: DO J = 1, N
505: DO I = 1, L
506: B( I, J ) = B( I, J ) - WORK( K-L+I, J )
507: END DO
508: END DO
509: *
510: * ---------------------------------------------------------------------------
1.7 bertrand 511: *
1.1 bertrand 512: ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN
513: *
514: * ---------------------------------------------------------------------------
515: *
516: * Let W = [ V ] (N-by-K)
517: * [ I ] (K-by-K)
518: *
519: * Form C H or C H**H where C = [ B A ] (B is M-by-N, A is M-by-K)
520: *
521: * H = I - W T W**H or H**H = I - W T**H W**H
522: *
523: * A = A - (A + B V) T or A = A - (A + B V) T**H
524: * B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
525: *
526: * ---------------------------------------------------------------------------
527: *
528: NP = MIN( L+1, N )
529: KP = MIN( K-L+1, K )
1.7 bertrand 530: *
1.1 bertrand 531: DO J = 1, L
532: DO I = 1, M
533: WORK( I, K-L+J ) = B( I, J )
534: END DO
535: END DO
536: CALL ZTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV,
537: $ WORK( 1, KP ), LDWORK )
1.7 bertrand 538: CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB,
1.1 bertrand 539: $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
1.7 bertrand 540: CALL ZGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
1.1 bertrand 541: $ V, LDV, ZERO, WORK, LDWORK )
1.7 bertrand 542: *
1.1 bertrand 543: DO J = 1, K
544: DO I = 1, M
545: WORK( I, J ) = WORK( I, J ) + A( I, J )
546: END DO
547: END DO
548: *
1.7 bertrand 549: CALL ZTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
1.1 bertrand 550: $ WORK, LDWORK )
1.7 bertrand 551: *
1.1 bertrand 552: DO J = 1, K
553: DO I = 1, M
554: A( I, J ) = A( I, J ) - WORK( I, J )
555: END DO
556: END DO
557: *
558: CALL ZGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK,
559: $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
560: CALL ZGEMM( 'N', 'C', M, L, K-L, -ONE, WORK, LDWORK,
561: $ V, LDV, ONE, B, LDB )
562: CALL ZTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, KP ), LDV,
563: $ WORK( 1, KP ), LDWORK )
564: DO J = 1, L
565: DO I = 1, M
566: B( I, J ) = B( I, J ) - WORK( I, K-L+J )
567: END DO
568: END DO
569: *
570: * ---------------------------------------------------------------------------
1.7 bertrand 571: *
1.1 bertrand 572: ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN
573: *
574: * ---------------------------------------------------------------------------
575: *
576: * Let W = [ I V ] ( I is K-by-K, V is K-by-M )
577: *
578: * Form H C or H**H C where C = [ A ] (K-by-N)
579: * [ B ] (M-by-N)
580: *
581: * H = I - W**H T W or H**H = I - W**H T**H W
582: *
583: * A = A - T (A + V B) or A = A - T**H (A + V B)
1.7 bertrand 584: * B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B)
1.1 bertrand 585: *
586: * ---------------------------------------------------------------------------
587: *
588: MP = MIN( M-L+1, M )
589: KP = MIN( L+1, K )
590: *
591: DO J = 1, N
592: DO I = 1, L
593: WORK( I, J ) = B( M-L+I, J )
594: END DO
1.7 bertrand 595: END DO
1.1 bertrand 596: CALL ZTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
597: $ WORK, LDB )
1.7 bertrand 598: CALL ZGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB,
1.1 bertrand 599: $ ONE, WORK, LDWORK )
1.7 bertrand 600: CALL ZGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV,
1.1 bertrand 601: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
602: *
603: DO J = 1, N
604: DO I = 1, K
605: WORK( I, J ) = WORK( I, J ) + A( I, J )
606: END DO
607: END DO
608: *
1.7 bertrand 609: CALL ZTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
1.1 bertrand 610: $ WORK, LDWORK )
611: *
612: DO J = 1, N
613: DO I = 1, K
614: A( I, J ) = A( I, J ) - WORK( I, J )
615: END DO
616: END DO
617: *
618: CALL ZGEMM( 'C', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
619: $ ONE, B, LDB )
1.7 bertrand 620: CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV,
1.1 bertrand 621: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
622: CALL ZTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, MP ), LDV,
623: $ WORK, LDWORK )
624: DO J = 1, N
625: DO I = 1, L
626: B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
627: END DO
628: END DO
629: *
630: * ---------------------------------------------------------------------------
1.7 bertrand 631: *
1.1 bertrand 632: ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN
633: *
634: * ---------------------------------------------------------------------------
635: *
636: * Let W = [ I V ] ( I is K-by-K, V is K-by-N )
637: *
638: * Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N)
639: *
640: * H = I - W**H T W or H**H = I - W**H T**H W
641: *
642: * A = A - (A + B V**H) T or A = A - (A + B V**H) T**H
643: * B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V
644: *
645: * ---------------------------------------------------------------------------
646: *
647: NP = MIN( N-L+1, N )
648: KP = MIN( L+1, K )
649: *
650: DO J = 1, L
651: DO I = 1, M
652: WORK( I, J ) = B( I, N-L+J )
653: END DO
654: END DO
655: CALL ZTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, NP ), LDV,
656: $ WORK, LDWORK )
657: CALL ZGEMM( 'N', 'C', M, L, N-L, ONE, B, LDB, V, LDV,
658: $ ONE, WORK, LDWORK )
1.7 bertrand 659: CALL ZGEMM( 'N', 'C', M, K-L, N, ONE, B, LDB,
1.1 bertrand 660: $ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK )
661: *
662: DO J = 1, K
663: DO I = 1, M
664: WORK( I, J ) = WORK( I, J ) + A( I, J )
665: END DO
666: END DO
667: *
1.7 bertrand 668: CALL ZTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
1.1 bertrand 669: $ WORK, LDWORK )
670: *
671: DO J = 1, K
672: DO I = 1, M
673: A( I, J ) = A( I, J ) - WORK( I, J )
674: END DO
675: END DO
676: *
1.7 bertrand 677: CALL ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
1.1 bertrand 678: $ V, LDV, ONE, B, LDB )
679: CALL ZGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
1.7 bertrand 680: $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )
1.1 bertrand 681: CALL ZTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV,
682: $ WORK, LDWORK )
683: DO J = 1, L
684: DO I = 1, M
685: B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
686: END DO
687: END DO
688: *
689: * ---------------------------------------------------------------------------
1.7 bertrand 690: *
1.1 bertrand 691: ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN
692: *
693: * ---------------------------------------------------------------------------
694: *
695: * Let W = [ V I ] ( I is K-by-K, V is K-by-M )
696: *
697: * Form H C or H**H C where C = [ B ] (M-by-N)
698: * [ A ] (K-by-N)
699: *
700: * H = I - W**H T W or H**H = I - W**H T**H W
701: *
702: * A = A - T (A + V B) or A = A - T**H (A + V B)
1.7 bertrand 703: * B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B)
1.1 bertrand 704: *
705: * ---------------------------------------------------------------------------
706: *
707: MP = MIN( L+1, M )
708: KP = MIN( K-L+1, K )
709: *
710: DO J = 1, N
711: DO I = 1, L
712: WORK( K-L+I, J ) = B( I, J )
713: END DO
714: END DO
715: CALL ZTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV,
716: $ WORK( KP, 1 ), LDWORK )
717: CALL ZGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV,
718: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
719: CALL ZGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB,
720: $ ZERO, WORK, LDWORK )
721: *
722: DO J = 1, N
723: DO I = 1, K
724: WORK( I, J ) = WORK( I, J ) + A( I, J )
725: END DO
726: END DO
727: *
728: CALL ZTRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT,
729: $ WORK, LDWORK )
730: *
731: DO J = 1, N
732: DO I = 1, K
733: A( I, J ) = A( I, J ) - WORK( I, J )
734: END DO
735: END DO
736: *
737: CALL ZGEMM( 'C', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV,
738: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
1.7 bertrand 739: CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V, LDV,
1.1 bertrand 740: $ WORK, LDWORK, ONE, B, LDB )
741: CALL ZTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( KP, 1 ), LDV,
1.7 bertrand 742: $ WORK( KP, 1 ), LDWORK )
1.1 bertrand 743: DO J = 1, N
744: DO I = 1, L
745: B( I, J ) = B( I, J ) - WORK( K-L+I, J )
746: END DO
747: END DO
748: *
749: * ---------------------------------------------------------------------------
1.7 bertrand 750: *
1.1 bertrand 751: ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN
752: *
753: * ---------------------------------------------------------------------------
754: *
755: * Let W = [ V I ] ( I is K-by-K, V is K-by-N )
756: *
757: * Form C H or C H**H where C = [ B A ] (A is M-by-K, B is M-by-N)
758: *
759: * H = I - W**H T W or H**H = I - W**H T**H W
760: *
761: * A = A - (A + B V**H) T or A = A - (A + B V**H) T**H
762: * B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V
763: *
764: * ---------------------------------------------------------------------------
765: *
766: NP = MIN( L+1, N )
767: KP = MIN( K-L+1, K )
768: *
769: DO J = 1, L
770: DO I = 1, M
771: WORK( I, K-L+J ) = B( I, J )
772: END DO
773: END DO
774: CALL ZTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( KP, 1 ), LDV,
775: $ WORK( 1, KP ), LDWORK )
776: CALL ZGEMM( 'N', 'C', M, L, N-L, ONE, B( 1, NP ), LDB,
777: $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
778: CALL ZGEMM( 'N', 'C', M, K-L, N, ONE, B, LDB, V, LDV,
1.7 bertrand 779: $ ZERO, WORK, LDWORK )
1.1 bertrand 780: *
781: DO J = 1, K
782: DO I = 1, M
783: WORK( I, J ) = WORK( I, J ) + A( I, J )
784: END DO
785: END DO
786: *
1.7 bertrand 787: CALL ZTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
1.1 bertrand 788: $ WORK, LDWORK )
789: *
790: DO J = 1, K
791: DO I = 1, M
792: A( I, J ) = A( I, J ) - WORK( I, J )
793: END DO
794: END DO
795: *
1.7 bertrand 796: CALL ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
1.1 bertrand 797: $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
1.7 bertrand 798: CALL ZGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,
1.1 bertrand 799: $ V, LDV, ONE, B, LDB )
800: CALL ZTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV,
801: $ WORK( 1, KP ), LDWORK )
802: DO J = 1, L
803: DO I = 1, M
804: B( I, J ) = B( I, J ) - WORK( I, K-L+J )
805: END DO
806: END DO
807: *
808: END IF
809: *
810: RETURN
811: *
812: * End of ZTPRFB
813: *
814: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztprfb.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack