Annotation of rpl/lapack/lapack/ztprfb.f, revision 1.5
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: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
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">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
22: * V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER DIRECT, SIDE, STOREV, TRANS
26: * INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
27: * ..
28: * .. Array Arguments ..
29: * COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ),
30: * $ V( LDV, * ), WORK( LDWORK, * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
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
41: *> blocks A and B, either from the left or right.
42: *>
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
83: *> The number of rows of the matrix B.
84: *> M >= 0.
85: *> \endverbatim
86: *>
87: *> \param[in] N
88: *> \verbatim
89: *> N is INTEGER
90: *> The number of columns of the matrix B.
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
98: *> reflectors whose product defines the block reflector.
99: *> K >= 0.
100: *> \endverbatim
101: *>
102: *> \param[in] L
103: *> \verbatim
104: *> L is INTEGER
105: *> The order of the trapezoidal part of V.
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
132: *> block reflector.
133: *> \endverbatim
134: *>
135: *> \param[in] LDT
136: *> \verbatim
137: *> LDT is INTEGER
138: *> The leading dimension of the array T.
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.
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 Futher Details.
149: *> \endverbatim
150: *>
151: *> \param[in] LDA
152: *> \verbatim
153: *> LDA is INTEGER
154: *> The leading dimension of the array A.
155: *> If SIDE = 'L', LDC >= max(1,K);
156: *> If SIDE = 'R', LDC >= max(1,M).
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
170: *> The leading dimension of the array B.
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.
185: *> If SIDE = 'L', LDWORK >= K;
186: *> if SIDE = 'R', LDWORK >= M.
187: *> \endverbatim
188: *
189: * Authors:
190: * ========
191: *
192: *> \author Univ. of Tennessee
193: *> \author Univ. of California Berkeley
194: *> \author Univ. of Colorado Denver
195: *> \author NAG Ltd.
196: *
1.3 bertrand 197: *> \date September 2012
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.
207: *> The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K,
208: *> and if SIDE = 'L', A is of size K-by-N.
209: *>
210: *> If SIDE = 'R' and DIRECT = 'F', C = [A B].
211: *>
212: *> If SIDE = 'L' and DIRECT = 'F', C = [A]
213: *> [B].
214: *>
215: *> If SIDE = 'R' and DIRECT = 'B', C = [B A].
216: *>
217: *> If SIDE = 'L' and DIRECT = 'B', C = [B]
218: *> [A].
219: *>
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
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]
238: *>
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: * =====================================================================
251: SUBROUTINE ZTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
252: $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
253: *
1.3 bertrand 254: * -- LAPACK auxiliary routine (version 3.4.2) --
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.3 bertrand 257: * September 2012
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 ..
264: COMPLEX*16 A( LDA, * ), B( LDB, * ), T( LDT, * ),
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: * ---------------------------------------------------------------------------
328: *
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)
342: * B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
343: *
344: * ---------------------------------------------------------------------------
345: *
346: MP = MIN( M-L+1, M )
347: KP = MIN( L+1, K )
348: *
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,
355: $ WORK, LDWORK )
356: CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V, LDV, B, LDB,
357: $ ONE, WORK, LDWORK )
358: CALL ZGEMM( 'C', 'N', K-L, N, M, ONE, V( 1, KP ), LDV,
359: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
360: *
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: *
367: CALL ZTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
368: $ WORK, LDWORK )
369: *
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,
379: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
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: * ---------------------------------------------------------------------------
389: *
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 )
408: *
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 )
416: CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB,
417: $ V, LDV, ONE, WORK, LDWORK )
418: CALL ZGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
419: $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK )
420: *
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: *
427: CALL ZTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
428: $ WORK, LDWORK )
429: *
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: * ---------------------------------------------------------------------------
449: *
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)
463: * B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
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 )
478: CALL ZGEMM( 'C', 'N', L, N, M-L, ONE, V( MP, KP ), LDV,
479: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
480: CALL ZGEMM( 'C', 'N', K-L, N, M, ONE, V, LDV,
481: $ B, LDB, ZERO, WORK, LDWORK )
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: *
489: CALL ZTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT,
490: $ WORK, LDWORK )
491: *
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: *
498: CALL ZGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV,
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: * ---------------------------------------------------------------------------
511: *
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 )
530: *
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 )
538: CALL ZGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB,
539: $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
540: CALL ZGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
541: $ V, LDV, ZERO, WORK, LDWORK )
542: *
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: *
549: CALL ZTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
550: $ WORK, LDWORK )
551: *
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: * ---------------------------------------------------------------------------
571: *
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)
584: * B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B)
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
595: END DO
596: CALL ZTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
597: $ WORK, LDB )
598: CALL ZGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB,
599: $ ONE, WORK, LDWORK )
600: CALL ZGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV,
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: *
609: CALL ZTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
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 )
620: CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV,
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: * ---------------------------------------------------------------------------
631: *
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 )
659: CALL ZGEMM( 'N', 'C', M, K-L, N, ONE, B, LDB,
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: *
668: CALL ZTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
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: *
677: CALL ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
678: $ V, LDV, ONE, B, LDB )
679: CALL ZGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
680: $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )
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: * ---------------------------------------------------------------------------
690: *
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)
703: * B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B)
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 )
739: CALL ZGEMM( 'C', 'N', L, N, K-L, -ONE, V, LDV,
740: $ WORK, LDWORK, ONE, B, LDB )
741: CALL ZTRMM( 'L', 'U', 'C', 'N', L, N, ONE, V( KP, 1 ), LDV,
742: $ WORK( KP, 1 ), LDWORK )
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: * ---------------------------------------------------------------------------
750: *
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,
779: $ ZERO, WORK, LDWORK )
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: *
787: CALL ZTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
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: *
796: CALL ZGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
797: $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
798: CALL ZGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,
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