LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 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 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 *> \date September 2012
198 *
199 *> \ingroup complexOTHERauxiliary
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 ctprfb( 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  COMPLEX A( lda, * ), B( ldb, * ), T( ldt, * ),
265  $ v( ldv, * ), work( ldwork, * )
266 * ..
267 *
268 * ==========================================================================
269 *
270 * .. Parameters ..
271  COMPLEX ONE, ZERO
272  parameter( one = (1.0,0.0), zero = (0.0,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 cgemm, ctrmm
284 * ..
285 * .. Intrinsic Functions ..
286  INTRINSIC conjg
287 * ..
288 * .. Executable Statements ..
289 *
290 * Quick return if possible
291 *
292  IF( m.LE.0 .OR. n.LE.0 .OR. k.LE.0 .OR. l.LT.0 ) RETURN
293 *
294  IF( lsame( storev, 'C' ) ) THEN
295  column = .true.
296  row = .false.
297  ELSE IF ( lsame( storev, 'R' ) ) THEN
298  column = .false.
299  row = .true.
300  ELSE
301  column = .false.
302  row = .false.
303  END IF
304 *
305  IF( lsame( side, 'L' ) ) THEN
306  left = .true.
307  right = .false.
308  ELSE IF( lsame( side, 'R' ) ) THEN
309  left = .false.
310  right = .true.
311  ELSE
312  left = .false.
313  right = .false.
314  END IF
315 *
316  IF( lsame( direct, 'F' ) ) THEN
317  forward = .true.
318  backward = .false.
319  ELSE IF( lsame( direct, 'B' ) ) THEN
320  forward = .false.
321  backward = .true.
322  ELSE
323  forward = .false.
324  backward = .false.
325  END IF
326 *
327 * ---------------------------------------------------------------------------
328 *
329  IF( column .AND. forward .AND. left ) THEN
330 *
331 * ---------------------------------------------------------------------------
332 *
333 * Let W = [ I ] (K-by-K)
334 * [ V ] (M-by-K)
335 *
336 * Form H C or H**H C where C = [ A ] (K-by-N)
337 * [ B ] (M-by-N)
338 *
339 * H = I - W T W**H or H**H = I - W T**H W**H
340 *
341 * A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
342 * B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
343 *
344 * ---------------------------------------------------------------------------
345 *
346  mp = min( m-l+1, m )
347  kp = min( l+1, k )
348 *
349  DO j = 1, n
350  DO i = 1, l
351  work( i, j ) = b( m-l+i, j )
352  END DO
353  END DO
354  CALL ctrmm( 'L', 'U', 'C', 'N', l, n, one, v( mp, 1 ), ldv,
355  $ work, ldwork )
356  CALL cgemm( 'C', 'N', l, n, m-l, one, v, ldv, b, ldb,
357  $ one, work, ldwork )
358  CALL cgemm( 'C', 'N', k-l, n, m, one, v( 1, kp ), ldv,
359  $ b, ldb, zero, work( kp, 1 ), ldwork )
360 *
361  DO j = 1, n
362  DO i = 1, k
363  work( i, j ) = work( i, j ) + a( i, j )
364  END DO
365  END DO
366 *
367  CALL ctrmm( 'L', 'U', trans, 'N', k, n, one, t, ldt,
368  $ work, ldwork )
369 *
370  DO j = 1, n
371  DO i = 1, k
372  a( i, j ) = a( i, j ) - work( i, j )
373  END DO
374  END DO
375 *
376  CALL cgemm( 'N', 'N', m-l, n, k, -one, v, ldv, work, ldwork,
377  $ one, b, ldb )
378  CALL cgemm( 'N', 'N', l, n, k-l, -one, v( mp, kp ), ldv,
379  $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
380  CALL ctrmm( 'L', 'U', 'N', 'N', l, n, one, v( mp, 1 ), ldv,
381  $ work, ldwork )
382  DO j = 1, n
383  DO i = 1, l
384  b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
385  END DO
386  END DO
387 *
388 * ---------------------------------------------------------------------------
389 *
390  ELSE IF( column .AND. forward .AND. right ) THEN
391 *
392 * ---------------------------------------------------------------------------
393 *
394 * Let W = [ I ] (K-by-K)
395 * [ V ] (N-by-K)
396 *
397 * Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N)
398 *
399 * H = I - W T W**H or H**H = I - W T**H W**H
400 *
401 * A = A - (A + B V) T or A = A - (A + B V) T**H
402 * B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
403 *
404 * ---------------------------------------------------------------------------
405 *
406  np = min( n-l+1, n )
407  kp = min( l+1, k )
408 *
409  DO j = 1, l
410  DO i = 1, m
411  work( i, j ) = b( i, n-l+j )
412  END DO
413  END DO
414  CALL ctrmm( 'R', 'U', 'N', 'N', m, l, one, v( np, 1 ), ldv,
415  $ work, ldwork )
416  CALL cgemm( 'N', 'N', m, l, n-l, one, b, ldb,
417  $ v, ldv, one, work, ldwork )
418  CALL cgemm( 'N', 'N', m, k-l, n, one, b, ldb,
419  $ v( 1, kp ), ldv, zero, work( 1, kp ), ldwork )
420 *
421  DO j = 1, k
422  DO i = 1, m
423  work( i, j ) = work( i, j ) + a( i, j )
424  END DO
425  END DO
426 *
427  CALL ctrmm( 'R', 'U', trans, 'N', m, k, one, t, ldt,
428  $ work, ldwork )
429 *
430  DO j = 1, k
431  DO i = 1, m
432  a( i, j ) = a( i, j ) - work( i, j )
433  END DO
434  END DO
435 *
436  CALL cgemm( 'N', 'C', m, n-l, k, -one, work, ldwork,
437  $ v, ldv, one, b, ldb )
438  CALL cgemm( 'N', 'C', m, l, k-l, -one, work( 1, kp ), ldwork,
439  $ v( np, kp ), ldv, one, b( 1, np ), ldb )
440  CALL ctrmm( 'R', 'U', 'C', 'N', m, l, one, v( np, 1 ), ldv,
441  $ work, ldwork )
442  DO j = 1, l
443  DO i = 1, m
444  b( i, n-l+j ) = b( i, n-l+j ) - work( i, j )
445  END DO
446  END DO
447 *
448 * ---------------------------------------------------------------------------
449 *
450  ELSE IF( column .AND. backward .AND. left ) THEN
451 *
452 * ---------------------------------------------------------------------------
453 *
454 * Let W = [ V ] (M-by-K)
455 * [ I ] (K-by-K)
456 *
457 * Form H C or H**H C where C = [ B ] (M-by-N)
458 * [ A ] (K-by-N)
459 *
460 * H = I - W T W**H or H**H = I - W T**H W**H
461 *
462 * A = A - T (A + V**H B) or A = A - T**H (A + V**H B)
463 * B = B - V T (A + V**H B) or B = B - V T**H (A + V**H B)
464 *
465 * ---------------------------------------------------------------------------
466 *
467  mp = min( l+1, m )
468  kp = min( k-l+1, k )
469 *
470  DO j = 1, n
471  DO i = 1, l
472  work( k-l+i, j ) = b( i, j )
473  END DO
474  END DO
475 *
476  CALL ctrmm( 'L', 'L', 'C', 'N', l, n, one, v( 1, kp ), ldv,
477  $ work( kp, 1 ), ldwork )
478  CALL cgemm( 'C', 'N', l, n, m-l, one, v( mp, kp ), ldv,
479  $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
480  CALL cgemm( 'C', 'N', k-l, n, m, one, v, ldv,
481  $ b, ldb, zero, work, ldwork )
482 *
483  DO j = 1, n
484  DO i = 1, k
485  work( i, j ) = work( i, j ) + a( i, j )
486  END DO
487  END DO
488 *
489  CALL ctrmm( 'L', 'L', trans, 'N', k, n, one, t, ldt,
490  $ work, ldwork )
491 *
492  DO j = 1, n
493  DO i = 1, k
494  a( i, j ) = a( i, j ) - work( i, j )
495  END DO
496  END DO
497 *
498  CALL cgemm( 'N', 'N', m-l, n, k, -one, v( mp, 1 ), ldv,
499  $ work, ldwork, one, b( mp, 1 ), ldb )
500  CALL cgemm( 'N', 'N', l, n, k-l, -one, v, ldv,
501  $ work, ldwork, one, b, ldb )
502  CALL ctrmm( 'L', 'L', 'N', 'N', l, n, one, v( 1, kp ), ldv,
503  $ work( kp, 1 ), ldwork )
504  DO j = 1, n
505  DO i = 1, l
506  b( i, j ) = b( i, j ) - work( k-l+i, j )
507  END DO
508  END DO
509 *
510 * ---------------------------------------------------------------------------
511 *
512  ELSE IF( column .AND. backward .AND. right ) THEN
513 *
514 * ---------------------------------------------------------------------------
515 *
516 * Let W = [ V ] (N-by-K)
517 * [ I ] (K-by-K)
518 *
519 * Form C H or C H**H where C = [ B A ] (B is M-by-N, A is M-by-K)
520 *
521 * H = I - W T W**H or H**H = I - W T**H W**H
522 *
523 * A = A - (A + B V) T or A = A - (A + B V) T**H
524 * B = B - (A + B V) T V**H or B = B - (A + B V) T**H V**H
525 *
526 * ---------------------------------------------------------------------------
527 *
528  np = min( l+1, n )
529  kp = min( k-l+1, k )
530 *
531  DO j = 1, l
532  DO i = 1, m
533  work( i, k-l+j ) = b( i, j )
534  END DO
535  END DO
536  CALL ctrmm( 'R', 'L', 'N', 'N', m, l, one, v( 1, kp ), ldv,
537  $ work( 1, kp ), ldwork )
538  CALL cgemm( 'N', 'N', m, l, n-l, one, b( 1, np ), ldb,
539  $ v( np, kp ), ldv, one, work( 1, kp ), ldwork )
540  CALL cgemm( 'N', 'N', m, k-l, n, one, b, ldb,
541  $ v, ldv, zero, work, ldwork )
542 *
543  DO j = 1, k
544  DO i = 1, m
545  work( i, j ) = work( i, j ) + a( i, j )
546  END DO
547  END DO
548 *
549  CALL ctrmm( 'R', 'L', trans, 'N', m, k, one, t, ldt,
550  $ work, ldwork )
551 *
552  DO j = 1, k
553  DO i = 1, m
554  a( i, j ) = a( i, j ) - work( i, j )
555  END DO
556  END DO
557 *
558  CALL cgemm( 'N', 'C', m, n-l, k, -one, work, ldwork,
559  $ v( np, 1 ), ldv, one, b( 1, np ), ldb )
560  CALL cgemm( 'N', 'C', m, l, k-l, -one, work, ldwork,
561  $ v, ldv, one, b, ldb )
562  CALL ctrmm( 'R', 'L', 'C', 'N', m, l, one, v( 1, kp ), ldv,
563  $ work( 1, kp ), ldwork )
564  DO j = 1, l
565  DO i = 1, m
566  b( i, j ) = b( i, j ) - work( i, k-l+j )
567  END DO
568  END DO
569 *
570 * ---------------------------------------------------------------------------
571 *
572  ELSE IF( row .AND. forward .AND. left ) THEN
573 *
574 * ---------------------------------------------------------------------------
575 *
576 * Let W = [ I V ] ( I is K-by-K, V is K-by-M )
577 *
578 * Form H C or H**H C where C = [ A ] (K-by-N)
579 * [ B ] (M-by-N)
580 *
581 * H = I - W**H T W or H**H = I - W**H T**H W
582 *
583 * A = A - T (A + V B) or A = A - T**H (A + V B)
584 * B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B)
585 *
586 * ---------------------------------------------------------------------------
587 *
588  mp = min( m-l+1, m )
589  kp = min( l+1, k )
590 *
591  DO j = 1, n
592  DO i = 1, l
593  work( i, j ) = b( m-l+i, j )
594  END DO
595  END DO
596  CALL ctrmm( 'L', 'L', 'N', 'N', l, n, one, v( 1, mp ), ldv,
597  $ work, ldb )
598  CALL cgemm( 'N', 'N', l, n, m-l, one, v, ldv,b, ldb,
599  $ one, work, ldwork )
600  CALL cgemm( 'N', 'N', k-l, n, m, one, v( kp, 1 ), ldv,
601  $ b, ldb, zero, work( kp, 1 ), ldwork )
602 *
603  DO j = 1, n
604  DO i = 1, k
605  work( i, j ) = work( i, j ) + a( i, j )
606  END DO
607  END DO
608 *
609  CALL ctrmm( 'L', 'U', trans, 'N', k, n, one, t, ldt,
610  $ work, ldwork )
611 *
612  DO j = 1, n
613  DO i = 1, k
614  a( i, j ) = a( i, j ) - work( i, j )
615  END DO
616  END DO
617 *
618  CALL cgemm( 'C', 'N', m-l, n, k, -one, v, ldv, work, ldwork,
619  $ one, b, ldb )
620  CALL cgemm( 'C', 'N', l, n, k-l, -one, v( kp, mp ), ldv,
621  $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
622  CALL ctrmm( 'L', 'L', 'C', 'N', l, n, one, v( 1, mp ), ldv,
623  $ work, ldwork )
624  DO j = 1, n
625  DO i = 1, l
626  b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
627  END DO
628  END DO
629 *
630 * ---------------------------------------------------------------------------
631 *
632  ELSE IF( row .AND. forward .AND. right ) THEN
633 *
634 * ---------------------------------------------------------------------------
635 *
636 * Let W = [ I V ] ( I is K-by-K, V is K-by-N )
637 *
638 * Form C H or C H**H where C = [ A B ] (A is M-by-K, B is M-by-N)
639 *
640 * H = I - W**H T W or H**H = I - W**H T**H W
641 *
642 * A = A - (A + B V**H) T or A = A - (A + B V**H) T**H
643 * B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V
644 *
645 * ---------------------------------------------------------------------------
646 *
647  np = min( n-l+1, n )
648  kp = min( l+1, k )
649 *
650  DO j = 1, l
651  DO i = 1, m
652  work( i, j ) = b( i, n-l+j )
653  END DO
654  END DO
655  CALL ctrmm( 'R', 'L', 'C', 'N', m, l, one, v( 1, np ), ldv,
656  $ work, ldwork )
657  CALL cgemm( 'N', 'C', m, l, n-l, one, b, ldb, v, ldv,
658  $ one, work, ldwork )
659  CALL cgemm( 'N', 'C', m, k-l, n, one, b, ldb,
660  $ v( kp, 1 ), ldv, zero, work( 1, kp ), ldwork )
661 *
662  DO j = 1, k
663  DO i = 1, m
664  work( i, j ) = work( i, j ) + a( i, j )
665  END DO
666  END DO
667 *
668  CALL ctrmm( 'R', 'U', trans, 'N', m, k, one, t, ldt,
669  $ work, ldwork )
670 *
671  DO j = 1, k
672  DO i = 1, m
673  a( i, j ) = a( i, j ) - work( i, j )
674  END DO
675  END DO
676 *
677  CALL cgemm( 'N', 'N', m, n-l, k, -one, work, ldwork,
678  $ v, ldv, one, b, ldb )
679  CALL cgemm( 'N', 'N', m, l, k-l, -one, work( 1, kp ), ldwork,
680  $ v( kp, np ), ldv, one, b( 1, np ), ldb )
681  CALL ctrmm( 'R', 'L', 'N', 'N', m, l, one, v( 1, np ), ldv,
682  $ work, ldwork )
683  DO j = 1, l
684  DO i = 1, m
685  b( i, n-l+j ) = b( i, n-l+j ) - work( i, j )
686  END DO
687  END DO
688 *
689 * ---------------------------------------------------------------------------
690 *
691  ELSE IF( row .AND. backward .AND. left ) THEN
692 *
693 * ---------------------------------------------------------------------------
694 *
695 * Let W = [ V I ] ( I is K-by-K, V is K-by-M )
696 *
697 * Form H C or H**H C where C = [ B ] (M-by-N)
698 * [ A ] (K-by-N)
699 *
700 * H = I - W**H T W or H**H = I - W**H T**H W
701 *
702 * A = A - T (A + V B) or A = A - T**H (A + V B)
703 * B = B - V**H T (A + V B) or B = B - V**H T**H (A + V B)
704 *
705 * ---------------------------------------------------------------------------
706 *
707  mp = min( l+1, m )
708  kp = min( k-l+1, k )
709 *
710  DO j = 1, n
711  DO i = 1, l
712  work( k-l+i, j ) = b( i, j )
713  END DO
714  END DO
715  CALL ctrmm( 'L', 'U', 'N', 'N', l, n, one, v( kp, 1 ), ldv,
716  $ work( kp, 1 ), ldwork )
717  CALL cgemm( 'N', 'N', l, n, m-l, one, v( kp, mp ), ldv,
718  $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
719  CALL cgemm( 'N', 'N', k-l, n, m, one, v, ldv, b, ldb,
720  $ zero, work, ldwork )
721 *
722  DO j = 1, n
723  DO i = 1, k
724  work( i, j ) = work( i, j ) + a( i, j )
725  END DO
726  END DO
727 *
728  CALL ctrmm( 'L', 'L ', trans, 'N', k, n, one, t, ldt,
729  $ work, ldwork )
730 *
731  DO j = 1, n
732  DO i = 1, k
733  a( i, j ) = a( i, j ) - work( i, j )
734  END DO
735  END DO
736 *
737  CALL cgemm( 'C', 'N', m-l, n, k, -one, v( 1, mp ), ldv,
738  $ work, ldwork, one, b( mp, 1 ), ldb )
739  CALL cgemm( 'C', 'N', l, n, k-l, -one, v, ldv,
740  $ work, ldwork, one, b, ldb )
741  CALL ctrmm( 'L', 'U', 'C', 'N', l, n, one, v( kp, 1 ), ldv,
742  $ work( kp, 1 ), ldwork )
743  DO j = 1, n
744  DO i = 1, l
745  b( i, j ) = b( i, j ) - work( k-l+i, j )
746  END DO
747  END DO
748 *
749 * ---------------------------------------------------------------------------
750 *
751  ELSE IF( row .AND. backward .AND. right ) THEN
752 *
753 * ---------------------------------------------------------------------------
754 *
755 * Let W = [ V I ] ( I is K-by-K, V is K-by-N )
756 *
757 * Form C H or C H**H where C = [ B A ] (A is M-by-K, B is M-by-N)
758 *
759 * H = I - W**H T W or H**H = I - W**H T**H W
760 *
761 * A = A - (A + B V**H) T or A = A - (A + B V**H) T**H
762 * B = B - (A + B V**H) T V or B = B - (A + B V**H) T**H V
763 *
764 * ---------------------------------------------------------------------------
765 *
766  np = min( l+1, n )
767  kp = min( k-l+1, k )
768 *
769  DO j = 1, l
770  DO i = 1, m
771  work( i, k-l+j ) = b( i, j )
772  END DO
773  END DO
774  CALL ctrmm( 'R', 'U', 'C', 'N', m, l, one, v( kp, 1 ), ldv,
775  $ work( 1, kp ), ldwork )
776  CALL cgemm( 'N', 'C', m, l, n-l, one, b( 1, np ), ldb,
777  $ v( kp, np ), ldv, one, work( 1, kp ), ldwork )
778  CALL cgemm( 'N', 'C', m, k-l, n, one, b, ldb, v, ldv,
779  $ zero, work, ldwork )
780 *
781  DO j = 1, k
782  DO i = 1, m
783  work( i, j ) = work( i, j ) + a( i, j )
784  END DO
785  END DO
786 *
787  CALL ctrmm( 'R', 'L', trans, 'N', m, k, one, t, ldt,
788  $ work, ldwork )
789 *
790  DO j = 1, k
791  DO i = 1, m
792  a( i, j ) = a( i, j ) - work( i, j )
793  END DO
794  END DO
795 *
796  CALL cgemm( 'N', 'N', m, n-l, k, -one, work, ldwork,
797  $ v( 1, np ), ldv, one, b( 1, np ), ldb )
798  CALL cgemm( 'N', 'N', m, l, k-l , -one, work, ldwork,
799  $ v, ldv, one, b, ldb )
800  CALL ctrmm( 'R', 'U', 'N', 'N', m, l, one, v( kp, 1 ), ldv,
801  $ work( 1, kp ), ldwork )
802  DO j = 1, l
803  DO i = 1, m
804  b( i, j ) = b( i, j ) - work( i, k-l+j )
805  END DO
806  END DO
807 *
808  END IF
809 *
810  RETURN
811 *
812 * End of CTPRFB
813 *
814  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:253
ctrmm
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:179
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189