File:
[local] /
rpl /
lapack /
lapack /
dtprfb.f
Revision
1.5:
download - view:
text,
annotated -
select for diffs -
revision graph
Mon Jan 27 09:28:29 2014 UTC (10 years, 3 months ago) by
bertrand
Branches:
MAIN
CVS tags:
rpl-4_1_24,
rpl-4_1_23,
rpl-4_1_22,
rpl-4_1_21,
rpl-4_1_20,
rpl-4_1_19,
rpl-4_1_18,
rpl-4_1_17,
HEAD
Cohérence.
1: *> \brief \b DTPRFB applies a real or complex "triangular-pentagonal" blocked reflector to a real or complex 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 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 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: *> \date September 2012
198: *
199: *> \ingroup doubleOTHERauxiliary
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 DTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
252: $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
253: *
254: * -- LAPACK auxiliary routine (version 3.4.2) --
255: * -- LAPACK is a software package provided by Univ. of Tennessee, --
256: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
257: * September 2012
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: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ),
265: $ V( LDV, * ), WORK( LDWORK, * )
266: * ..
267: *
268: * ==========================================================================
269: *
270: * .. Parameters ..
271: DOUBLE PRECISION ONE, ZERO
272: PARAMETER ( ONE = 1.0, ZERO = 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 DGEMM, DTRMM
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: * ---------------------------------------------------------------------------
325: *
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**T C where C = [ A ] (K-by-N)
334: * [ B ] (M-by-N)
335: *
336: * H = I - W T W**T or H**T = I - W T**T W**T
337: *
338: * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
339: * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
340: *
341: * ---------------------------------------------------------------------------
342: *
343: MP = MIN( M-L+1, M )
344: KP = MIN( L+1, K )
345: *
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 DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( MP, 1 ), LDV,
352: $ WORK, LDWORK )
353: CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V, LDV, B, LDB,
354: $ ONE, WORK, LDWORK )
355: CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V( 1, KP ), LDV,
356: $ B, LDB, ZERO, WORK( KP, 1 ), LDWORK )
357: *
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: *
364: CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
365: $ WORK, LDWORK )
366: *
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 DGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
374: $ ONE, B, LDB )
375: CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV,
376: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
377: CALL DTRMM( '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: * ---------------------------------------------------------------------------
386: *
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**T where C = [ A B ] (A is M-by-K, B is M-by-N)
395: *
396: * H = I - W T W**T or H**T = I - W T**T W**T
397: *
398: * A = A - (A + B V) T or A = A - (A + B V) T**T
399: * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
400: *
401: * ---------------------------------------------------------------------------
402: *
403: NP = MIN( N-L+1, N )
404: KP = MIN( L+1, K )
405: *
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 DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV,
412: $ WORK, LDWORK )
413: CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB,
414: $ V, LDV, ONE, WORK, LDWORK )
415: CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
416: $ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK )
417: *
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: *
424: CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
425: $ WORK, LDWORK )
426: *
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 DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
434: $ V, LDV, ONE, B, LDB )
435: CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
436: $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
437: CALL DTRMM( 'R', 'U', 'T', '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: * ---------------------------------------------------------------------------
446: *
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**T C where C = [ B ] (M-by-N)
455: * [ A ] (K-by-N)
456: *
457: * H = I - W T W**T or H**T = I - W T**T W**T
458: *
459: * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
460: * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
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 DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, KP ), LDV,
474: $ WORK( KP, 1 ), LDWORK )
475: CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V( MP, KP ), LDV,
476: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
477: CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V, LDV,
478: $ B, LDB, ZERO, WORK, LDWORK )
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: *
486: CALL DTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT,
487: $ WORK, LDWORK )
488: *
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: *
495: CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV,
496: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
497: CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV,
498: $ WORK, LDWORK, ONE, B, LDB )
499: CALL DTRMM( '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: * ---------------------------------------------------------------------------
508: *
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**T where C = [ B A ] (B is M-by-N, A is M-by-K)
517: *
518: * H = I - W T W**T or H**T = I - W T**T W**T
519: *
520: * A = A - (A + B V) T or A = A - (A + B V) T**T
521: * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
522: *
523: * ---------------------------------------------------------------------------
524: *
525: NP = MIN( L+1, N )
526: KP = MIN( K-L+1, K )
527: *
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 DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV,
534: $ WORK( 1, KP ), LDWORK )
535: CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB,
536: $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
537: CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB,
538: $ V, LDV, ZERO, WORK, LDWORK )
539: *
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: *
546: CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
547: $ WORK, LDWORK )
548: *
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 DGEMM( 'N', 'T', M, N-L, K, -ONE, WORK, LDWORK,
556: $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
557: CALL DGEMM( 'N', 'T', M, L, K-L, -ONE, WORK, LDWORK,
558: $ V, LDV, ONE, B, LDB )
559: CALL DTRMM( 'R', 'L', 'T', '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: * ---------------------------------------------------------------------------
568: *
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**T C where C = [ A ] (K-by-N)
576: * [ B ] (M-by-N)
577: *
578: * H = I - W**T T W or H**T = I - W**T T**T W
579: *
580: * A = A - T (A + V B) or A = A - T**T (A + V B)
581: * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
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
592: END DO
593: CALL DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV,
594: $ WORK, LDB )
595: CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB,
596: $ ONE, WORK, LDWORK )
597: CALL DGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV,
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: *
606: CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT,
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 DGEMM( 'T', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
616: $ ONE, B, LDB )
617: CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV,
618: $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
619: CALL DTRMM( 'L', 'L', 'T', '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: * ---------------------------------------------------------------------------
628: *
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**T where C = [ A B ] (A is M-by-K, B is M-by-N)
636: *
637: * H = I - W**T T W or H**T = I - W**T T**T W
638: *
639: * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
640: * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T 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 DTRMM( 'R', 'L', 'T', 'N', M, L, ONE, V( 1, NP ), LDV,
653: $ WORK, LDWORK )
654: CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B, LDB, V, LDV,
655: $ ONE, WORK, LDWORK )
656: CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB,
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: *
665: CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT,
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: *
674: CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
675: $ V, LDV, ONE, B, LDB )
676: CALL DGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
677: $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )
678: CALL DTRMM( '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: * ---------------------------------------------------------------------------
687: *
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**T C where C = [ B ] (M-by-N)
695: * [ A ] (K-by-N)
696: *
697: * H = I - W**T T W or H**T = I - W**T T**T W
698: *
699: * A = A - T (A + V B) or A = A - T**T (A + V B)
700: * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
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 DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( KP, 1 ), LDV,
713: $ WORK( KP, 1 ), LDWORK )
714: CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V( KP, MP ), LDV,
715: $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
716: CALL DGEMM( '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 DTRMM( '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 DGEMM( 'T', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV,
735: $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
736: CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V, LDV,
737: $ WORK, LDWORK, ONE, B, LDB )
738: CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( KP, 1 ), LDV,
739: $ WORK( KP, 1 ), LDWORK )
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: * ---------------------------------------------------------------------------
747: *
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**T where C = [ B A ] (A is M-by-K, B is M-by-N)
755: *
756: * H = I - W**T T W or H**T = I - W**T T**T W
757: *
758: * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
759: * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T 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 DTRMM( 'R', 'U', 'T', 'N', M, L, ONE, V( KP, 1 ), LDV,
772: $ WORK( 1, KP ), LDWORK )
773: CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B( 1, NP ), LDB,
774: $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
775: CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, V, LDV,
776: $ ZERO, WORK, LDWORK )
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: *
784: CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT,
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: *
793: CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK,
794: $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
795: CALL DGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK,
796: $ V, LDV, ONE, B, LDB )
797: CALL DTRMM( '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 DTPRFB
810: *
811: 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