OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dtprfb.f
Go to the documentation of this file.
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 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.AND..AND. ELSE IF( COLUMN BACKWARD 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.AND..AND. ELSE IF( COLUMN BACKWARD 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
dtprfb
subroutine dtprfb(side, trans, direct, storev, m, n, k, l, v, ldv, t, ldt, a, lda, b, ldb, work, ldwork)
DTPRFB applies a real or complex "triangular-pentagonal" blocked reflector to a real or complex matri...
Definition dtprfb.f:251
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187
dtrmm
subroutine dtrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRMM
Definition dtrmm.f:177
min
#define min(a, b)
Definition macros.h:20