LAPACK 3.12.0
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 "triangular-pentagonal" 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 249 of file stprfb.f.

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*
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: