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