1: *> \brief \b DTPRFB applies a real "triangular-pentagonal" block reflector to a real matrix, which is composed of two blocks.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DTPRFB + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtprfb.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtprfb.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtprfb.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DTPRFB( 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: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
30: * $ V( LDV, * ), WORK( LDWORK, * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DTPRFB applies a real "triangular-pentagonal" block reflector H or its
40: *> transpose H**T to a real 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**T from the Left
52: *> = 'R': apply H or H**T from the Right
53: *> \endverbatim
54: *>
55: *> \param[in] TRANS
56: *> \verbatim
57: *> TRANS is CHARACTER*1
58: *> = 'N': apply H (No transpose)
59: *> = 'T': apply H**T (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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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**T*C or C*H or C*H**T. See Further 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', LDA >= max(1,K);
156: *> If SIDE = 'R', LDA >= max(1,M).
157: *> \endverbatim
158: *>
159: *> \param[in,out] B
160: *> \verbatim
161: *> B is DOUBLE PRECISION 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**T*C or C*H or C*H**T. 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 DOUBLE PRECISION 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: *
197: *> \ingroup doubleOTHERauxiliary
198: *
199: *> \par Further Details:
200: * =====================
201: *>
202: *> \verbatim
203: *>
204: *> The matrix C is a composite matrix formed from blocks A and B.
205: *> The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K,
206: *> and if SIDE = 'L', A is of size K-by-N.
207: *>
208: *> If SIDE = 'R' and DIRECT = 'F', C = [A B].
209: *>
210: *> If SIDE = 'L' and DIRECT = 'F', C = [A]
211: *> [B].
212: *>
213: *> If SIDE = 'R' and DIRECT = 'B', C = [B A].
214: *>
215: *> If SIDE = 'L' and DIRECT = 'B', C = [B]
216: *> [A].
217: *>
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
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]
236: *>
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: * =====================================================================
249: SUBROUTINE DTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
250: $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
251: *
252: * -- LAPACK auxiliary routine --
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 ..
261: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
262: $ V( LDV, * ), WORK( LDWORK, * )
263: * ..
264: *
265: * ==========================================================================
266: *
267: * .. Parameters ..
268: DOUBLE PRECISION ONE, ZERO
269: PARAMETER ( ONE = 1.0, ZERO = 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 DGEMM, DTRMM
281: * ..
282: * .. Executable Statements ..
283: *
284: * Quick return if possible
285: *
286: IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 .OR. L.LT.0 ) RETURN
287: *
288: IF( LSAME( STOREV, 'C' ) ) THEN
289: COLUMN = .TRUE.
290: ROW = .FALSE.
291: ELSE IF ( LSAME( STOREV, 'R' ) ) THEN
292: COLUMN = .FALSE.
293: ROW = .TRUE.
294: ELSE
295: COLUMN = .FALSE.
296: ROW = .FALSE.
297: END IF
298: *
299: IF( LSAME( SIDE, 'L' ) ) THEN
300: LEFT = .TRUE.
301: RIGHT = .FALSE.
302: ELSE IF( LSAME( SIDE, 'R' ) ) THEN
303: LEFT = .FALSE.
304: RIGHT = .TRUE.
305: ELSE
306: LEFT = .FALSE.
307: RIGHT = .FALSE.
308: END IF
309: *
310: IF( LSAME( DIRECT, 'F' ) ) THEN
311: FORWARD = .TRUE.
312: BACKWARD = .FALSE.
313: ELSE IF( LSAME( DIRECT, 'B' ) ) THEN
314: FORWARD = .FALSE.
315: BACKWARD = .TRUE.
316: ELSE
317: FORWARD = .FALSE.
318: BACKWARD = .FALSE.
319: END IF
320: *
321: * ---------------------------------------------------------------------------
322: *
323: IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN
324: *
325: * ---------------------------------------------------------------------------
326: *
327: * Let W = [ I ] (K-by-K)
328: * [ V ] (M-by-K)
329: *
330: * Form H C or H**T C where C = [ A ] (K-by-N)
331: * [ B ] (M-by-N)
332: *
333: * H = I - W T W**T or H**T = I - W T**T W**T
334: *
335: * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
336: * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
337: *
338: * ---------------------------------------------------------------------------
339: *
340: MP = MIN( M-L+1, M )
341: KP = MIN( L+1, K )
342: *
343: DO J = 1, N
344: DO I = 1, L
345: WORK( I, J ) = B( M-L+I, J )
346: END DO
347: END DO
348: CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( MP, 1 ), LDV,
349: $ WORK, LDWORK )
350: CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V, LDV, B, LDB,
351: $ ONE, WORK, LDWORK )
352: CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V( 1, KP ), LDV,
353: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
354: *
355: DO J = 1, N
356: DO I = 1, K
357: WORK( I, J ) = WORK( I, J ) + A( I, J )
358: END DO
359: END DO
360: *
361: CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
362: $ WORK, LDWORK )
363: *
364: DO J = 1, N
365: DO I = 1, K
366: A( I, J ) = A( I, J ) - WORK( I, J )
367: END DO
368: END DO
369: *
370: CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
371: $ ONE, B, LDB )
372: CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV,
373: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
374: CALL DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV,
375: $ WORK, LDWORK )
376: DO J = 1, N
377: DO I = 1, L
378: B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
379: END DO
380: END DO
381: *
382: * ---------------------------------------------------------------------------
383: *
384: ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN
385: *
386: * ---------------------------------------------------------------------------
387: *
388: * Let W = [ I ] (K-by-K)
389: * [ V ] (N-by-K)
390: *
391: * Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N)
392: *
393: * H = I - W T W**T or H**T = I - W T**T W**T
394: *
395: * A = A - (A + B V) T or A = A - (A + B V) T**T
396: * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
397: *
398: * ---------------------------------------------------------------------------
399: *
400: NP = MIN( N-L+1, N )
401: KP = MIN( L+1, K )
402: *
403: DO J = 1, L
404: DO I = 1, M
405: WORK( I, J ) = B( I, N-L+J )
406: END DO
407: END DO
408: CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV,
409: $ WORK, LDWORK )
410: CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB,
411: $ V, LDV, ONE, WORK, LDWORK )
412: CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
413: $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK )
414: *
415: DO J = 1, K
416: DO I = 1, M
417: WORK( I, J ) = WORK( I, J ) + A( I, J )
418: END DO
419: END DO
420: *
421: CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
422: $ WORK, LDWORK )
423: *
424: DO J = 1, K
425: DO I = 1, M
426: A( I, J ) = A( I, J ) - WORK( I, J )
427: END DO
428: END DO
429: *
430: CALL DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
431: $ V, LDV, ONE, B, LDB )
432: CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
433: $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
434: CALL DTRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( NP, 1 ), LDV,
435: $ WORK, LDWORK )
436: DO J = 1, L
437: DO I = 1, M
438: B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
439: END DO
440: END DO
441: *
442: * ---------------------------------------------------------------------------
443: *
444: ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN
445: *
446: * ---------------------------------------------------------------------------
447: *
448: * Let W = [ V ] (M-by-K)
449: * [ I ] (K-by-K)
450: *
451: * Form H C or H**T C where C = [ B ] (M-by-N)
452: * [ A ] (K-by-N)
453: *
454: * H = I - W T W**T or H**T = I - W T**T W**T
455: *
456: * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
457: * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
458: *
459: * ---------------------------------------------------------------------------
460: *
461: MP = MIN( L+1, M )
462: KP = MIN( K-L+1, K )
463: *
464: DO J = 1, N
465: DO I = 1, L
466: WORK( K-L+I, J ) = B( I, J )
467: END DO
468: END DO
469: *
470: CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, KP ), LDV,
471: $ WORK( KP, 1 ), LDWORK )
472: CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V( MP, KP ), LDV,
473: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
474: CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V, LDV,
475: $ B, LDB, ZERO, WORK, LDWORK )
476: *
477: DO J = 1, N
478: DO I = 1, K
479: WORK( I, J ) = WORK( I, J ) + A( I, J )
480: END DO
481: END DO
482: *
483: CALL DTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT,
484: $ WORK, LDWORK )
485: *
486: DO J = 1, N
487: DO I = 1, K
488: A( I, J ) = A( I, J ) - WORK( I, J )
489: END DO
490: END DO
491: *
492: CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV,
493: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
494: CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV,
495: $ WORK, LDWORK, ONE, B, LDB )
496: CALL DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, KP ), LDV,
497: $ WORK( KP, 1 ), LDWORK )
498: DO J = 1, N
499: DO I = 1, L
500: B( I, J ) = B( I, J ) - WORK( K-L+I, J )
501: END DO
502: END DO
503: *
504: * ---------------------------------------------------------------------------
505: *
506: ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN
507: *
508: * ---------------------------------------------------------------------------
509: *
510: * Let W = [ V ] (N-by-K)
511: * [ I ] (K-by-K)
512: *
513: * Form C H or C H**T where C = [ B A ] (B is M-by-N, A is M-by-K)
514: *
515: * H = I - W T W**T or H**T = I - W T**T W**T
516: *
517: * A = A - (A + B V) T or A = A - (A + B V) T**T
518: * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
519: *
520: * ---------------------------------------------------------------------------
521: *
522: NP = MIN( L+1, N )
523: KP = MIN( K-L+1, K )
524: *
525: DO J = 1, L
526: DO I = 1, M
527: WORK( I, K-L+J ) = B( I, J )
528: END DO
529: END DO
530: CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV,
531: $ WORK( 1, KP ), LDWORK )
532: CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB,
533: $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
534: CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
535: $ V, LDV, ZERO, WORK, LDWORK )
536: *
537: DO J = 1, K
538: DO I = 1, M
539: WORK( I, J ) = WORK( I, J ) + A( I, J )
540: END DO
541: END DO
542: *
543: CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
544: $ WORK, LDWORK )
545: *
546: DO J = 1, K
547: DO I = 1, M
548: A( I, J ) = A( I, J ) - WORK( I, J )
549: END DO
550: END DO
551: *
552: CALL DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
553: $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
554: CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK, LDWORK,
555: $ V, LDV, ONE, B, LDB )
556: CALL DTRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, KP ), LDV,
557: $ WORK( 1, KP ), LDWORK )
558: DO J = 1, L
559: DO I = 1, M
560: B( I, J ) = B( I, J ) - WORK( I, K-L+J )
561: END DO
562: END DO
563: *
564: * ---------------------------------------------------------------------------
565: *
566: ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN
567: *
568: * ---------------------------------------------------------------------------
569: *
570: * Let W = [ I V ] ( I is K-by-K, V is K-by-M )
571: *
572: * Form H C or H**T C where C = [ A ] (K-by-N)
573: * [ B ] (M-by-N)
574: *
575: * H = I - W**T T W or H**T = I - W**T T**T W
576: *
577: * A = A - T (A + V B) or A = A - T**T (A + V B)
578: * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
579: *
580: * ---------------------------------------------------------------------------
581: *
582: MP = MIN( M-L+1, M )
583: KP = MIN( L+1, K )
584: *
585: DO J = 1, N
586: DO I = 1, L
587: WORK( I, J ) = B( M-L+I, J )
588: END DO
589: END DO
590: CALL DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
591: $ WORK, LDB )
592: CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB,
593: $ ONE, WORK, LDWORK )
594: CALL DGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV,
595: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
596: *
597: DO J = 1, N
598: DO I = 1, K
599: WORK( I, J ) = WORK( I, J ) + A( I, J )
600: END DO
601: END DO
602: *
603: CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
604: $ WORK, LDWORK )
605: *
606: DO J = 1, N
607: DO I = 1, K
608: A( I, J ) = A( I, J ) - WORK( I, J )
609: END DO
610: END DO
611: *
612: CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
613: $ ONE, B, LDB )
614: CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV,
615: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
616: CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, MP ), LDV,
617: $ WORK, LDWORK )
618: DO J = 1, N
619: DO I = 1, L
620: B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
621: END DO
622: END DO
623: *
624: * ---------------------------------------------------------------------------
625: *
626: ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN
627: *
628: * ---------------------------------------------------------------------------
629: *
630: * Let W = [ I V ] ( I is K-by-K, V is K-by-N )
631: *
632: * Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N)
633: *
634: * H = I - W**T T W or H**T = I - W**T T**T W
635: *
636: * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
637: * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V
638: *
639: * ---------------------------------------------------------------------------
640: *
641: NP = MIN( N-L+1, N )
642: KP = MIN( L+1, K )
643: *
644: DO J = 1, L
645: DO I = 1, M
646: WORK( I, J ) = B( I, N-L+J )
647: END DO
648: END DO
649: CALL DTRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, NP ), LDV,
650: $ WORK, LDWORK )
651: CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B, LDB, V, LDV,
652: $ ONE, WORK, LDWORK )
653: CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB,
654: $ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK )
655: *
656: DO J = 1, K
657: DO I = 1, M
658: WORK( I, J ) = WORK( I, J ) + A( I, J )
659: END DO
660: END DO
661: *
662: CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
663: $ WORK, LDWORK )
664: *
665: DO J = 1, K
666: DO I = 1, M
667: A( I, J ) = A( I, J ) - WORK( I, J )
668: END DO
669: END DO
670: *
671: CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
672: $ V, LDV, ONE, B, LDB )
673: CALL DGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
674: $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )
675: CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV,
676: $ WORK, LDWORK )
677: DO J = 1, L
678: DO I = 1, M
679: B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
680: END DO
681: END DO
682: *
683: * ---------------------------------------------------------------------------
684: *
685: ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN
686: *
687: * ---------------------------------------------------------------------------
688: *
689: * Let W = [ V I ] ( I is K-by-K, V is K-by-M )
690: *
691: * Form H C or H**T C where C = [ B ] (M-by-N)
692: * [ A ] (K-by-N)
693: *
694: * H = I - W**T T W or H**T = I - W**T T**T W
695: *
696: * A = A - T (A + V B) or A = A - T**T (A + V B)
697: * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
698: *
699: * ---------------------------------------------------------------------------
700: *
701: MP = MIN( L+1, M )
702: KP = MIN( K-L+1, K )
703: *
704: DO J = 1, N
705: DO I = 1, L
706: WORK( K-L+I, J ) = B( I, J )
707: END DO
708: END DO
709: CALL DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV,
710: $ WORK( KP, 1 ), LDWORK )
711: CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV,
712: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
713: CALL DGEMM( 'N', 'N', K-L, N, M, ONE, V, LDV, B, LDB,
714: $ ZERO, WORK, LDWORK )
715: *
716: DO J = 1, N
717: DO I = 1, K
718: WORK( I, J ) = WORK( I, J ) + A( I, J )
719: END DO
720: END DO
721: *
722: CALL DTRMM( 'L', 'L ', TRANS, 'N', K, N, ONE, T, LDT,
723: $ WORK, LDWORK )
724: *
725: DO J = 1, N
726: DO I = 1, K
727: A( I, J ) = A( I, J ) - WORK( I, J )
728: END DO
729: END DO
730: *
731: CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV,
732: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
733: CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V, LDV,
734: $ WORK, LDWORK, ONE, B, LDB )
735: CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( KP, 1 ), LDV,
736: $ WORK( KP, 1 ), LDWORK )
737: DO J = 1, N
738: DO I = 1, L
739: B( I, J ) = B( I, J ) - WORK( K-L+I, J )
740: END DO
741: END DO
742: *
743: * ---------------------------------------------------------------------------
744: *
745: ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN
746: *
747: * ---------------------------------------------------------------------------
748: *
749: * Let W = [ V I ] ( I is K-by-K, V is K-by-N )
750: *
751: * Form C H or C H**T where C = [ B A ] (A is M-by-K, B is M-by-N)
752: *
753: * H = I - W**T T W or H**T = I - W**T T**T W
754: *
755: * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
756: * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V
757: *
758: * ---------------------------------------------------------------------------
759: *
760: NP = MIN( L+1, N )
761: KP = MIN( K-L+1, K )
762: *
763: DO J = 1, L
764: DO I = 1, M
765: WORK( I, K-L+J ) = B( I, J )
766: END DO
767: END DO
768: CALL DTRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( KP, 1 ), LDV,
769: $ WORK( 1, KP ), LDWORK )
770: CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B( 1, NP ), LDB,
771: $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
772: CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, V, LDV,
773: $ ZERO, WORK, LDWORK )
774: *
775: DO J = 1, K
776: DO I = 1, M
777: WORK( I, J ) = WORK( I, J ) + A( I, J )
778: END DO
779: END DO
780: *
781: CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
782: $ WORK, LDWORK )
783: *
784: DO J = 1, K
785: DO I = 1, M
786: A( I, J ) = A( I, J ) - WORK( I, J )
787: END DO
788: END DO
789: *
790: CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
791: $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
792: CALL DGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,
793: $ V, LDV, ONE, B, LDB )
794: CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV,
795: $ WORK( 1, KP ), LDWORK )
796: DO J = 1, L
797: DO I = 1, M
798: B( I, J ) = B( I, J ) - WORK( I, K-L+J )
799: END DO
800: END DO
801: *
802: END IF
803: *
804: RETURN
805: *
806: * End of DTPRFB
807: *
808: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtprfb.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack