LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
stprfb.f
1 *> \brief \b STPRFB applies a real "triangular-pentagonal" block reflector to a real 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 STPRFB + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stprfb.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stprfb.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stprfb.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE STPRFB( 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 * REAL A( LDA, * ), B( LDB, * ), T( LDT, * ),
30 * $ V( LDV, * ), WORK( LDWORK, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> STPRFB 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 REAL 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 REAL 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 REAL array, dimension
145 *> (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'
146 *> On entry, the K-by-N or M-by-K matrix A.
147 *> On exit, A is overwritten by the corresponding block of
148 *> H*C or H**T*C or C*H or C*H**T. See Further Details.
149 *> \endverbatim
150 *>
151 *> \param[in] LDA
152 *> \verbatim
153 *> LDA is INTEGER
154 *> The leading dimension of the array A.
155 *> If SIDE = 'L', LDA >= max(1,K);
156 *> If SIDE = 'R', LDA >= max(1,M).
157 *> \endverbatim
158 *>
159 *> \param[in,out] B
160 *> \verbatim
161 *> B is REAL 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 REAL 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 tprfb
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 stprfb( 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  REAL A( lda, * ), B( ldb, * ), T( ldt, * ),
262  $ v( ldv, * ), work( ldwork, * )
263 * ..
264 *
265 * ==========================================================================
266 *
267 * .. Parameters ..
268  REAL ONE, ZERO
269  parameter( one = 1.0, zero = 0.0 )
270 * ..
271 * .. Local Scalars ..
272  INTEGER I, J, MP, NP, KP
273  LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW
274 * ..
275 * .. External Functions ..
276  LOGICAL LSAME
277  EXTERNAL lsame
278 * ..
279 * .. External Subroutines ..
280  EXTERNAL sgemm, strmm
281 * ..
282 * .. Executable Statements ..
283 *
284 * Quick return if possible
285 *
286  IF( m.LE.0 .OR. n.LE.0 .OR. k.LE.0 .OR. l.LT.0 ) RETURN
287 *
288  IF( lsame( storev, 'C' ) ) THEN
289  column = .true.
290  row = .false.
291  ELSE IF ( lsame( storev, 'R' ) ) THEN
292  column = .false.
293  row = .true.
294  ELSE
295  column = .false.
296  row = .false.
297  END IF
298 *
299  IF( lsame( side, 'L' ) ) THEN
300  left = .true.
301  right = .false.
302  ELSE IF( lsame( side, 'R' ) ) THEN
303  left = .false.
304  right = .true.
305  ELSE
306  left = .false.
307  right = .false.
308  END IF
309 *
310  IF( lsame( direct, 'F' ) ) THEN
311  forward = .true.
312  backward = .false.
313  ELSE IF( lsame( direct, 'B' ) ) THEN
314  forward = .false.
315  backward = .true.
316  ELSE
317  forward = .false.
318  backward = .false.
319  END IF
320 *
321 * ---------------------------------------------------------------------------
322 *
323  IF( column .AND. forward .AND. left ) THEN
324 *
325 * ---------------------------------------------------------------------------
326 *
327 * Let W = [ I ] (K-by-K)
328 * [ V ] (M-by-K)
329 *
330 * Form H C or H**T C where C = [ A ] (K-by-N)
331 * [ B ] (M-by-N)
332 *
333 * H = I - W T W**T or H**T = I - W T**T W**T
334 *
335 * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
336 * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
337 *
338 * ---------------------------------------------------------------------------
339 *
340  mp = min( m-l+1, m )
341  kp = min( l+1, k )
342 *
343  DO j = 1, n
344  DO i = 1, l
345  work( i, j ) = b( m-l+i, j )
346  END DO
347  END DO
348  CALL strmm( 'L', 'U', 'T', 'N', l, n, one, v( mp, 1 ), ldv,
349  $ work, ldwork )
350  CALL sgemm( 'T', 'N', l, n, m-l, one, v, ldv, b, ldb,
351  $ one, work, ldwork )
352  CALL sgemm( 'T', 'N', k-l, n, m, one, v( 1, kp ), ldv,
353  $ b, ldb, zero, work( kp, 1 ), ldwork )
354 *
355  DO j = 1, n
356  DO i = 1, k
357  work( i, j ) = work( i, j ) + a( i, j )
358  END DO
359  END DO
360 *
361  CALL strmm( 'L', 'U', trans, 'N', k, n, one, t, ldt,
362  $ work, ldwork )
363 *
364  DO j = 1, n
365  DO i = 1, k
366  a( i, j ) = a( i, j ) - work( i, j )
367  END DO
368  END DO
369 *
370  CALL sgemm( 'N', 'N', m-l, n, k, -one, v, ldv, work, ldwork,
371  $ one, b, ldb )
372  CALL sgemm( 'N', 'N', l, n, k-l, -one, v( mp, kp ), ldv,
373  $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
374  CALL strmm( 'L', 'U', 'N', 'N', l, n, one, v( mp, 1 ), ldv,
375  $ work, ldwork )
376  DO j = 1, n
377  DO i = 1, l
378  b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
379  END DO
380  END DO
381 *
382 * ---------------------------------------------------------------------------
383 *
384  ELSE IF( column .AND. forward .AND. right ) THEN
385 *
386 * ---------------------------------------------------------------------------
387 *
388 * Let W = [ I ] (K-by-K)
389 * [ V ] (N-by-K)
390 *
391 * Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N)
392 *
393 * H = I - W T W**T or H**T = I - W T**T W**T
394 *
395 * A = A - (A + B V) T or A = A - (A + B V) T**T
396 * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
397 *
398 * ---------------------------------------------------------------------------
399 *
400  np = min( n-l+1, n )
401  kp = min( l+1, k )
402 *
403  DO j = 1, l
404  DO i = 1, m
405  work( i, j ) = b( i, n-l+j )
406  END DO
407  END DO
408  CALL strmm( 'R', 'U', 'N', 'N', m, l, one, v( np, 1 ), ldv,
409  $ work, ldwork )
410  CALL sgemm( 'N', 'N', m, l, n-l, one, b, ldb,
411  $ v, ldv, one, work, ldwork )
412  CALL sgemm( 'N', 'N', m, k-l, n, one, b, ldb,
413  $ v( 1, kp ), ldv, zero, work( 1, kp ), ldwork )
414 *
415  DO j = 1, k
416  DO i = 1, m
417  work( i, j ) = work( i, j ) + a( i, j )
418  END DO
419  END DO
420 *
421  CALL strmm( 'R', 'U', trans, 'N', m, k, one, t, ldt,
422  $ work, ldwork )
423 *
424  DO j = 1, k
425  DO i = 1, m
426  a( i, j ) = a( i, j ) - work( i, j )
427  END DO
428  END DO
429 *
430  CALL sgemm( 'N', 'T', m, n-l, k, -one, work, ldwork,
431  $ v, ldv, one, b, ldb )
432  CALL sgemm( 'N', 'T', m, l, k-l, -one, work( 1, kp ), ldwork,
433  $ v( np, kp ), ldv, one, b( 1, np ), ldb )
434  CALL strmm( 'R', 'U', 'T', 'N', m, l, one, v( np, 1 ), ldv,
435  $ work, ldwork )
436  DO j = 1, l
437  DO i = 1, m
438  b( i, n-l+j ) = b( i, n-l+j ) - work( i, j )
439  END DO
440  END DO
441 *
442 * ---------------------------------------------------------------------------
443 *
444  ELSE IF( column .AND. backward .AND. left ) THEN
445 *
446 * ---------------------------------------------------------------------------
447 *
448 * Let W = [ V ] (M-by-K)
449 * [ I ] (K-by-K)
450 *
451 * Form H C or H**T C where C = [ B ] (M-by-N)
452 * [ A ] (K-by-N)
453 *
454 * H = I - W T W**T or H**T = I - W T**T W**T
455 *
456 * A = A - T (A + V**T B) or A = A - T**T (A + V**T B)
457 * B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B)
458 *
459 * ---------------------------------------------------------------------------
460 *
461  mp = min( l+1, m )
462  kp = min( k-l+1, k )
463 *
464  DO j = 1, n
465  DO i = 1, l
466  work( k-l+i, j ) = b( i, j )
467  END DO
468  END DO
469 *
470  CALL strmm( 'L', 'L', 'T', 'N', l, n, one, v( 1, kp ), ldv,
471  $ work( kp, 1 ), ldwork )
472  CALL sgemm( 'T', 'N', l, n, m-l, one, v( mp, kp ), ldv,
473  $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
474  CALL sgemm( 'T', 'N', k-l, n, m, one, v, ldv,
475  $ b, ldb, zero, work, ldwork )
476 *
477  DO j = 1, n
478  DO i = 1, k
479  work( i, j ) = work( i, j ) + a( i, j )
480  END DO
481  END DO
482 *
483  CALL strmm( 'L', 'L', trans, 'N', k, n, one, t, ldt,
484  $ work, ldwork )
485 *
486  DO j = 1, n
487  DO i = 1, k
488  a( i, j ) = a( i, j ) - work( i, j )
489  END DO
490  END DO
491 *
492  CALL sgemm( 'N', 'N', m-l, n, k, -one, v( mp, 1 ), ldv,
493  $ work, ldwork, one, b( mp, 1 ), ldb )
494  CALL sgemm( 'N', 'N', l, n, k-l, -one, v, ldv,
495  $ work, ldwork, one, b, ldb )
496  CALL strmm( 'L', 'L', 'N', 'N', l, n, one, v( 1, kp ), ldv,
497  $ work( kp, 1 ), ldwork )
498  DO j = 1, n
499  DO i = 1, l
500  b( i, j ) = b( i, j ) - work( k-l+i, j )
501  END DO
502  END DO
503 *
504 * ---------------------------------------------------------------------------
505 *
506  ELSE IF( column .AND. backward .AND. right ) THEN
507 *
508 * ---------------------------------------------------------------------------
509 *
510 * Let W = [ V ] (N-by-K)
511 * [ I ] (K-by-K)
512 *
513 * Form C H or C H**T where C = [ B A ] (B is M-by-N, A is M-by-K)
514 *
515 * H = I - W T W**T or H**T = I - W T**T W**T
516 *
517 * A = A - (A + B V) T or A = A - (A + B V) T**T
518 * B = B - (A + B V) T V**T or B = B - (A + B V) T**T V**T
519 *
520 * ---------------------------------------------------------------------------
521 *
522  np = min( l+1, n )
523  kp = min( k-l+1, k )
524 *
525  DO j = 1, l
526  DO i = 1, m
527  work( i, k-l+j ) = b( i, j )
528  END DO
529  END DO
530  CALL strmm( 'R', 'L', 'N', 'N', m, l, one, v( 1, kp ), ldv,
531  $ work( 1, kp ), ldwork )
532  CALL sgemm( 'N', 'N', m, l, n-l, one, b( 1, np ), ldb,
533  $ v( np, kp ), ldv, one, work( 1, kp ), ldwork )
534  CALL sgemm( 'N', 'N', m, k-l, n, one, b, ldb,
535  $ v, ldv, zero, work, ldwork )
536 *
537  DO j = 1, k
538  DO i = 1, m
539  work( i, j ) = work( i, j ) + a( i, j )
540  END DO
541  END DO
542 *
543  CALL strmm( 'R', 'L', trans, 'N', m, k, one, t, ldt,
544  $ work, ldwork )
545 *
546  DO j = 1, k
547  DO i = 1, m
548  a( i, j ) = a( i, j ) - work( i, j )
549  END DO
550  END DO
551 *
552  CALL sgemm( 'N', 'T', m, n-l, k, -one, work, ldwork,
553  $ v( np, 1 ), ldv, one, b( 1, np ), ldb )
554  CALL sgemm( 'N', 'T', m, l, k-l, -one, work, ldwork,
555  $ v, ldv, one, b, ldb )
556  CALL strmm( 'R', 'L', 'T', 'N', m, l, one, v( 1, kp ), ldv,
557  $ work( 1, kp ), ldwork )
558  DO j = 1, l
559  DO i = 1, m
560  b( i, j ) = b( i, j ) - work( i, k-l+j )
561  END DO
562  END DO
563 *
564 * ---------------------------------------------------------------------------
565 *
566  ELSE IF( row .AND. forward .AND. left ) THEN
567 *
568 * ---------------------------------------------------------------------------
569 *
570 * Let W = [ I V ] ( I is K-by-K, V is K-by-M )
571 *
572 * Form H C or H**T C where C = [ A ] (K-by-N)
573 * [ B ] (M-by-N)
574 *
575 * H = I - W**T T W or H**T = I - W**T T**T W
576 *
577 * A = A - T (A + V B) or A = A - T**T (A + V B)
578 * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
579 *
580 * ---------------------------------------------------------------------------
581 *
582  mp = min( m-l+1, m )
583  kp = min( l+1, k )
584 *
585  DO j = 1, n
586  DO i = 1, l
587  work( i, j ) = b( m-l+i, j )
588  END DO
589  END DO
590  CALL strmm( 'L', 'L', 'N', 'N', l, n, one, v( 1, mp ), ldv,
591  $ work, ldb )
592  CALL sgemm( 'N', 'N', l, n, m-l, one, v, ldv,b, ldb,
593  $ one, work, ldwork )
594  CALL sgemm( 'N', 'N', k-l, n, m, one, v( kp, 1 ), ldv,
595  $ b, ldb, zero, work( kp, 1 ), ldwork )
596 *
597  DO j = 1, n
598  DO i = 1, k
599  work( i, j ) = work( i, j ) + a( i, j )
600  END DO
601  END DO
602 *
603  CALL strmm( 'L', 'U', trans, 'N', k, n, one, t, ldt,
604  $ work, ldwork )
605 *
606  DO j = 1, n
607  DO i = 1, k
608  a( i, j ) = a( i, j ) - work( i, j )
609  END DO
610  END DO
611 *
612  CALL sgemm( 'T', 'N', m-l, n, k, -one, v, ldv, work, ldwork,
613  $ one, b, ldb )
614  CALL sgemm( 'T', 'N', l, n, k-l, -one, v( kp, mp ), ldv,
615  $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
616  CALL strmm( 'L', 'L', 'T', 'N', l, n, one, v( 1, mp ), ldv,
617  $ work, ldwork )
618  DO j = 1, n
619  DO i = 1, l
620  b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
621  END DO
622  END DO
623 *
624 * ---------------------------------------------------------------------------
625 *
626  ELSE IF( row .AND. forward .AND. right ) THEN
627 *
628 * ---------------------------------------------------------------------------
629 *
630 * Let W = [ I V ] ( I is K-by-K, V is K-by-N )
631 *
632 * Form C H or C H**T where C = [ A B ] (A is M-by-K, B is M-by-N)
633 *
634 * H = I - W**T T W or H**T = I - W**T T**T W
635 *
636 * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
637 * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V
638 *
639 * ---------------------------------------------------------------------------
640 *
641  np = min( n-l+1, n )
642  kp = min( l+1, k )
643 *
644  DO j = 1, l
645  DO i = 1, m
646  work( i, j ) = b( i, n-l+j )
647  END DO
648  END DO
649  CALL strmm( 'R', 'L', 'T', 'N', m, l, one, v( 1, np ), ldv,
650  $ work, ldwork )
651  CALL sgemm( 'N', 'T', m, l, n-l, one, b, ldb, v, ldv,
652  $ one, work, ldwork )
653  CALL sgemm( 'N', 'T', m, k-l, n, one, b, ldb,
654  $ v( kp, 1 ), ldv, zero, work( 1, kp ), ldwork )
655 *
656  DO j = 1, k
657  DO i = 1, m
658  work( i, j ) = work( i, j ) + a( i, j )
659  END DO
660  END DO
661 *
662  CALL strmm( 'R', 'U', trans, 'N', m, k, one, t, ldt,
663  $ work, ldwork )
664 *
665  DO j = 1, k
666  DO i = 1, m
667  a( i, j ) = a( i, j ) - work( i, j )
668  END DO
669  END DO
670 *
671  CALL sgemm( 'N', 'N', m, n-l, k, -one, work, ldwork,
672  $ v, ldv, one, b, ldb )
673  CALL sgemm( 'N', 'N', m, l, k-l, -one, work( 1, kp ), ldwork,
674  $ v( kp, np ), ldv, one, b( 1, np ), ldb )
675  CALL strmm( 'R', 'L', 'N', 'N', m, l, one, v( 1, np ), ldv,
676  $ work, ldwork )
677  DO j = 1, l
678  DO i = 1, m
679  b( i, n-l+j ) = b( i, n-l+j ) - work( i, j )
680  END DO
681  END DO
682 *
683 * ---------------------------------------------------------------------------
684 *
685  ELSE IF( row .AND. backward .AND. left ) THEN
686 *
687 * ---------------------------------------------------------------------------
688 *
689 * Let W = [ V I ] ( I is K-by-K, V is K-by-M )
690 *
691 * Form H C or H**T C where C = [ B ] (M-by-N)
692 * [ A ] (K-by-N)
693 *
694 * H = I - W**T T W or H**T = I - W**T T**T W
695 *
696 * A = A - T (A + V B) or A = A - T**T (A + V B)
697 * B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B)
698 *
699 * ---------------------------------------------------------------------------
700 *
701  mp = min( l+1, m )
702  kp = min( k-l+1, k )
703 *
704  DO j = 1, n
705  DO i = 1, l
706  work( k-l+i, j ) = b( i, j )
707  END DO
708  END DO
709  CALL strmm( 'L', 'U', 'N', 'N', l, n, one, v( kp, 1 ), ldv,
710  $ work( kp, 1 ), ldwork )
711  CALL sgemm( 'N', 'N', l, n, m-l, one, v( kp, mp ), ldv,
712  $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
713  CALL sgemm( 'N', 'N', k-l, n, m, one, v, ldv, b, ldb,
714  $ zero, work, ldwork )
715 *
716  DO j = 1, n
717  DO i = 1, k
718  work( i, j ) = work( i, j ) + a( i, j )
719  END DO
720  END DO
721 *
722  CALL strmm( 'L', 'L ', trans, 'N', k, n, one, t, ldt,
723  $ work, ldwork )
724 *
725  DO j = 1, n
726  DO i = 1, k
727  a( i, j ) = a( i, j ) - work( i, j )
728  END DO
729  END DO
730 *
731  CALL sgemm( 'T', 'N', m-l, n, k, -one, v( 1, mp ), ldv,
732  $ work, ldwork, one, b( mp, 1 ), ldb )
733  CALL sgemm( 'T', 'N', l, n, k-l, -one, v, ldv,
734  $ work, ldwork, one, b, ldb )
735  CALL strmm( 'L', 'U', 'T', 'N', l, n, one, v( kp, 1 ), ldv,
736  $ work( kp, 1 ), ldwork )
737  DO j = 1, n
738  DO i = 1, l
739  b( i, j ) = b( i, j ) - work( k-l+i, j )
740  END DO
741  END DO
742 *
743 * ---------------------------------------------------------------------------
744 *
745  ELSE IF( row .AND. backward .AND. right ) THEN
746 *
747 * ---------------------------------------------------------------------------
748 *
749 * Let W = [ V I ] ( I is K-by-K, V is K-by-N )
750 *
751 * Form C H or C H**T where C = [ B A ] (A is M-by-K, B is M-by-N)
752 *
753 * H = I - W**T T W or H**T = I - W**T T**T W
754 *
755 * A = A - (A + B V**T) T or A = A - (A + B V**T) T**T
756 * B = B - (A + B V**T) T V or B = B - (A + B V**T) T**T V
757 *
758 * ---------------------------------------------------------------------------
759 *
760  np = min( l+1, n )
761  kp = min( k-l+1, k )
762 *
763  DO j = 1, l
764  DO i = 1, m
765  work( i, k-l+j ) = b( i, j )
766  END DO
767  END DO
768  CALL strmm( 'R', 'U', 'T', 'N', m, l, one, v( kp, 1 ), ldv,
769  $ work( 1, kp ), ldwork )
770  CALL sgemm( 'N', 'T', m, l, n-l, one, b( 1, np ), ldb,
771  $ v( kp, np ), ldv, one, work( 1, kp ), ldwork )
772  CALL sgemm( 'N', 'T', m, k-l, n, one, b, ldb, v, ldv,
773  $ zero, work, ldwork )
774 *
775  DO j = 1, k
776  DO i = 1, m
777  work( i, j ) = work( i, j ) + a( i, j )
778  END DO
779  END DO
780 *
781  CALL strmm( 'R', 'L', trans, 'N', m, k, one, t, ldt,
782  $ work, ldwork )
783 *
784  DO j = 1, k
785  DO i = 1, m
786  a( i, j ) = a( i, j ) - work( i, j )
787  END DO
788  END DO
789 *
790  CALL sgemm( 'N', 'N', m, n-l, k, -one, work, ldwork,
791  $ v( 1, np ), ldv, one, b( 1, np ), ldb )
792  CALL sgemm( 'N', 'N', m, l, k-l , -one, work, ldwork,
793  $ v, ldv, one, b, ldb )
794  CALL strmm( 'R', 'U', 'N', 'N', m, l, one, v( kp, 1 ), ldv,
795  $ work( 1, kp ), ldwork )
796  DO j = 1, l
797  DO i = 1, m
798  b( i, j ) = b( i, j ) - work( i, k-l+j )
799  END DO
800  END DO
801 *
802  END IF
803 *
804  RETURN
805 *
806 * End of STPRFB
807 *
808  END
strmm
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
Definition: strmm.f:177
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
stprfb
subroutine stprfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK)
STPRFB applies a real "triangular-pentagonal" block reflector to a real matrix, which is composed of ...
Definition: stprfb.f:251