LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ctprfb.f
Go to the documentation of this file.
1*> \brief \b CTPRFB applies a complex "triangular-pentagonal" block reflector to a 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*> Download CTPRFB + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctprfb.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctprfb.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctprfb.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
20* V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
21*
22* .. Scalar Arguments ..
23* CHARACTER DIRECT, SIDE, STOREV, TRANS
24* INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
25* ..
26* .. Array Arguments ..
27* COMPLEX A( LDA, * ), B( LDB, * ), T( LDT, * ),
28* $ V( LDV, * ), WORK( LDWORK, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> CTPRFB applies a complex "triangular-pentagonal" block reflector H or its
38*> conjugate transpose H**H to a complex matrix C, which is composed of two
39*> blocks A and B, either from the left or right.
40*>
41*> \endverbatim
42*
43* Arguments:
44* ==========
45*
46*> \param[in] SIDE
47*> \verbatim
48*> SIDE is CHARACTER*1
49*> = 'L': apply H or H**H from the Left
50*> = 'R': apply H or H**H from the Right
51*> \endverbatim
52*>
53*> \param[in] TRANS
54*> \verbatim
55*> TRANS is CHARACTER*1
56*> = 'N': apply H (No transpose)
57*> = 'C': apply H**H (Conjugate transpose)
58*> \endverbatim
59*>
60*> \param[in] DIRECT
61*> \verbatim
62*> DIRECT is CHARACTER*1
63*> Indicates how H is formed from a product of elementary
64*> reflectors
65*> = 'F': H = H(1) H(2) . . . H(k) (Forward)
66*> = 'B': H = H(k) . . . H(2) H(1) (Backward)
67*> \endverbatim
68*>
69*> \param[in] STOREV
70*> \verbatim
71*> STOREV is CHARACTER*1
72*> Indicates how the vectors which define the elementary
73*> reflectors are stored:
74*> = 'C': Columns
75*> = 'R': Rows
76*> \endverbatim
77*>
78*> \param[in] M
79*> \verbatim
80*> M is INTEGER
81*> The number of rows of the matrix B.
82*> M >= 0.
83*> \endverbatim
84*>
85*> \param[in] N
86*> \verbatim
87*> N is INTEGER
88*> The number of columns of the matrix B.
89*> N >= 0.
90*> \endverbatim
91*>
92*> \param[in] K
93*> \verbatim
94*> K is INTEGER
95*> The order of the matrix T, i.e. the number of elementary
96*> reflectors whose product defines the block reflector.
97*> K >= 0.
98*> \endverbatim
99*>
100*> \param[in] L
101*> \verbatim
102*> L is INTEGER
103*> The order of the trapezoidal part of V.
104*> K >= L >= 0. See Further Details.
105*> \endverbatim
106*>
107*> \param[in] V
108*> \verbatim
109*> V is COMPLEX array, dimension
110*> (LDV,K) if STOREV = 'C'
111*> (LDV,M) if STOREV = 'R' and SIDE = 'L'
112*> (LDV,N) if STOREV = 'R' and SIDE = 'R'
113*> The pentagonal matrix V, which contains the elementary reflectors
114*> H(1), H(2), ..., H(K). See Further Details.
115*> \endverbatim
116*>
117*> \param[in] LDV
118*> \verbatim
119*> LDV is INTEGER
120*> The leading dimension of the array V.
121*> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
122*> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
123*> if STOREV = 'R', LDV >= K.
124*> \endverbatim
125*>
126*> \param[in] T
127*> \verbatim
128*> T is COMPLEX array, dimension (LDT,K)
129*> The triangular K-by-K matrix T in the representation of the
130*> block reflector.
131*> \endverbatim
132*>
133*> \param[in] LDT
134*> \verbatim
135*> LDT is INTEGER
136*> The leading dimension of the array T.
137*> LDT >= K.
138*> \endverbatim
139*>
140*> \param[in,out] A
141*> \verbatim
142*> A is COMPLEX array, dimension
143*> (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'
144*> On entry, the K-by-N or M-by-K matrix A.
145*> On exit, A is overwritten by the corresponding block of
146*> H*C or H**H*C or C*H or C*H**H. See Further Details.
147*> \endverbatim
148*>
149*> \param[in] LDA
150*> \verbatim
151*> LDA is INTEGER
152*> The leading dimension of the array A.
153*> If SIDE = 'L', LDA >= max(1,K);
154*> If SIDE = 'R', LDA >= max(1,M).
155*> \endverbatim
156*>
157*> \param[in,out] B
158*> \verbatim
159*> B is COMPLEX array, dimension (LDB,N)
160*> On entry, the M-by-N matrix B.
161*> On exit, B is overwritten by the corresponding block of
162*> H*C or H**H*C or C*H or C*H**H. See Further Details.
163*> \endverbatim
164*>
165*> \param[in] LDB
166*> \verbatim
167*> LDB is INTEGER
168*> The leading dimension of the array B.
169*> LDB >= max(1,M).
170*> \endverbatim
171*>
172*> \param[out] WORK
173*> \verbatim
174*> WORK is COMPLEX array, dimension
175*> (LDWORK,N) if SIDE = 'L',
176*> (LDWORK,K) if SIDE = 'R'.
177*> \endverbatim
178*>
179*> \param[in] LDWORK
180*> \verbatim
181*> LDWORK is INTEGER
182*> The leading dimension of the array WORK.
183*> If SIDE = 'L', LDWORK >= K;
184*> if SIDE = 'R', LDWORK >= M.
185*> \endverbatim
186*
187* Authors:
188* ========
189*
190*> \author Univ. of Tennessee
191*> \author Univ. of California Berkeley
192*> \author Univ. of Colorado Denver
193*> \author NAG Ltd.
194*
195*> \ingroup tprfb
196*
197*> \par Further Details:
198* =====================
199*>
200*> \verbatim
201*>
202*> The matrix C is a composite matrix formed from blocks A and B.
203*> The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K,
204*> and if SIDE = 'L', A is of size K-by-N.
205*>
206*> If SIDE = 'R' and DIRECT = 'F', C = [A B].
207*>
208*> If SIDE = 'L' and DIRECT = 'F', C = [A]
209*> [B].
210*>
211*> If SIDE = 'R' and DIRECT = 'B', C = [B A].
212*>
213*> If SIDE = 'L' and DIRECT = 'B', C = [B]
214*> [A].
215*>
216*> The pentagonal matrix V is composed of a rectangular block V1 and a
217*> trapezoidal block V2. The size of the trapezoidal block is determined by
218*> the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular;
219*> if L=0, there is no trapezoidal block, thus V = V1 is rectangular.
220*>
221*> If DIRECT = 'F' and STOREV = 'C': V = [V1]
222*> [V2]
223*> - V2 is upper trapezoidal (first L rows of K-by-K upper triangular)
224*>
225*> If DIRECT = 'F' and STOREV = 'R': V = [V1 V2]
226*>
227*> - V2 is lower trapezoidal (first L columns of K-by-K lower triangular)
228*>
229*> If DIRECT = 'B' and STOREV = 'C': V = [V2]
230*> [V1]
231*> - V2 is lower trapezoidal (last L rows of K-by-K lower triangular)
232*>
233*> If DIRECT = 'B' and STOREV = 'R': V = [V2 V1]
234*>
235*> - V2 is upper trapezoidal (last L columns of K-by-K upper triangular)
236*>
237*> If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K.
238*>
239*> If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K.
240*>
241*> If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L.
242*>
243*> If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.
244*> \endverbatim
245*>
246* =====================================================================
247 SUBROUTINE ctprfb( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
248 $ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK )
249*
250* -- LAPACK auxiliary routine --
251* -- LAPACK is a software package provided by Univ. of Tennessee, --
252* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253*
254* .. Scalar Arguments ..
255 CHARACTER DIRECT, SIDE, STOREV, TRANS
256 INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
257* ..
258* .. Array Arguments ..
259 COMPLEX A( LDA, * ), B( LDB, * ), T( LDT, * ),
260 $ v( ldv, * ), work( ldwork, * )
261* ..
262*
263* ==========================================================================
264*
265* .. Parameters ..
266 COMPLEX ONE, ZERO
267 parameter( one = (1.0,0.0), zero = (0.0,0.0) )
268* ..
269* .. Local Scalars ..
270 INTEGER I, J, MP, NP, KP
271 LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW
272* ..
273* .. External Functions ..
274 LOGICAL LSAME
275 EXTERNAL lsame
276* ..
277* .. External Subroutines ..
278 EXTERNAL cgemm, ctrmm
279* ..
280* .. Intrinsic Functions ..
281 INTRINSIC conjg
282* ..
283* .. Executable Statements ..
284*
285* Quick return if possible
286*
287 IF( m.LE.0 .OR. n.LE.0 .OR. k.LE.0 .OR. l.LT.0 ) RETURN
288*
289 IF( lsame( storev, 'C' ) ) THEN
290 column = .true.
291 row = .false.
292 ELSE IF ( lsame( storev, 'R' ) ) THEN
293 column = .false.
294 row = .true.
295 ELSE
296 column = .false.
297 row = .false.
298 END IF
299*
300 IF( lsame( side, 'L' ) ) THEN
301 left = .true.
302 right = .false.
303 ELSE IF( lsame( side, 'R' ) ) THEN
304 left = .false.
305 right = .true.
306 ELSE
307 left = .false.
308 right = .false.
309 END IF
310*
311 IF( lsame( direct, 'F' ) ) THEN
312 forward = .true.
313 backward = .false.
314 ELSE IF( lsame( direct, 'B' ) ) THEN
315 forward = .false.
316 backward = .true.
317 ELSE
318 forward = .false.
319 backward = .false.
320 END IF
321*
322* ---------------------------------------------------------------------------
323*
324 IF( column .AND. forward .AND. left ) THEN
325*
326* ---------------------------------------------------------------------------
327*
328* Let W = [ I ] (K-by-K)
329* [ V ] (M-by-K)
330*
331* Form H C or H**H C where C = [ A ] (K-by-N)
332* [ B ] (M-by-N)
333*
334* H = I - W T W**H or H**H = I - W T**H W**H
335*
336* A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
337* B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
338*
339* ---------------------------------------------------------------------------
340*
341 mp = min( m-l+1, m )
342 kp = min( l+1, k )
343*
344 DO j = 1, n
345 DO i = 1, l
346 work( i, j ) = b( m-l+i, j )
347 END DO
348 END DO
349 CALL ctrmm( 'L', 'U', 'C', 'N', l, n, one, v( mp, 1 ), ldv,
350 $ work, ldwork )
351 CALL cgemm( 'C', 'N', l, n, m-l, one, v, ldv, b, ldb,
352 $ one, work, ldwork )
353 CALL cgemm( 'C', 'N', k-l, n, m, one, v( 1, kp ), ldv,
354 $ b, ldb, zero, work( kp, 1 ), ldwork )
355*
356 DO j = 1, n
357 DO i = 1, k
358 work( i, j ) = work( i, j ) + a( i, j )
359 END DO
360 END DO
361*
362 CALL ctrmm( 'L', 'U', trans, 'N', k, n, one, t, ldt,
363 $ work, ldwork )
364*
365 DO j = 1, n
366 DO i = 1, k
367 a( i, j ) = a( i, j ) - work( i, j )
368 END DO
369 END DO
370*
371 CALL cgemm( 'N', 'N', m-l, n, k, -one, v, ldv, work, ldwork,
372 $ one, b, ldb )
373 CALL cgemm( 'N', 'N', l, n, k-l, -one, v( mp, kp ), ldv,
374 $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
375 CALL ctrmm( 'L', 'U', 'N', 'N', l, n, one, v( mp, 1 ), ldv,
376 $ work, ldwork )
377 DO j = 1, n
378 DO i = 1, l
379 b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
380 END DO
381 END DO
382*
383* ---------------------------------------------------------------------------
384*
385 ELSE IF( column .AND. forward .AND. right ) THEN
386*
387* ---------------------------------------------------------------------------
388*
389* Let W = [ I ] (K-by-K)
390* [ V ] (N-by-K)
391*
392* Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N)
393*
394* H = I - W T W**H or H**H = I - W T**H W**H
395*
396* A = A - (A + B V) T or A = A - (A + B V) T**H
397* B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
398*
399* ---------------------------------------------------------------------------
400*
401 np = min( n-l+1, n )
402 kp = min( l+1, k )
403*
404 DO j = 1, l
405 DO i = 1, m
406 work( i, j ) = b( i, n-l+j )
407 END DO
408 END DO
409 CALL ctrmm( 'R', 'U', 'N', 'N', m, l, one, v( np, 1 ), ldv,
410 $ work, ldwork )
411 CALL cgemm( 'N', 'N', m, l, n-l, one, b, ldb,
412 $ v, ldv, one, work, ldwork )
413 CALL cgemm( 'N', 'N', m, k-l, n, one, b, ldb,
414 $ v( 1, kp ), ldv, zero, work( 1, kp ), ldwork )
415*
416 DO j = 1, k
417 DO i = 1, m
418 work( i, j ) = work( i, j ) + a( i, j )
419 END DO
420 END DO
421*
422 CALL ctrmm( 'R', 'U', trans, 'N', m, k, one, t, ldt,
423 $ work, ldwork )
424*
425 DO j = 1, k
426 DO i = 1, m
427 a( i, j ) = a( i, j ) - work( i, j )
428 END DO
429 END DO
430*
431 CALL cgemm( 'N', 'C', m, n-l, k, -one, work, ldwork,
432 $ v, ldv, one, b, ldb )
433 CALL cgemm( 'N', 'C', m, l, k-l, -one, work( 1, kp ),
434 $ ldwork,
435 $ v( np, kp ), ldv, one, b( 1, np ), ldb )
436 CALL ctrmm( 'R', 'U', 'C', 'N', m, l, one, v( np, 1 ), ldv,
437 $ work, ldwork )
438 DO j = 1, l
439 DO i = 1, m
440 b( i, n-l+j ) = b( i, n-l+j ) - work( i, j )
441 END DO
442 END DO
443*
444* ---------------------------------------------------------------------------
445*
446 ELSE IF( column .AND. backward .AND. left ) THEN
447*
448* ---------------------------------------------------------------------------
449*
450* Let W = [ V ] (M-by-K)
451* [ I ] (K-by-K)
452*
453* Form H C or H**H C where C = [ B ] (M-by-N)
454* [ A ] (K-by-N)
455*
456* H = I - W T W**H or H**H = I - W T**H W**H
457*
458* A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
459* B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
460*
461* ---------------------------------------------------------------------------
462*
463 mp = min( l+1, m )
464 kp = min( k-l+1, k )
465*
466 DO j = 1, n
467 DO i = 1, l
468 work( k-l+i, j ) = b( i, j )
469 END DO
470 END DO
471*
472 CALL ctrmm( 'L', 'L', 'C', 'N', l, n, one, v( 1, kp ), ldv,
473 $ work( kp, 1 ), ldwork )
474 CALL cgemm( 'C', 'N', l, n, m-l, one, v( mp, kp ), ldv,
475 $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
476 CALL cgemm( 'C', 'N', k-l, n, m, one, v, ldv,
477 $ b, ldb, zero, work, ldwork )
478*
479 DO j = 1, n
480 DO i = 1, k
481 work( i, j ) = work( i, j ) + a( i, j )
482 END DO
483 END DO
484*
485 CALL ctrmm( 'L', 'L', trans, 'N', k, n, one, t, ldt,
486 $ work, ldwork )
487*
488 DO j = 1, n
489 DO i = 1, k
490 a( i, j ) = a( i, j ) - work( i, j )
491 END DO
492 END DO
493*
494 CALL cgemm( 'N', 'N', m-l, n, k, -one, v( mp, 1 ), ldv,
495 $ work, ldwork, one, b( mp, 1 ), ldb )
496 CALL cgemm( 'N', 'N', l, n, k-l, -one, v, ldv,
497 $ work, ldwork, one, b, ldb )
498 CALL ctrmm( 'L', 'L', 'N', 'N', l, n, one, v( 1, kp ), ldv,
499 $ work( kp, 1 ), ldwork )
500 DO j = 1, n
501 DO i = 1, l
502 b( i, j ) = b( i, j ) - work( k-l+i, j )
503 END DO
504 END DO
505*
506* ---------------------------------------------------------------------------
507*
508 ELSE IF( column .AND. backward .AND. right ) THEN
509*
510* ---------------------------------------------------------------------------
511*
512* Let W = [ V ] (N-by-K)
513* [ I ] (K-by-K)
514*
515* Form C H or C H**H where C = [ B A ] (B is M-by-N, A is M-by-K)
516*
517* H = I - W T W**H or H**H = I - W T**H W**H
518*
519* A = A - (A + B V) T or A = A - (A + B V) T**H
520* B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
521*
522* ---------------------------------------------------------------------------
523*
524 np = min( l+1, n )
525 kp = min( k-l+1, k )
526*
527 DO j = 1, l
528 DO i = 1, m
529 work( i, k-l+j ) = b( i, j )
530 END DO
531 END DO
532 CALL ctrmm( 'R', 'L', 'N', 'N', m, l, one, v( 1, kp ), ldv,
533 $ work( 1, kp ), ldwork )
534 CALL cgemm( 'N', 'N', m, l, n-l, one, b( 1, np ), ldb,
535 $ v( np, kp ), ldv, one, work( 1, kp ), ldwork )
536 CALL cgemm( 'N', 'N', m, k-l, n, one, b, ldb,
537 $ v, ldv, zero, work, ldwork )
538*
539 DO j = 1, k
540 DO i = 1, m
541 work( i, j ) = work( i, j ) + a( i, j )
542 END DO
543 END DO
544*
545 CALL ctrmm( 'R', 'L', trans, 'N', m, k, one, t, ldt,
546 $ work, ldwork )
547*
548 DO j = 1, k
549 DO i = 1, m
550 a( i, j ) = a( i, j ) - work( i, j )
551 END DO
552 END DO
553*
554 CALL cgemm( 'N', 'C', m, n-l, k, -one, work, ldwork,
555 $ v( np, 1 ), ldv, one, b( 1, np ), ldb )
556 CALL cgemm( 'N', 'C', m, l, k-l, -one, work, ldwork,
557 $ v, ldv, one, b, ldb )
558 CALL ctrmm( 'R', 'L', 'C', 'N', m, l, one, v( 1, kp ), ldv,
559 $ work( 1, kp ), ldwork )
560 DO j = 1, l
561 DO i = 1, m
562 b( i, j ) = b( i, j ) - work( i, k-l+j )
563 END DO
564 END DO
565*
566* ---------------------------------------------------------------------------
567*
568 ELSE IF( row .AND. forward .AND. left ) THEN
569*
570* ---------------------------------------------------------------------------
571*
572* Let W = [ I V ] ( I is K-by-K, V is K-by-M )
573*
574* Form H C or H**H C where C = [ A ] (K-by-N)
575* [ B ] (M-by-N)
576*
577* H = I - W**H T W or H**H = I - W**H T**H W
578*
579* A = A - T (A + V B) or A = A - T**H (A + V B)
580* B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B)
581*
582* ---------------------------------------------------------------------------
583*
584 mp = min( m-l+1, m )
585 kp = min( l+1, k )
586*
587 DO j = 1, n
588 DO i = 1, l
589 work( i, j ) = b( m-l+i, j )
590 END DO
591 END DO
592 CALL ctrmm( 'L', 'L', 'N', 'N', l, n, one, v( 1, mp ), ldv,
593 $ work, ldb )
594 CALL cgemm( 'N', 'N', l, n, m-l, one, v, ldv,b, ldb,
595 $ one, work, ldwork )
596 CALL cgemm( 'N', 'N', k-l, n, m, one, v( kp, 1 ), ldv,
597 $ b, ldb, zero, work( kp, 1 ), ldwork )
598*
599 DO j = 1, n
600 DO i = 1, k
601 work( i, j ) = work( i, j ) + a( i, j )
602 END DO
603 END DO
604*
605 CALL ctrmm( 'L', 'U', trans, 'N', k, n, one, t, ldt,
606 $ work, ldwork )
607*
608 DO j = 1, n
609 DO i = 1, k
610 a( i, j ) = a( i, j ) - work( i, j )
611 END DO
612 END DO
613*
614 CALL cgemm( 'C', 'N', m-l, n, k, -one, v, ldv, work, ldwork,
615 $ one, b, ldb )
616 CALL cgemm( 'C', 'N', l, n, k-l, -one, v( kp, mp ), ldv,
617 $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
618 CALL ctrmm( 'L', 'L', 'C', 'N', l, n, one, v( 1, mp ), ldv,
619 $ work, ldwork )
620 DO j = 1, n
621 DO i = 1, l
622 b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
623 END DO
624 END DO
625*
626* ---------------------------------------------------------------------------
627*
628 ELSE IF( row .AND. forward .AND. right ) THEN
629*
630* ---------------------------------------------------------------------------
631*
632* Let W = [ I V ] ( I is K-by-K, V is K-by-N )
633*
634* Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N)
635*
636* H = I - W**H T W or H**H = I - W**H T**H W
637*
638* A = A - (A + B V**H) T or A = A - (A + B V**H) T**H
639* B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V
640*
641* ---------------------------------------------------------------------------
642*
643 np = min( n-l+1, n )
644 kp = min( l+1, k )
645*
646 DO j = 1, l
647 DO i = 1, m
648 work( i, j ) = b( i, n-l+j )
649 END DO
650 END DO
651 CALL ctrmm( 'R', 'L', 'C', 'N', m, l, one, v( 1, np ), ldv,
652 $ work, ldwork )
653 CALL cgemm( 'N', 'C', m, l, n-l, one, b, ldb, v, ldv,
654 $ one, work, ldwork )
655 CALL cgemm( 'N', 'C', m, k-l, n, one, b, ldb,
656 $ v( kp, 1 ), ldv, zero, work( 1, kp ), ldwork )
657*
658 DO j = 1, k
659 DO i = 1, m
660 work( i, j ) = work( i, j ) + a( i, j )
661 END DO
662 END DO
663*
664 CALL ctrmm( 'R', 'U', trans, 'N', m, k, one, t, ldt,
665 $ work, ldwork )
666*
667 DO j = 1, k
668 DO i = 1, m
669 a( i, j ) = a( i, j ) - work( i, j )
670 END DO
671 END DO
672*
673 CALL cgemm( 'N', 'N', m, n-l, k, -one, work, ldwork,
674 $ v, ldv, one, b, ldb )
675 CALL cgemm( 'N', 'N', m, l, k-l, -one, work( 1, kp ),
676 $ 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 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**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 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**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
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
ctprfb
subroutine ctprfb(side, trans, direct, storev, m, n, k, l, v, ldv, t, ldt, a, lda, b, ldb, work, ldwork)
CTPRFB applies a complex "triangular-pentagonal" block reflector to a complex matrix,...
Definition ctprfb.f:249
ctrmm
subroutine ctrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
CTRMM
Definition ctrmm.f:177