LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ stprfb()

subroutine stprfb ( character side,
character trans,
character direct,
character storev,
integer m,
integer n,
integer k,
integer l,
real, dimension( ldv, * ) v,
integer ldv,
real, dimension( ldt, * ) t,
integer ldt,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( ldwork, * ) work,
integer ldwork )

STPRFB applies a real "triangular-pentagonal" block reflector to a real matrix, which is composed of two blocks.

Download STPRFB + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> STPRFB applies a real  block reflector H or its
!> transpose H**T to a real matrix C, which is composed of two
!> blocks A and B, either from the left or right.
!>
!>
Parameters
[in]SIDE
!>          SIDE is CHARACTER*1
!>          = 'L': apply H or H**T from the Left
!>          = 'R': apply H or H**T from the Right
!>
[in]TRANS
!>          TRANS is CHARACTER*1
!>          = 'N': apply H (No transpose)
!>          = 'T': apply H**T (Transpose)
!>
[in]DIRECT
!>          DIRECT is CHARACTER*1
!>          Indicates how H is formed from a product of elementary
!>          reflectors
!>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
!>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
!>
[in]STOREV
!>          STOREV is CHARACTER*1
!>          Indicates how the vectors which define the elementary
!>          reflectors are stored:
!>          = 'C': Columns
!>          = 'R': Rows
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix B.
!>          M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix B.
!>          N >= 0.
!>
[in]K
!>          K is INTEGER
!>          The order of the matrix T, i.e. the number of elementary
!>          reflectors whose product defines the block reflector.
!>          K >= 0.
!>
[in]L
!>          L is INTEGER
!>          The order of the trapezoidal part of V.
!>          K >= L >= 0.  See Further Details.
!>
[in]V
!>          V is REAL array, dimension
!>                                (LDV,K) if STOREV = 'C'
!>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
!>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
!>          The pentagonal matrix V, which contains the elementary reflectors
!>          H(1), H(2), ..., H(K).  See Further Details.
!>
[in]LDV
!>          LDV is INTEGER
!>          The leading dimension of the array V.
!>          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
!>          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
!>          if STOREV = 'R', LDV >= K.
!>
[in]T
!>          T is REAL array, dimension (LDT,K)
!>          The triangular K-by-K matrix T in the representation of the
!>          block reflector.
!>
[in]LDT
!>          LDT is INTEGER
!>          The leading dimension of the array T.
!>          LDT >= K.
!>
[in,out]A
!>          A is REAL array, dimension
!>          (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R'
!>          On entry, the K-by-N or M-by-K matrix A.
!>          On exit, A is overwritten by the corresponding block of
!>          H*C or H**T*C or C*H or C*H**T.  See Further Details.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.
!>          If SIDE = 'L', LDA >= max(1,K);
!>          If SIDE = 'R', LDA >= max(1,M).
!>
[in,out]B
!>          B is REAL array, dimension (LDB,N)
!>          On entry, the M-by-N matrix B.
!>          On exit, B is overwritten by the corresponding block of
!>          H*C or H**T*C or C*H or C*H**T.  See Further Details.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.
!>          LDB >= max(1,M).
!>
[out]WORK
!>          WORK is REAL array, dimension
!>          (LDWORK,N) if SIDE = 'L',
!>          (LDWORK,K) if SIDE = 'R'.
!>
[in]LDWORK
!>          LDWORK is INTEGER
!>          The leading dimension of the array WORK.
!>          If SIDE = 'L', LDWORK >= K;
!>          if SIDE = 'R', LDWORK >= M.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  The matrix C is a composite matrix formed from blocks A and B.
!>  The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K,
!>  and if SIDE = 'L', A is of size K-by-N.
!>
!>  If SIDE = 'R' and DIRECT = 'F', C = [A B].
!>
!>  If SIDE = 'L' and DIRECT = 'F', C = [A]
!>                                      [B].
!>
!>  If SIDE = 'R' and DIRECT = 'B', C = [B A].
!>
!>  If SIDE = 'L' and DIRECT = 'B', C = [B]
!>                                      [A].
!>
!>  The pentagonal matrix V is composed of a rectangular block V1 and a
!>  trapezoidal block V2.  The size of the trapezoidal block is determined by
!>  the parameter L, where 0<=L<=K.  If L=K, the V2 block of V is triangular;
!>  if L=0, there is no trapezoidal block, thus V = V1 is rectangular.
!>
!>  If DIRECT = 'F' and STOREV = 'C':  V = [V1]
!>                                         [V2]
!>     - V2 is upper trapezoidal (first L rows of K-by-K upper triangular)
!>
!>  If DIRECT = 'F' and STOREV = 'R':  V = [V1 V2]
!>
!>     - V2 is lower trapezoidal (first L columns of K-by-K lower triangular)
!>
!>  If DIRECT = 'B' and STOREV = 'C':  V = [V2]
!>                                         [V1]
!>     - V2 is lower trapezoidal (last L rows of K-by-K lower triangular)
!>
!>  If DIRECT = 'B' and STOREV = 'R':  V = [V2 V1]
!>
!>     - V2 is upper trapezoidal (last L columns of K-by-K upper triangular)
!>
!>  If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K.
!>
!>  If STOREV = 'C' and SIDE = 'R', V is N-by-K with V2 L-by-K.
!>
!>  If STOREV = 'R' and SIDE = 'L', V is K-by-M with V2 K-by-L.
!>
!>  If STOREV = 'R' and SIDE = 'R', V is K-by-N with V2 K-by-L.
!>

Definition at line 247 of file stprfb.f.

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*
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
strmm
subroutine strmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRMM
Definition strmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: