LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 Futher Details.
149 *> \endverbatim
150 *>
151 *> \param[in] LDA
152 *> \verbatim
153 *> LDA is INTEGER
154 *> The leading dimension of the array A.
155 *> If SIDE = 'L', LDC >= max(1,K);
156 *> If SIDE = 'R', LDC >= max(1,M).
157 *> \endverbatim
158 *>
159 *> \param[in,out] B
160 *> \verbatim
161 *> B is DOUBLE PRECISION array, dimension (LDB,N)
162 *> On entry, the M-by-N matrix B.
163 *> On exit, B is overwritten by the corresponding block of
164 *> H*C or H**T*C or C*H or C*H**T. See Further Details.
165 *> \endverbatim
166 *>
167 *> \param[in] LDB
168 *> \verbatim
169 *> LDB is INTEGER
170 *> The leading dimension of the array B.
171 *> LDB >= max(1,M).
172 *> \endverbatim
173 *>
174 *> \param[out] WORK
175 *> \verbatim
176 *> WORK is DOUBLE PRECISION array, dimension
177 *> (LDWORK,N) if SIDE = 'L',
178 *> (LDWORK,K) if SIDE = 'R'.
179 *> \endverbatim
180 *>
181 *> \param[in] LDWORK
182 *> \verbatim
183 *> LDWORK is INTEGER
184 *> The leading dimension of the array WORK.
185 *> If SIDE = 'L', LDWORK >= K;
186 *> if SIDE = 'R', LDWORK >= M.
187 *> \endverbatim
188 *
189 * Authors:
190 * ========
191 *
192 *> \author Univ. of Tennessee
193 *> \author Univ. of California Berkeley
194 *> \author Univ. of Colorado Denver
195 *> \author NAG Ltd.
196 *
197 *> \date September 2012
198 *
199 *> \ingroup doubleOTHERauxiliary
200 *
201 *> \par Further Details:
202 * =====================
203 *>
204 *> \verbatim
205 *>
206 *> The matrix C is a composite matrix formed from blocks A and B.
207 *> The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K,
208 *> and if SIDE = 'L', A is of size K-by-N.
209 *>
210 *> If SIDE = 'R' and DIRECT = 'F', C = [A B].
211 *>
212 *> If SIDE = 'L' and DIRECT = 'F', C = [A]
213 *> [B].
214 *>
215 *> If SIDE = 'R' and DIRECT = 'B', C = [B A].
216 *>
217 *> If SIDE = 'L' and DIRECT = 'B', C = [B]
218 *> [A].
219 *>
220 *> The pentagonal matrix V is composed of a rectangular block V1 and a
221 *> trapezoidal block V2. The size of the trapezoidal block is determined by
222 *> the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular;
223 *> if L=0, there is no trapezoidal block, thus V = V1 is rectangular.
224 *>
225 *> If DIRECT = 'F' and STOREV = 'C': V = [V1]
226 *> [V2]
227 *> - V2 is upper trapezoidal (first L rows of K-by-K upper triangular)
228 *>
229 *> If DIRECT = 'F' and STOREV = 'R': V = [V1 V2]
230 *>
231 *> - V2 is lower trapezoidal (first L columns of K-by-K lower triangular)
232 *>
233 *> If DIRECT = 'B' and STOREV = 'C': V = [V2]
234 *> [V1]
235 *> - V2 is lower trapezoidal (last L rows of K-by-K lower triangular)
236 *>
237 *> If DIRECT = 'B' and STOREV = 'R': V = [V2 V1]
238 *>
239 *> - V2 is upper trapezoidal (last L columns of K-by-K upper triangular)
240 *>
241 *> If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K.
242 *>
243 *> If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K.
244 *>
245 *> If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L.
246 *>
247 *> If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.
248 *> \endverbatim
249 *>
250 * =====================================================================
251  SUBROUTINE dtprfb( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
252  $ v, ldv, t, ldt, a, lda, b, ldb, work, ldwork )
253 *
254 * -- LAPACK auxiliary routine (version 3.4.2) --
255 * -- LAPACK is a software package provided by Univ. of Tennessee, --
256 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
257 * September 2012
258 *
259 * .. Scalar Arguments ..
260  CHARACTER DIRECT, SIDE, STOREV, TRANS
261  INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
262 * ..
263 * .. Array Arguments ..
264  DOUBLE PRECISION A( lda, * ), B( ldb, * ), T( ldt, * ),
265  $ v( ldv, * ), work( ldwork, * )
266 * ..
267 *
268 * ==========================================================================
269 *
270 * .. Parameters ..
271  DOUBLE PRECISION ONE, ZERO
272  parameter( one = 1.0, zero = 0.0 )
273 * ..
274 * .. Local Scalars ..
275  INTEGER I, J, MP, NP, KP
276  LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW
277 * ..
278 * .. External Functions ..
279  LOGICAL LSAME
280  EXTERNAL lsame
281 * ..
282 * .. External Subroutines ..
283  EXTERNAL dgemm, dtrmm
284 * ..
285 * .. Executable Statements ..
286 *
287 * Quick return if possible
288 *
289  IF( m.LE.0 .OR. n.LE.0 .OR. k.LE.0 .OR. l.LT.0 ) RETURN
290 *
291  IF( lsame( storev, 'C' ) ) THEN
292  column = .true.
293  row = .false.
294  ELSE IF ( lsame( storev, 'R' ) ) THEN
295  column = .false.
296  row = .true.
297  ELSE
298  column = .false.
299  row = .false.
300  END IF
301 *
302  IF( lsame( side, 'L' ) ) THEN
303  left = .true.
304  right = .false.
305  ELSE IF( lsame( side, 'R' ) ) THEN
306  left = .false.
307  right = .true.
308  ELSE
309  left = .false.
310  right = .false.
311  END IF
312 *
313  IF( lsame( direct, 'F' ) ) THEN
314  forward = .true.
315  backward = .false.
316  ELSE IF( lsame( direct, 'B' ) ) THEN
317  forward = .false.
318  backward = .true.
319  ELSE
320  forward = .false.
321  backward = .false.
322  END IF
323 *
324 * ---------------------------------------------------------------------------
325 *
326  IF( column .AND. forward .AND. left ) THEN
327 *
328 * ---------------------------------------------------------------------------
329 *
330 * Let W = [ I ] (K-by-K)
331 * [ V ] (M-by-K)
332 *
333 * Form H C or H**T C where C = [ A ] (K-by-N)
334 * [ B ] (M-by-N)
335 *
336 * H = I - W T W**T or H**T = I - W T**T W**T
337 *
338 * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
339 * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
340 *
341 * ---------------------------------------------------------------------------
342 *
343  mp = min( m-l+1, m )
344  kp = min( l+1, k )
345 *
346  DO j = 1, n
347  DO i = 1, l
348  work( i, j ) = b( m-l+i, j )
349  END DO
350  END DO
351  CALL dtrmm( 'L', 'U', 'T', 'N', l, n, one, v( mp, 1 ), ldv,
352  $ work, ldwork )
353  CALL dgemm( 'T', 'N', l, n, m-l, one, v, ldv, b, ldb,
354  $ one, work, ldwork )
355  CALL dgemm( 'T', 'N', k-l, n, m, one, v( 1, kp ), ldv,
356  $ b, ldb, zero, work( kp, 1 ), ldwork )
357 *
358  DO j = 1, n
359  DO i = 1, k
360  work( i, j ) = work( i, j ) + a( i, j )
361  END DO
362  END DO
363 *
364  CALL dtrmm( 'L', 'U', trans, 'N', k, n, one, t, ldt,
365  $ work, ldwork )
366 *
367  DO j = 1, n
368  DO i = 1, k
369  a( i, j ) = a( i, j ) - work( i, j )
370  END DO
371  END DO
372 *
373  CALL dgemm( 'N', 'N', m-l, n, k, -one, v, ldv, work, ldwork,
374  $ one, b, ldb )
375  CALL dgemm( 'N', 'N', l, n, k-l, -one, v( mp, kp ), ldv,
376  $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
377  CALL dtrmm( 'L', 'U', 'N', 'N', l, n, one, v( mp, 1 ), ldv,
378  $ work, ldwork )
379  DO j = 1, n
380  DO i = 1, l
381  b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
382  END DO
383  END DO
384 *
385 * ---------------------------------------------------------------------------
386 *
387  ELSE IF( column .AND. forward .AND. right ) THEN
388 *
389 * ---------------------------------------------------------------------------
390 *
391 * Let W = [ I ] (K-by-K)
392 * [ V ] (N-by-K)
393 *
394 * Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N)
395 *
396 * H = I - W T W**T or H**T = I - W T**T W**T
397 *
398 * A = A - (A + B V) T or A = A - (A + B V) T**T
399 * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
400 *
401 * ---------------------------------------------------------------------------
402 *
403  np = min( n-l+1, n )
404  kp = min( l+1, k )
405 *
406  DO j = 1, l
407  DO i = 1, m
408  work( i, j ) = b( i, n-l+j )
409  END DO
410  END DO
411  CALL dtrmm( 'R', 'U', 'N', 'N', m, l, one, v( np, 1 ), ldv,
412  $ work, ldwork )
413  CALL dgemm( 'N', 'N', m, l, n-l, one, b, ldb,
414  $ v, ldv, one, work, ldwork )
415  CALL dgemm( 'N', 'N', m, k-l, n, one, b, ldb,
416  $ v( 1, kp ), ldv, zero, work( 1, kp ), ldwork )
417 *
418  DO j = 1, k
419  DO i = 1, m
420  work( i, j ) = work( i, j ) + a( i, j )
421  END DO
422  END DO
423 *
424  CALL dtrmm( 'R', 'U', trans, 'N', m, k, one, t, ldt,
425  $ work, ldwork )
426 *
427  DO j = 1, k
428  DO i = 1, m
429  a( i, j ) = a( i, j ) - work( i, j )
430  END DO
431  END DO
432 *
433  CALL dgemm( 'N', 'T', m, n-l, k, -one, work, ldwork,
434  $ v, ldv, one, b, ldb )
435  CALL dgemm( 'N', 'T', m, l, k-l, -one, work( 1, kp ), ldwork,
436  $ v( np, kp ), ldv, one, b( 1, np ), ldb )
437  CALL dtrmm( 'R', 'U', 'T', 'N', m, l, one, v( np, 1 ), ldv,
438  $ work, ldwork )
439  DO j = 1, l
440  DO i = 1, m
441  b( i, n-l+j ) = b( i, n-l+j ) - work( i, j )
442  END DO
443  END DO
444 *
445 * ---------------------------------------------------------------------------
446 *
447  ELSE IF( column .AND. backward .AND. left ) THEN
448 *
449 * ---------------------------------------------------------------------------
450 *
451 * Let W = [ V ] (M-by-K)
452 * [ I ] (K-by-K)
453 *
454 * Form H C or H**T C where C = [ B ] (M-by-N)
455 * [ A ] (K-by-N)
456 *
457 * H = I - W T W**T or H**T = I - W T**T W**T
458 *
459 * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
460 * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
461 *
462 * ---------------------------------------------------------------------------
463 *
464  mp = min( l+1, m )
465  kp = min( k-l+1, k )
466 *
467  DO j = 1, n
468  DO i = 1, l
469  work( k-l+i, j ) = b( i, j )
470  END DO
471  END DO
472 *
473  CALL dtrmm( 'L', 'L', 'T', 'N', l, n, one, v( 1, kp ), ldv,
474  $ work( kp, 1 ), ldwork )
475  CALL dgemm( 'T', 'N', l, n, m-l, one, v( mp, kp ), ldv,
476  $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
477  CALL dgemm( 'T', 'N', k-l, n, m, one, v, ldv,
478  $ b, ldb, zero, work, ldwork )
479 *
480  DO j = 1, n
481  DO i = 1, k
482  work( i, j ) = work( i, j ) + a( i, j )
483  END DO
484  END DO
485 *
486  CALL dtrmm( 'L', 'L', trans, 'N', k, n, one, t, ldt,
487  $ work, ldwork )
488 *
489  DO j = 1, n
490  DO i = 1, k
491  a( i, j ) = a( i, j ) - work( i, j )
492  END DO
493  END DO
494 *
495  CALL dgemm( 'N', 'N', m-l, n, k, -one, v( mp, 1 ), ldv,
496  $ work, ldwork, one, b( mp, 1 ), ldb )
497  CALL dgemm( 'N', 'N', l, n, k-l, -one, v, ldv,
498  $ work, ldwork, one, b, ldb )
499  CALL dtrmm( 'L', 'L', 'N', 'N', l, n, one, v( 1, kp ), ldv,
500  $ work( kp, 1 ), ldwork )
501  DO j = 1, n
502  DO i = 1, l
503  b( i, j ) = b( i, j ) - work( k-l+i, j )
504  END DO
505  END DO
506 *
507 * ---------------------------------------------------------------------------
508 *
509  ELSE IF( column .AND. backward .AND. right ) THEN
510 *
511 * ---------------------------------------------------------------------------
512 *
513 * Let W = [ V ] (N-by-K)
514 * [ I ] (K-by-K)
515 *
516 * Form C H or C H**T where C = [ B A ] (B is M-by-N, A is M-by-K)
517 *
518 * H = I - W T W**T or H**T = I - W T**T W**T
519 *
520 * A = A - (A + B V) T or A = A - (A + B V) T**T
521 * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
522 *
523 * ---------------------------------------------------------------------------
524 *
525  np = min( l+1, n )
526  kp = min( k-l+1, k )
527 *
528  DO j = 1, l
529  DO i = 1, m
530  work( i, k-l+j ) = b( i, j )
531  END DO
532  END DO
533  CALL dtrmm( 'R', 'L', 'N', 'N', m, l, one, v( 1, kp ), ldv,
534  $ work( 1, kp ), ldwork )
535  CALL dgemm( 'N', 'N', m, l, n-l, one, b( 1, np ), ldb,
536  $ v( np, kp ), ldv, one, work( 1, kp ), ldwork )
537  CALL dgemm( 'N', 'N', m, k-l, n, one, b, ldb,
538  $ v, ldv, zero, work, ldwork )
539 *
540  DO j = 1, k
541  DO i = 1, m
542  work( i, j ) = work( i, j ) + a( i, j )
543  END DO
544  END DO
545 *
546  CALL dtrmm( 'R', 'L', trans, 'N', m, k, one, t, ldt,
547  $ work, ldwork )
548 *
549  DO j = 1, k
550  DO i = 1, m
551  a( i, j ) = a( i, j ) - work( i, j )
552  END DO
553  END DO
554 *
555  CALL dgemm( 'N', 'T', m, n-l, k, -one, work, ldwork,
556  $ v( np, 1 ), ldv, one, b( 1, np ), ldb )
557  CALL dgemm( 'N', 'T', m, l, k-l, -one, work, ldwork,
558  $ v, ldv, one, b, ldb )
559  CALL dtrmm( 'R', 'L', 'T', 'N', m, l, one, v( 1, kp ), ldv,
560  $ work( 1, kp ), ldwork )
561  DO j = 1, l
562  DO i = 1, m
563  b( i, j ) = b( i, j ) - work( i, k-l+j )
564  END DO
565  END DO
566 *
567 * ---------------------------------------------------------------------------
568 *
569  ELSE IF( row .AND. forward .AND. left ) THEN
570 *
571 * ---------------------------------------------------------------------------
572 *
573 * Let W = [ I V ] ( I is K-by-K, V is K-by-M )
574 *
575 * Form H C or H**T C where C = [ A ] (K-by-N)
576 * [ B ] (M-by-N)
577 *
578 * H = I - W**T T W or H**T = I - W**T T**T W
579 *
580 * A = A - T (A + V B) or A = A - T**T (A + V B)
581 * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
582 *
583 * ---------------------------------------------------------------------------
584 *
585  mp = min( m-l+1, m )
586  kp = min( l+1, k )
587 *
588  DO j = 1, n
589  DO i = 1, l
590  work( i, j ) = b( m-l+i, j )
591  END DO
592  END DO
593  CALL dtrmm( 'L', 'L', 'N', 'N', l, n, one, v( 1, mp ), ldv,
594  $ work, ldb )
595  CALL dgemm( 'N', 'N', l, n, m-l, one, v, ldv,b, ldb,
596  $ one, work, ldwork )
597  CALL dgemm( 'N', 'N', k-l, n, m, one, v( kp, 1 ), ldv,
598  $ b, ldb, zero, work( kp, 1 ), ldwork )
599 *
600  DO j = 1, n
601  DO i = 1, k
602  work( i, j ) = work( i, j ) + a( i, j )
603  END DO
604  END DO
605 *
606  CALL dtrmm( 'L', 'U', trans, 'N', k, n, one, t, ldt,
607  $ work, ldwork )
608 *
609  DO j = 1, n
610  DO i = 1, k
611  a( i, j ) = a( i, j ) - work( i, j )
612  END DO
613  END DO
614 *
615  CALL dgemm( 'T', 'N', m-l, n, k, -one, v, ldv, work, ldwork,
616  $ one, b, ldb )
617  CALL dgemm( 'T', 'N', l, n, k-l, -one, v( kp, mp ), ldv,
618  $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
619  CALL dtrmm( 'L', 'L', 'T', 'N', l, n, one, v( 1, mp ), ldv,
620  $ work, ldwork )
621  DO j = 1, n
622  DO i = 1, l
623  b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
624  END DO
625  END DO
626 *
627 * ---------------------------------------------------------------------------
628 *
629  ELSE IF( row .AND. forward .AND. right ) THEN
630 *
631 * ---------------------------------------------------------------------------
632 *
633 * Let W = [ I V ] ( I is K-by-K, V is K-by-N )
634 *
635 * Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N)
636 *
637 * H = I - W**T T W or H**T = I - W**T T**T W
638 *
639 * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
640 * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V
641 *
642 * ---------------------------------------------------------------------------
643 *
644  np = min( n-l+1, n )
645  kp = min( l+1, k )
646 *
647  DO j = 1, l
648  DO i = 1, m
649  work( i, j ) = b( i, n-l+j )
650  END DO
651  END DO
652  CALL dtrmm( 'R', 'L', 'T', 'N', m, l, one, v( 1, np ), ldv,
653  $ work, ldwork )
654  CALL dgemm( 'N', 'T', m, l, n-l, one, b, ldb, v, ldv,
655  $ one, work, ldwork )
656  CALL dgemm( 'N', 'T', m, k-l, n, one, b, ldb,
657  $ v( kp, 1 ), ldv, zero, work( 1, kp ), ldwork )
658 *
659  DO j = 1, k
660  DO i = 1, m
661  work( i, j ) = work( i, j ) + a( i, j )
662  END DO
663  END DO
664 *
665  CALL dtrmm( 'R', 'U', trans, 'N', m, k, one, t, ldt,
666  $ work, ldwork )
667 *
668  DO j = 1, k
669  DO i = 1, m
670  a( i, j ) = a( i, j ) - work( i, j )
671  END DO
672  END DO
673 *
674  CALL dgemm( 'N', 'N', m, n-l, k, -one, work, ldwork,
675  $ v, ldv, one, b, ldb )
676  CALL dgemm( 'N', 'N', m, l, k-l, -one, work( 1, kp ), ldwork,
677  $ v( kp, np ), ldv, one, b( 1, np ), ldb )
678  CALL dtrmm( 'R', 'L', 'N', 'N', m, l, one, v( 1, np ), ldv,
679  $ work, ldwork )
680  DO j = 1, l
681  DO i = 1, m
682  b( i, n-l+j ) = b( i, n-l+j ) - work( i, j )
683  END DO
684  END DO
685 *
686 * ---------------------------------------------------------------------------
687 *
688  ELSE IF( row .AND. backward .AND. left ) THEN
689 *
690 * ---------------------------------------------------------------------------
691 *
692 * Let W = [ V I ] ( I is K-by-K, V is K-by-M )
693 *
694 * Form H C or H**T C where C = [ B ] (M-by-N)
695 * [ A ] (K-by-N)
696 *
697 * H = I - W**T T W or H**T = I - W**T T**T W
698 *
699 * A = A - T (A + V B) or A = A - T**T (A + V B)
700 * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
701 *
702 * ---------------------------------------------------------------------------
703 *
704  mp = min( l+1, m )
705  kp = min( k-l+1, k )
706 *
707  DO j = 1, n
708  DO i = 1, l
709  work( k-l+i, j ) = b( i, j )
710  END DO
711  END DO
712  CALL dtrmm( 'L', 'U', 'N', 'N', l, n, one, v( kp, 1 ), ldv,
713  $ work( kp, 1 ), ldwork )
714  CALL dgemm( 'N', 'N', l, n, m-l, one, v( kp, mp ), ldv,
715  $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
716  CALL dgemm( 'N', 'N', k-l, n, m, one, v, ldv, b, ldb,
717  $ zero, work, ldwork )
718 *
719  DO j = 1, n
720  DO i = 1, k
721  work( i, j ) = work( i, j ) + a( i, j )
722  END DO
723  END DO
724 *
725  CALL dtrmm( 'L', 'L ', trans, 'N', k, n, one, t, ldt,
726  $ work, ldwork )
727 *
728  DO j = 1, n
729  DO i = 1, k
730  a( i, j ) = a( i, j ) - work( i, j )
731  END DO
732  END DO
733 *
734  CALL dgemm( 'T', 'N', m-l, n, k, -one, v( 1, mp ), ldv,
735  $ work, ldwork, one, b( mp, 1 ), ldb )
736  CALL dgemm( 'T', 'N', l, n, k-l, -one, v, ldv,
737  $ work, ldwork, one, b, ldb )
738  CALL dtrmm( 'L', 'U', 'T', 'N', l, n, one, v( kp, 1 ), ldv,
739  $ work( kp, 1 ), ldwork )
740  DO j = 1, n
741  DO i = 1, l
742  b( i, j ) = b( i, j ) - work( k-l+i, j )
743  END DO
744  END DO
745 *
746 * ---------------------------------------------------------------------------
747 *
748  ELSE IF( row .AND. backward .AND. right ) THEN
749 *
750 * ---------------------------------------------------------------------------
751 *
752 * Let W = [ V I ] ( I is K-by-K, V is K-by-N )
753 *
754 * Form C H or C H**T where C = [ B A ] (A is M-by-K, B is M-by-N)
755 *
756 * H = I - W**T T W or H**T = I - W**T T**T W
757 *
758 * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
759 * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V
760 *
761 * ---------------------------------------------------------------------------
762 *
763  np = min( l+1, n )
764  kp = min( k-l+1, k )
765 *
766  DO j = 1, l
767  DO i = 1, m
768  work( i, k-l+j ) = b( i, j )
769  END DO
770  END DO
771  CALL dtrmm( 'R', 'U', 'T', 'N', m, l, one, v( kp, 1 ), ldv,
772  $ work( 1, kp ), ldwork )
773  CALL dgemm( 'N', 'T', m, l, n-l, one, b( 1, np ), ldb,
774  $ v( kp, np ), ldv, one, work( 1, kp ), ldwork )
775  CALL dgemm( 'N', 'T', m, k-l, n, one, b, ldb, v, ldv,
776  $ zero, work, ldwork )
777 *
778  DO j = 1, k
779  DO i = 1, m
780  work( i, j ) = work( i, j ) + a( i, j )
781  END DO
782  END DO
783 *
784  CALL dtrmm( 'R', 'L', trans, 'N', m, k, one, t, ldt,
785  $ work, ldwork )
786 *
787  DO j = 1, k
788  DO i = 1, m
789  a( i, j ) = a( i, j ) - work( i, j )
790  END DO
791  END DO
792 *
793  CALL dgemm( 'N', 'N', m, n-l, k, -one, work, ldwork,
794  $ v( 1, np ), ldv, one, b( 1, np ), ldb )
795  CALL dgemm( 'N', 'N', m, l, k-l , -one, work, ldwork,
796  $ v, ldv, one, b, ldb )
797  CALL dtrmm( 'R', 'U', 'N', 'N', m, l, one, v( kp, 1 ), ldv,
798  $ work( 1, kp ), ldwork )
799  DO j = 1, l
800  DO i = 1, m
801  b( i, j ) = b( i, j ) - work( i, k-l+j )
802  END DO
803  END DO
804 *
805  END IF
806 *
807  RETURN
808 *
809 * End of DTPRFB
810 *
811  END
dtrmm
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRMM
Definition: dtrmm.f:179
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
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:253