OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ctprfb.f
Go to the documentation of this file.
1*> \brief \b CTPRFB 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 CTPRFB + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctprfb.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctprfb.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctprfb.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CTPRFB( 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 A( LDA, * ), B( LDB, * ), T( LDT, * ),
30* $ V( LDV, * ), WORK( LDWORK, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> CTPRFB 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 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 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 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 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 COMPLEX 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 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 complexOTHERauxiliary
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 ctprfb( 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 COMPLEX A( LDA, * ), B( LDB, * ), T( LDT, * ),
262 $ v( ldv, * ), work( ldwork, * )
263* ..
264*
265* ==========================================================================
266*
267* .. Parameters ..
268 COMPLEX ONE, ZERO
269 parameter( one = (1.0,0.0), zero = (0.0,0.0) )
270* ..
271* .. Local Scalars ..
272 INTEGER I, J, MP, NP, KP
273 LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW
274* ..
275* .. External Functions ..
276 LOGICAL LSAME
277 EXTERNAL lsame
278* ..
279* .. External Subroutines ..
280 EXTERNAL cgemm, ctrmm
281* ..
282* .. Intrinsic Functions ..
283 INTRINSIC conjg
284* ..
285* .. Executable Statements ..
286*
287* Quick return if possible
288*
289 IF( m.LE.0 .OR. n.LE.0 .OR. k.LE.0 .OR. l.LT.0 ) RETURN
290*
291 IF( lsame( storev, 'C' ) ) THEN
292 column = .true.
293 row = .false.
294 ELSE IF ( lsame( storev, 'R' ) ) THEN
295 column = .false.
296 row = .true.
297 ELSE
298 column = .false.
299 row = .false.
300 END IF
301*
302 IF( lsame( side, 'L' ) ) THEN
303 left = .true.
304 right = .false.
305 ELSE IF( lsame( side, 'R' ) ) THEN
306 left = .false.
307 right = .true.
308 ELSE
309 left = .false.
310 right = .false.
311 END IF
312*
313 IF( lsame( direct, 'F' ) ) THEN
314 forward = .true.
315 backward = .false.
316 ELSE IF( lsame( direct, 'B' ) ) THEN
317 forward = .false.
318 backward = .true.
319 ELSE
320 forward = .false.
321 backward = .false.
322 END IF
323*
324* ---------------------------------------------------------------------------
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**H C where C = [ A ] (K-by-N)
334* [ B ] (M-by-N)
335*
336* H = I - W T W**H or H**H = I - W T**H W**H
337*
338* A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
339* B = B - V T (A + V**H B) or B = B - V T**H (A + V**H 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 ctrmm( 'L', 'U', 'C', 'N', l, n, one, v( mp, 1 ), ldv,
352 $ work, ldwork )
353 CALL cgemm( 'C', 'N', l, n, m-l, one, v, ldv, b, ldb,
354 $ one, work, ldwork )
355 CALL cgemm( 'C', '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 ctrmm( '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 cgemm( 'N', 'N', m-l, n, k, -one, v, ldv, work, ldwork,
374 $ one, b, ldb )
375 CALL cgemm( 'N', 'N', l, n, k-l, -one, v( mp, kp ), ldv,
376 $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
377 CALL ctrmm( '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**H where C = [ A B ] (A is M-by-K, B is M-by-N)
395*
396* H = I - W T W**H or H**H = I - W T**H W**H
397*
398* A = A - (A + B V) T or A = A - (A + B V) T**H
399* B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
400*
401* ---------------------------------------------------------------------------
402*
403 np = min( n-l+1, n )
404 kp = min( l+1, k )
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 ctrmm( 'R', 'U', 'N', 'N', m, l, one, v( np, 1 ), ldv,
412 $ work, ldwork )
413 CALL cgemm( 'N', 'n', M, L, N-L, ONE, B, LDB,
414 $ V, LDV, ONE, WORK, LDWORK )
415 CALL CGEMM( '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 CTRMM( '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 CGEMM( 'n', 'c', M, N-L, K, -ONE, WORK, LDWORK,
434 $ V, LDV, ONE, B, LDB )
435 CALL CGEMM( 'n', 'c', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
436 $ V( NP, KP ), LDV, ONE, B( 1, NP ), LDB )
437 CALL CTRMM( 'r', 'u', 'c', 'n', M, L, ONE, V( NP, 1 ), LDV,
438 $ WORK, LDWORK )
439 DO J = 1, L
440 DO I = 1, M
441 B( I, N-L+J ) = B( I, N-L+J ) - WORK( I, J )
442 END DO
443 END DO
444*
445* ---------------------------------------------------------------------------
446*
447.AND..AND. ELSE IF( COLUMN BACKWARD LEFT ) THEN
448*
449* ---------------------------------------------------------------------------
450*
451* Let W = [ V ] (M-by-K)
452* [ I ] (K-by-K)
453*
454* Form H C or H**H C where C = [ B ] (M-by-N)
455* [ A ] (K-by-N)
456*
457* H = I - W T W**H or H**H = I - W T**H W**H
458*
459* A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
460* B = B - V T (A + V**H B) or B = B - V T**H (A + V**H 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 CTRMM( 'l', 'l', 'c', 'n', L, N, ONE, V( 1, KP ), LDV,
474 $ WORK( KP, 1 ), LDWORK )
475 CALL CGEMM( 'c', 'n', L, N, M-L, ONE, V( MP, KP ), LDV,
476 $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
477 CALL CGEMM( 'c', '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 CTRMM( '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 CGEMM( 'n', 'n', M-L, N, K, -ONE, V( MP, 1 ), LDV,
496 $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
497 CALL CGEMM( 'n', 'n', L, N, K-L, -ONE, V, LDV,
498 $ WORK, LDWORK, ONE, B, LDB )
499 CALL CTRMM( '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.AND..AND. ELSE IF( COLUMN BACKWARD RIGHT ) THEN
510*
511* ---------------------------------------------------------------------------
512*
513* Let W = [ V ] (N-by-K)
514* [ I ] (K-by-K)
515*
516* Form C H or C H**H where C = [ B A ] (B is M-by-N, A is M-by-K)
517*
518* H = I - W T W**H or H**H = I - W T**H W**H
519*
520* A = A - (A + B V) T or A = A - (A + B V) T**H
521* B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
522*
523* ---------------------------------------------------------------------------
524*
525 NP = MIN( L+1, N )
526 KP = MIN( K-L+1, K )
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 CTRMM( 'r', 'l', 'n', 'n', M, L, ONE, V( 1, KP ), LDV,
534 $ WORK( 1, KP ), LDWORK )
535 CALL CGEMM( 'n', 'n', M, L, N-L, ONE, B( 1, NP ), LDB,
536 $ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK )
537 CALL CGEMM( '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 CTRMM( '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 CGEMM( 'n', 'c', M, N-L, K, -ONE, WORK, LDWORK,
556 $ V( NP, 1 ), LDV, ONE, B( 1, NP ), LDB )
557 CALL CGEMM( 'n', 'c', M, L, K-L, -ONE, WORK, LDWORK,
558 $ V, LDV, ONE, B, LDB )
559 CALL CTRMM( 'r', 'l', 'c', 'n', M, L, ONE, V( 1, KP ), LDV,
560 $ WORK( 1, KP ), LDWORK )
561 DO J = 1, L
562 DO I = 1, M
563 B( I, J ) = B( I, J ) - WORK( I, K-L+J )
564 END DO
565 END DO
566*
567* ---------------------------------------------------------------------------
568*
569.AND..AND. ELSE IF( ROW FORWARD LEFT ) THEN
570*
571* ---------------------------------------------------------------------------
572*
573* Let W = [ I V ] ( I is K-by-K, V is K-by-M )
574*
575* Form H C or H**H C where C = [ A ] (K-by-N)
576* [ B ] (M-by-N)
577*
578* H = I - W**H T W or H**H = I - W**H T**H W
579*
580* A = A - T (A + V B) or A = A - T**H (A + V B)
581* B = B - V**H T (A + V B) or B = B - V**H T**H (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 CTRMM( 'l', 'l', 'n', 'n', L, N, ONE, V( 1, MP ), LDV,
594 $ WORK, LDB )
595 CALL CGEMM( 'n', 'n', L, N, M-L, ONE, V, LDV,B, LDB,
596 $ ONE, WORK, LDWORK )
597 CALL CGEMM( '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 CTRMM( '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 CGEMM( 'c', 'n', M-L, N, K, -ONE, V, LDV, WORK, LDWORK,
616 $ ONE, B, LDB )
617 CALL CGEMM( 'c', 'n', L, N, K-L, -ONE, V( KP, MP ), LDV,
618 $ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB )
619 CALL CTRMM( 'l', 'l', 'c', 'n', L, N, ONE, V( 1, MP ), LDV,
620 $ WORK, LDWORK )
621 DO J = 1, N
622 DO I = 1, L
623 B( M-L+I, J ) = B( M-L+I, J ) - WORK( I, J )
624 END DO
625 END DO
626*
627* ---------------------------------------------------------------------------
628*
629.AND..AND. ELSE IF( ROW FORWARD RIGHT ) THEN
630*
631* ---------------------------------------------------------------------------
632*
633* Let W = [ I V ] ( I is K-by-K, V is K-by-N )
634*
635* Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N)
636*
637* H = I - W**H T W or H**H = I - W**H T**H W
638*
639* A = A - (A + B V**H) T or A = A - (A + B V**H) T**H
640* B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V
641*
642* ---------------------------------------------------------------------------
643*
644 NP = MIN( N-L+1, N )
645 KP = MIN( L+1, K )
646*
647 DO J = 1, L
648 DO I = 1, M
649 WORK( I, J ) = B( I, N-L+J )
650 END DO
651 END DO
652 CALL CTRMM( 'r', 'l', 'c', 'n', M, L, ONE, V( 1, NP ), LDV,
653 $ WORK, LDWORK )
654 CALL CGEMM( 'n', 'c', M, L, N-L, ONE, B, LDB, V, LDV,
655 $ ONE, WORK, LDWORK )
656 CALL CGEMM( 'n', 'c', 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 CTRMM( '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 CGEMM( 'n', 'n', M, N-L, K, -ONE, WORK, LDWORK,
675 $ V, LDV, ONE, B, LDB )
676 CALL CGEMM( 'n', 'n', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK,
677 $ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB )
678 CALL CTRMM( '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.AND..AND. ELSE IF( ROW BACKWARD LEFT ) THEN
689*
690* ---------------------------------------------------------------------------
691*
692* Let W = [ V I ] ( I is K-by-K, V is K-by-M )
693*
694* Form H C or H**H C where C = [ B ] (M-by-N)
695* [ A ] (K-by-N)
696*
697* H = I - W**H T W or H**H = I - W**H T**H W
698*
699* A = A - T (A + V B) or A = A - T**H (A + V B)
700* B = B - V**H T (A + V B) or B = B - V**H T**H (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 CTRMM( 'l', 'u', 'n', 'n', L, N, ONE, V( KP, 1 ), LDV,
713 $ WORK( KP, 1 ), LDWORK )
714 CALL CGEMM( 'n', 'n', L, N, M-L, ONE, V( KP, MP ), LDV,
715 $ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK )
716 CALL CGEMM( '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 CTRMM( '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 CGEMM( 'c', 'n', M-L, N, K, -ONE, V( 1, MP ), LDV,
735 $ WORK, LDWORK, ONE, B( MP, 1 ), LDB )
736 CALL CGEMM( 'c', 'n', L, N, K-L, -ONE, V, LDV,
737 $ WORK, LDWORK, ONE, B, LDB )
738 CALL CTRMM( 'l', 'u', 'c', '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.AND..AND. ELSE IF( ROW BACKWARD RIGHT ) THEN
749*
750* ---------------------------------------------------------------------------
751*
752* Let W = [ V I ] ( I is K-by-K, V is K-by-N )
753*
754* Form C H or C H**H where C = [ B A ] (A is M-by-K, B is M-by-N)
755*
756* H = I - W**H T W or H**H = I - W**H T**H W
757*
758* A = A - (A + B V**H) T or A = A - (A + B V**H) T**H
759* B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V
760*
761* ---------------------------------------------------------------------------
762*
763 NP = MIN( L+1, N )
764 KP = MIN( K-L+1, K )
765*
766 DO J = 1, L
767 DO I = 1, M
768 WORK( I, K-L+J ) = B( I, J )
769 END DO
770 END DO
771 CALL CTRMM( 'r', 'u', 'c', 'n', M, L, ONE, V( KP, 1 ), LDV,
772 $ WORK( 1, KP ), LDWORK )
773 CALL CGEMM( 'n', 'c', M, L, N-L, ONE, B( 1, NP ), LDB,
774 $ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK )
775 CALL CGEMM( 'n', 'c', 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 CTRMM( '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 CGEMM( 'n', 'n', M, N-L, K, -ONE, WORK, LDWORK,
794 $ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB )
795 CALL CGEMM( 'n', 'n', M, L, K-L , -ONE, WORK, LDWORK,
796 $ V, LDV, ONE, B, LDB )
797 CALL CTRMM( '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 CTPRFB
810*
811 END
ctprfb
subroutine ctprfb(side, trans, direct, storev, m, n, k, l, v, ldv, t, ldt, a, lda, b, ldb, work, ldwork)
CTPRFB applies a real or complex "triangular-pentagonal" blocked reflector to a real or complex matri...
Definition ctprfb.f:251
ctrmm
subroutine ctrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRMM
Definition ctrmm.f:177
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:187
min
#define min(a, b)
Definition macros.h:20