Annotation of rpl/lapack/lapack/ztprfb.f, revision 1.11
1.11 ! bertrand 1: *> \brief \b ZTPRFB applies a complex "triangular-pentagonal" block reflector to a 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.10 bertrand 155: *> If SIDE = 'L', LDA >= max(1,K);
156: *> If SIDE = 'R', LDA >= 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: *
197: *> \ingroup complex16OTHERauxiliary
198: *
199: *> \par Further Details:
200: * =====================
201: *>
202: *> \verbatim
203: *>
204: *> The matrix C is a composite matrix formed from blocks A and B.
1.7 bertrand 205: *> The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K,
1.1 bertrand 206: *> and if SIDE = 'L', A is of size K-by-N.
207: *>
208: *> If SIDE = 'R' and DIRECT = 'F', C = [A B].
209: *>
1.7 bertrand 210: *> If SIDE = 'L' and DIRECT = 'F', C = [A]
1.1 bertrand 211: *> [B].
212: *>
213: *> If SIDE = 'R' and DIRECT = 'B', C = [B A].
214: *>
215: *> If SIDE = 'L' and DIRECT = 'B', C = [B]
1.7 bertrand 216: *> [A].
1.1 bertrand 217: *>
1.7 bertrand 218: *> The pentagonal matrix V is composed of a rectangular block V1 and a
219: *> trapezoidal block V2. The size of the trapezoidal block is determined by
1.1 bertrand 220: *> the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular;
221: *> if L=0, there is no trapezoidal block, thus V = V1 is rectangular.
222: *>
223: *> If DIRECT = 'F' and STOREV = 'C': V = [V1]
224: *> [V2]
225: *> - V2 is upper trapezoidal (first L rows of K-by-K upper triangular)
226: *>
227: *> If DIRECT = 'F' and STOREV = 'R': V = [V1 V2]
228: *>
229: *> - V2 is lower trapezoidal (first L columns of K-by-K lower triangular)
230: *>
231: *> If DIRECT = 'B' and STOREV = 'C': V = [V2]
232: *> [V1]
233: *> - V2 is lower trapezoidal (last L rows of K-by-K lower triangular)
234: *>
235: *> If DIRECT = 'B' and STOREV = 'R': V = [V2 V1]
1.7 bertrand 236: *>
1.1 bertrand 237: *> - V2 is upper trapezoidal (last L columns of K-by-K upper triangular)
238: *>
239: *> If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K.
240: *>
241: *> If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K.
242: *>
243: *> If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L.
244: *>
245: *> If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.
246: *> \endverbatim
247: *>
248: * =====================================================================
1.7 bertrand 249: SUBROUTINE ZTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
1.1 bertrand 250: $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
251: *
1.11 ! bertrand 252: * -- LAPACK auxiliary routine --
1.1 bertrand 253: * -- LAPACK is a software package provided by Univ. of Tennessee, --
254: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
255: *
256: * .. Scalar Arguments ..
257: CHARACTER DIRECT, SIDE, STOREV, TRANS
258: INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
259: * ..
260: * .. Array Arguments ..
1.7 bertrand 261: COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ),
1.1 bertrand 262: $ V( LDV, * ), WORK( LDWORK, * )
263: * ..
264: *
265: * ==========================================================================
266: *
267: * .. Parameters ..
268: COMPLEX*16 ONE, ZERO
269: PARAMETER ( ONE = (1.0,0.0), ZERO = (0.0,0.0) )
270: * ..
271: * .. Local Scalars ..
272: INTEGER I, J, MP, NP, KP
273: LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW
274: * ..
275: * .. External Functions ..
276: LOGICAL LSAME
277: EXTERNAL LSAME
278: * ..
279: * .. External Subroutines ..
280: EXTERNAL ZGEMM, ZTRMM
281: * ..
282: * .. Intrinsic Functions ..
283: INTRINSIC CONJG
284: * ..
285: * .. Executable Statements ..
286: *
287: * Quick return if possible
288: *
289: IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN
290: *
291: IF( LSAME( STOREV, 'C' ) ) THEN
292: COLUMN = .TRUE.
293: ROW = .FALSE.
294: ELSE IF ( LSAME( STOREV, 'R' ) ) THEN
295: COLUMN = .FALSE.
296: ROW = .TRUE.
297: ELSE
298: COLUMN = .FALSE.
299: ROW = .FALSE.
300: END IF
301: *
302: IF( LSAME( SIDE, 'L' ) ) THEN
303: LEFT = .TRUE.
304: RIGHT = .FALSE.
305: ELSE IF( LSAME( SIDE, 'R' ) ) THEN
306: LEFT = .FALSE.
307: RIGHT = .TRUE.
308: ELSE
309: LEFT = .FALSE.
310: RIGHT = .FALSE.
311: END IF
312: *
313: IF( LSAME( DIRECT, 'F' ) ) THEN
314: FORWARD = .TRUE.
315: BACKWARD = .FALSE.
316: ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
317: FORWARD = .FALSE.
318: BACKWARD = .TRUE.
319: ELSE
320: FORWARD = .FALSE.
321: BACKWARD = .FALSE.
322: END IF
323: *
324: * ---------------------------------------------------------------------------
1.7 bertrand 325: *
1.1 bertrand 326: IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN
327: *
328: * ---------------------------------------------------------------------------
329: *
330: * Let W = [ I ] (K-by-K)
331: * [ V ] (M-by-K)
332: *
333: * Form H C or H**H C where C = [ A ] (K-by-N)
334: * [ B ] (M-by-N)
335: *
336: * H = I - W T W**H or H**H = I - W T**H W**H
337: *
338: * A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
1.7 bertrand 339: * B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
1.1 bertrand 340: *
341: * ---------------------------------------------------------------------------
342: *
343: MP = MIN( M-L+1, M )
344: KP = MIN( L+1, K )
1.7 bertrand 345: *
1.1 bertrand 346: DO J = 1, N
347: DO I = 1, L
348: WORK( I, J ) = B( M-L+I, J )
349: END DO
350: END DO
351: CALL ZTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( MP, 1 ), LDV,
1.7 bertrand 352: $ WORK, LDWORK )
353: CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V, LDV, B, LDB,
1.1 bertrand 354: $ ONE, WORK, LDWORK )
1.7 bertrand 355: CALL ZGEMM( 'C', 'N', K-L, N, M, ONE, V( 1, KP ), LDV,
1.1 bertrand 356: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
1.7 bertrand 357: *
1.1 bertrand 358: DO J = 1, N
359: DO I = 1, K
360: WORK( I, J ) = WORK( I, J ) + A( I, J )
361: END DO
362: END DO
363: *
1.7 bertrand 364: CALL ZTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
1.1 bertrand 365: $ WORK, LDWORK )
1.7 bertrand 366: *
1.1 bertrand 367: DO J = 1, N
368: DO I = 1, K
369: A( I, J ) = A( I, J ) - WORK( I, J )
370: END DO
371: END DO
372: *
373: CALL ZGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
374: $ ONE, B, LDB )
375: CALL ZGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV,
1.7 bertrand 376: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
1.1 bertrand 377: CALL ZTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV,
378: $ WORK, LDWORK )
379: DO J = 1, N
380: DO I = 1, L
381: B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
382: END DO
383: END DO
384: *
385: * ---------------------------------------------------------------------------
1.7 bertrand 386: *
1.1 bertrand 387: ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN
388: *
389: * ---------------------------------------------------------------------------
390: *
391: * Let W = [ I ] (K-by-K)
392: * [ V ] (N-by-K)
393: *
394: * Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N)
395: *
396: * H = I - W T W**H or H**H = I - W T**H W**H
397: *
398: * A = A - (A + B V) T or A = A - (A + B V) T**H
399: * B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
400: *
401: * ---------------------------------------------------------------------------
402: *
403: NP = MIN( N-L+1, N )
404: KP = MIN( L+1, K )
1.7 bertrand 405: *
1.1 bertrand 406: DO J = 1, L
407: DO I = 1, M
408: WORK( I, J ) = B( I, N-L+J )
409: END DO
410: END DO
411: CALL ZTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV,
412: $ WORK, LDWORK )
1.7 bertrand 413: CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB,
1.1 bertrand 414: $ V, LDV, ONE, WORK, LDWORK )
1.7 bertrand 415: CALL ZGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
1.1 bertrand 416: $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK )
1.7 bertrand 417: *
1.1 bertrand 418: DO J = 1, K
419: DO I = 1, M
420: WORK( I, J ) = WORK( I, J ) + A( I, J )
421: END DO
422: END DO
423: *
1.7 bertrand 424: CALL ZTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
1.1 bertrand 425: $ WORK, LDWORK )
1.7 bertrand 426: *
1.1 bertrand 427: DO J = 1, K
428: DO I = 1, M
429: A( I, J ) = A( I, J ) - WORK( I, J )
430: END DO
431: END DO
432: *
433: CALL ZGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK,
434: $ V, LDV, ONE, B, LDB )
435: CALL ZGEMM( 'N', 'C', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
436: $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
437: CALL ZTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( NP, 1 ), LDV,
438: $ WORK, LDWORK )
439: DO J = 1, L
440: DO I = 1, M
441: B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
442: END DO
443: END DO
444: *
445: * ---------------------------------------------------------------------------
1.7 bertrand 446: *
1.1 bertrand 447: ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN
448: *
449: * ---------------------------------------------------------------------------
450: *
451: * Let W = [ V ] (M-by-K)
452: * [ I ] (K-by-K)
453: *
454: * Form H C or H**H C where C = [ B ] (M-by-N)
455: * [ A ] (K-by-N)
456: *
457: * H = I - W T W**H or H**H = I - W T**H W**H
458: *
459: * A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
1.7 bertrand 460: * B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
1.1 bertrand 461: *
462: * ---------------------------------------------------------------------------
463: *
464: MP = MIN( L+1, M )
465: KP = MIN( K-L+1, K )
466: *
467: DO J = 1, N
468: DO I = 1, L
469: WORK( K-L+I, J ) = B( I, J )
470: END DO
471: END DO
472: *
473: CALL ZTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, KP ), LDV,
474: $ WORK( KP, 1 ), LDWORK )
1.7 bertrand 475: CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V( MP, KP ), LDV,
1.1 bertrand 476: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
477: CALL ZGEMM( 'C', 'N', K-L, N, M, ONE, V, LDV,
1.7 bertrand 478: $ B, LDB, ZERO, WORK, LDWORK )
1.1 bertrand 479: *
480: DO J = 1, N
481: DO I = 1, K
482: WORK( I, J ) = WORK( I, J ) + A( I, J )
483: END DO
484: END DO
485: *
1.7 bertrand 486: CALL ZTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT,
1.1 bertrand 487: $ WORK, LDWORK )
1.7 bertrand 488: *
1.1 bertrand 489: DO J = 1, N
490: DO I = 1, K
491: A( I, J ) = A( I, J ) - WORK( I, J )
492: END DO
493: END DO
494: *
1.7 bertrand 495: CALL ZGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV,
1.1 bertrand 496: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
497: CALL ZGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV,
498: $ WORK, LDWORK, ONE, B, LDB )
499: CALL ZTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV,
500: $ WORK( KP, 1 ), LDWORK )
501: DO J = 1, N
502: DO I = 1, L
503: B( I, J ) = B( I, J ) - WORK( K-L+I, J )
504: END DO
505: END DO
506: *
507: * ---------------------------------------------------------------------------
1.7 bertrand 508: *
1.1 bertrand 509: ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN
510: *
511: * ---------------------------------------------------------------------------
512: *
513: * Let W = [ V ] (N-by-K)
514: * [ I ] (K-by-K)
515: *
516: * Form C H or C H**H where C = [ B A ] (B is M-by-N, A is M-by-K)
517: *
518: * H = I - W T W**H or H**H = I - W T**H W**H
519: *
520: * A = A - (A + B V) T or A = A - (A + B V) T**H
521: * B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
522: *
523: * ---------------------------------------------------------------------------
524: *
525: NP = MIN( L+1, N )
526: KP = MIN( K-L+1, K )
1.7 bertrand 527: *
1.1 bertrand 528: DO J = 1, L
529: DO I = 1, M
530: WORK( I, K-L+J ) = B( I, J )
531: END DO
532: END DO
533: CALL ZTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV,
534: $ WORK( 1, KP ), LDWORK )
1.7 bertrand 535: CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB,
1.1 bertrand 536: $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
1.7 bertrand 537: CALL ZGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
1.1 bertrand 538: $ V, LDV, ZERO, WORK, LDWORK )
1.7 bertrand 539: *
1.1 bertrand 540: DO J = 1, K
541: DO I = 1, M
542: WORK( I, J ) = WORK( I, J ) + A( I, J )
543: END DO
544: END DO
545: *
1.7 bertrand 546: CALL ZTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
1.1 bertrand 547: $ WORK, LDWORK )
1.7 bertrand 548: *
1.1 bertrand 549: DO J = 1, K
550: DO I = 1, M
551: A( I, J ) = A( I, J ) - WORK( I, J )
552: END DO
553: END DO
554: *
555: CALL ZGEMM( 'N', 'C', M, N-L, K, -ONE, WORK, LDWORK,
556: $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
557: CALL ZGEMM( 'N', 'C', M, L, K-L, -ONE, WORK, LDWORK,
558: $ V, LDV, ONE, B, LDB )
559: CALL ZTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, KP ), LDV,
560: $ WORK( 1, KP ), LDWORK )
561: DO J = 1, L
562: DO I = 1, M
563: B( I, J ) = B( I, J ) - WORK( I, K-L+J )
564: END DO
565: END DO
566: *
567: * ---------------------------------------------------------------------------
1.7 bertrand 568: *
1.1 bertrand 569: ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN
570: *
571: * ---------------------------------------------------------------------------
572: *
573: * Let W = [ I V ] ( I is K-by-K, V is K-by-M )
574: *
575: * Form H C or H**H C where C = [ A ] (K-by-N)
576: * [ B ] (M-by-N)
577: *
578: * H = I - W**H T W or H**H = I - W**H T**H W
579: *
580: * A = A - T (A + V B) or A = A - T**H (A + V B)
1.7 bertrand 581: * B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B)
1.1 bertrand 582: *
583: * ---------------------------------------------------------------------------
584: *
585: MP = MIN( M-L+1, M )
586: KP = MIN( L+1, K )
587: *
588: DO J = 1, N
589: DO I = 1, L
590: WORK( I, J ) = B( M-L+I, J )
591: END DO
1.7 bertrand 592: END DO
1.1 bertrand 593: CALL ZTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
594: $ WORK, LDB )
1.7 bertrand 595: CALL ZGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB,
1.1 bertrand 596: $ ONE, WORK, LDWORK )
1.7 bertrand 597: CALL ZGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV,
1.1 bertrand 598: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
599: *
600: DO J = 1, N
601: DO I = 1, K
602: WORK( I, J ) = WORK( I, J ) + A( I, J )
603: END DO
604: END DO
605: *
1.7 bertrand 606: CALL ZTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
1.1 bertrand 607: $ WORK, LDWORK )
608: *
609: DO J = 1, N
610: DO I = 1, K
611: A( I, J ) = A( I, J ) - WORK( I, J )
612: END DO
613: END DO
614: *
615: CALL ZGEMM( 'C', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
616: $ ONE, B, LDB )
1.7 bertrand 617: CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV,
1.1 bertrand 618: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
619: CALL ZTRMM( 'L', 'L', 'C', 'N', L, N, ONE, V( 1, MP ), LDV,
620: $ WORK, LDWORK )
621: DO J = 1, N
622: DO I = 1, L
623: B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
624: END DO
625: END DO
626: *
627: * ---------------------------------------------------------------------------
1.7 bertrand 628: *
1.1 bertrand 629: ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN
630: *
631: * ---------------------------------------------------------------------------
632: *
633: * Let W = [ I V ] ( I is K-by-K, V is K-by-N )
634: *
635: * Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N)
636: *
637: * H = I - W**H T W or H**H = I - W**H T**H W
638: *
639: * A = A - (A + B V**H) T or A = A - (A + B V**H) T**H
640: * B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V
641: *
642: * ---------------------------------------------------------------------------
643: *
644: NP = MIN( N-L+1, N )
645: KP = MIN( L+1, K )
646: *
647: DO J = 1, L
648: DO I = 1, M
649: WORK( I, J ) = B( I, N-L+J )
650: END DO
651: END DO
652: CALL ZTRMM( 'R', 'L', 'C', 'N', M, L, ONE, V( 1, NP ), LDV,
653: $ WORK, LDWORK )
654: CALL ZGEMM( 'N', 'C', M, L, N-L, ONE, B, LDB, V, LDV,
655: $ ONE, WORK, LDWORK )
1.7 bertrand 656: CALL ZGEMM( 'N', 'C', M, K-L, N, ONE, B, LDB,
1.1 bertrand 657: $ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK )
658: *
659: DO J = 1, K
660: DO I = 1, M
661: WORK( I, J ) = WORK( I, J ) + A( I, J )
662: END DO
663: END DO
664: *
1.7 bertrand 665: CALL ZTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
1.1 bertrand 666: $ WORK, LDWORK )
667: *
668: DO J = 1, K
669: DO I = 1, M
670: A( I, J ) = A( I, J ) - WORK( I, J )
671: END DO
672: END DO
673: *
1.7 bertrand 674: CALL ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
1.1 bertrand 675: $ V, LDV, ONE, B, LDB )
676: CALL ZGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
1.7 bertrand 677: $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )
1.1 bertrand 678: CALL ZTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV,
679: $ WORK, LDWORK )
680: DO J = 1, L
681: DO I = 1, M
682: B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
683: END DO
684: END DO
685: *
686: * ---------------------------------------------------------------------------
1.7 bertrand 687: *
1.1 bertrand 688: ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN
689: *
690: * ---------------------------------------------------------------------------
691: *
692: * Let W = [ V I ] ( I is K-by-K, V is K-by-M )
693: *
694: * Form H C or H**H C where C = [ B ] (M-by-N)
695: * [ A ] (K-by-N)
696: *
697: * H = I - W**H T W or H**H = I - W**H T**H W
698: *
699: * A = A - T (A + V B) or A = A - T**H (A + V B)
1.7 bertrand 700: * B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B)
1.1 bertrand 701: *
702: * ---------------------------------------------------------------------------
703: *
704: MP = MIN( L+1, M )
705: KP = MIN( K-L+1, K )
706: *
707: DO J = 1, N
708: DO I = 1, L
709: WORK( K-L+I, J ) = B( I, J )
710: END DO
711: END DO
712: CALL ZTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV,
713: $ WORK( KP, 1 ), LDWORK )
714: CALL ZGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV,
715: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
716: CALL ZGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB,
717: $ ZERO, WORK, LDWORK )
718: *
719: DO J = 1, N
720: DO I = 1, K
721: WORK( I, J ) = WORK( I, J ) + A( I, J )
722: END DO
723: END DO
724: *
725: CALL ZTRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT,
726: $ WORK, LDWORK )
727: *
728: DO J = 1, N
729: DO I = 1, K
730: A( I, J ) = A( I, J ) - WORK( I, J )
731: END DO
732: END DO
733: *
734: CALL ZGEMM( 'C', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV,
735: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
1.7 bertrand 736: CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V, LDV,
1.1 bertrand 737: $ WORK, LDWORK, ONE, B, LDB )
738: CALL ZTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( KP, 1 ), LDV,
1.7 bertrand 739: $ WORK( KP, 1 ), LDWORK )
1.1 bertrand 740: DO J = 1, N
741: DO I = 1, L
742: B( I, J ) = B( I, J ) - WORK( K-L+I, J )
743: END DO
744: END DO
745: *
746: * ---------------------------------------------------------------------------
1.7 bertrand 747: *
1.1 bertrand 748: ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN
749: *
750: * ---------------------------------------------------------------------------
751: *
752: * Let W = [ V I ] ( I is K-by-K, V is K-by-N )
753: *
754: * Form C H or C H**H where C = [ B A ] (A is M-by-K, B is M-by-N)
755: *
756: * H = I - W**H T W or H**H = I - W**H T**H W
757: *
758: * A = A - (A + B V**H) T or A = A - (A + B V**H) T**H
759: * B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V
760: *
761: * ---------------------------------------------------------------------------
762: *
763: NP = MIN( L+1, N )
764: KP = MIN( K-L+1, K )
765: *
766: DO J = 1, L
767: DO I = 1, M
768: WORK( I, K-L+J ) = B( I, J )
769: END DO
770: END DO
771: CALL ZTRMM( 'R', 'U', 'C', 'N', M, L, ONE, V( KP, 1 ), LDV,
772: $ WORK( 1, KP ), LDWORK )
773: CALL ZGEMM( 'N', 'C', M, L, N-L, ONE, B( 1, NP ), LDB,
774: $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
775: CALL ZGEMM( 'N', 'C', M, K-L, N, ONE, B, LDB, V, LDV,
1.7 bertrand 776: $ ZERO, WORK, LDWORK )
1.1 bertrand 777: *
778: DO J = 1, K
779: DO I = 1, M
780: WORK( I, J ) = WORK( I, J ) + A( I, J )
781: END DO
782: END DO
783: *
1.7 bertrand 784: CALL ZTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
1.1 bertrand 785: $ WORK, LDWORK )
786: *
787: DO J = 1, K
788: DO I = 1, M
789: A( I, J ) = A( I, J ) - WORK( I, J )
790: END DO
791: END DO
792: *
1.7 bertrand 793: CALL ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
1.1 bertrand 794: $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
1.7 bertrand 795: CALL ZGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,
1.1 bertrand 796: $ V, LDV, ONE, B, LDB )
797: CALL ZTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV,
798: $ WORK( 1, KP ), LDWORK )
799: DO J = 1, L
800: DO I = 1, M
801: B( I, J ) = B( I, J ) - WORK( I, K-L+J )
802: END DO
803: END DO
804: *
805: END IF
806: *
807: RETURN
808: *
809: * End of ZTPRFB
810: *
811: 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