LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
stprfb.f
Go to the documentation of this file.
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
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
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
strmm
subroutine strmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRMM
Definition strmm.f:177