251 SUBROUTINE dtprfb( SIDE, TRANS, DIRECT, STOREV, M, N, K, L,
252 $ v, ldv, t, ldt, a, lda, b, ldb, work, ldwork )
260 CHARACTER DIRECT, SIDE, STOREV, TRANS
261 INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N
264 DOUBLE PRECISION A( lda, * ), B( ldb, * ), T( ldt, * ),
265 $ v( ldv, * ), work( ldwork, * )
271 DOUBLE PRECISION ONE, ZERO
272 parameter( one = 1.0, zero = 0.0 )
275 INTEGER I, J, MP, NP, KP
276 LOGICAL LEFT, FORWARD, COLUMN, RIGHT, BACKWARD, ROW
289 IF( m.LE.0 .OR. n.LE.0 .OR. k.LE.0 .OR. l.LT.0 )
RETURN
291 IF( lsame( storev,
'C' ) )
THEN
294 ELSE IF ( lsame( storev,
'R' ) )
THEN
302 IF( lsame( side,
'L' ) )
THEN
305 ELSE IF( lsame( side,
'R' ) )
THEN
313 IF( lsame( direct,
'F' ) )
THEN
316 ELSE IF( lsame( direct,
'B' ) )
THEN
326 IF( column .AND. forward .AND. left )
THEN
348 work( i, j ) = b( m-l+i, j )
351 CALL
dtrmm(
'L',
'U',
'T',
'N', l, n, one, v( mp, 1 ), ldv,
353 CALL
dgemm(
'T',
'N', l, n, m-l, one, v, ldv, b, ldb,
354 $ one, work, ldwork )
355 CALL
dgemm(
'T',
'N', k-l, n, m, one, v( 1, kp ), ldv,
356 $ b, ldb, zero, work( kp, 1 ), ldwork )
360 work( i, j ) = work( i, j ) + a( i, j )
364 CALL
dtrmm(
'L',
'U', trans,
'N', k, n, one, t, ldt,
369 a( i, j ) = a( i, j ) - work( i, j )
373 CALL
dgemm(
'N',
'N', m-l, n, k, -one, v, ldv, work, ldwork,
375 CALL
dgemm(
'N',
'N', l, n, k-l, -one, v( mp, kp ), ldv,
376 $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
377 CALL
dtrmm(
'L',
'U',
'N',
'N', l, n, one, v( mp, 1 ), ldv,
381 b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
387 ELSE IF( column .AND. forward .AND. right )
THEN
408 work( i, j ) = b( i, n-l+j )
411 CALL
dtrmm(
'R',
'U',
'N',
'N', m, l, one, v( np, 1 ), ldv,
413 CALL
dgemm(
'N',
'N', m, l, n-l, one, b, ldb,
414 $ v, ldv, one, work, ldwork )
415 CALL
dgemm(
'N',
'N', m, k-l, n, one, b, ldb,
416 $ v( 1, kp ), ldv, zero, work( 1, kp ), ldwork )
420 work( i, j ) = work( i, j ) + a( i, j )
424 CALL
dtrmm(
'R',
'U', trans,
'N', m, k, one, t, ldt,
429 a( i, j ) = a( i, j ) - work( i, j )
433 CALL
dgemm(
'N',
'T', m, n-l, k, -one, work, ldwork,
434 $ v, ldv, one, b, ldb )
435 CALL
dgemm(
'N',
'T', m, l, k-l, -one, work( 1, kp ), ldwork,
436 $ v( np, kp ), ldv, one, b( 1, np ), ldb )
437 CALL
dtrmm(
'R',
'U',
'T',
'N', m, l, one, v( np, 1 ), ldv,
441 b( i, n-l+j ) = b( i, n-l+j ) - work( i, j )
447 ELSE IF( column .AND. backward .AND. left )
THEN
469 work( k-l+i, j ) = b( i, j )
473 CALL
dtrmm(
'L',
'L',
'T',
'N', l, n, one, v( 1, kp ), ldv,
474 $ work( kp, 1 ), ldwork )
475 CALL
dgemm(
'T',
'N', l, n, m-l, one, v( mp, kp ), ldv,
476 $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
477 CALL
dgemm(
'T',
'N', k-l, n, m, one, v, ldv,
478 $ b, ldb, zero, work, ldwork )
482 work( i, j ) = work( i, j ) + a( i, j )
486 CALL
dtrmm(
'L',
'L', trans,
'N', k, n, one, t, ldt,
491 a( i, j ) = a( i, j ) - work( i, j )
495 CALL
dgemm(
'N',
'N', m-l, n, k, -one, v( mp, 1 ), ldv,
496 $ work, ldwork, one, b( mp, 1 ), ldb )
497 CALL
dgemm(
'N',
'N', l, n, k-l, -one, v, ldv,
498 $ work, ldwork, one, b, ldb )
499 CALL
dtrmm(
'L',
'L',
'N',
'N', l, n, one, v( 1, kp ), ldv,
500 $ work( kp, 1 ), ldwork )
503 b( i, j ) = b( i, j ) - work( k-l+i, j )
509 ELSE IF( column .AND. backward .AND. right )
THEN
530 work( i, k-l+j ) = b( i, j )
533 CALL
dtrmm(
'R',
'L',
'N',
'N', m, l, one, v( 1, kp ), ldv,
534 $ work( 1, kp ), ldwork )
535 CALL
dgemm(
'N',
'N', m, l, n-l, one, b( 1, np ), ldb,
536 $ v( np, kp ), ldv, one, work( 1, kp ), ldwork )
537 CALL
dgemm(
'N',
'N', m, k-l, n, one, b, ldb,
538 $ v, ldv, zero, work, ldwork )
542 work( i, j ) = work( i, j ) + a( i, j )
546 CALL
dtrmm(
'R',
'L', trans,
'N', m, k, one, t, ldt,
551 a( i, j ) = a( i, j ) - work( i, j )
555 CALL
dgemm(
'N',
'T', m, n-l, k, -one, work, ldwork,
556 $ v( np, 1 ), ldv, one, b( 1, np ), ldb )
557 CALL
dgemm(
'N',
'T', m, l, k-l, -one, work, ldwork,
558 $ v, ldv, one, b, ldb )
559 CALL
dtrmm(
'R',
'L',
'T',
'N', m, l, one, v( 1, kp ), ldv,
560 $ work( 1, kp ), ldwork )
563 b( i, j ) = b( i, j ) - work( i, k-l+j )
569 ELSE IF( row .AND. forward .AND. left )
THEN
590 work( i, j ) = b( m-l+i, j )
593 CALL
dtrmm(
'L',
'L',
'N',
'N', l, n, one, v( 1, mp ), ldv,
595 CALL
dgemm(
'N',
'N', l, n, m-l, one, v, ldv,b, ldb,
596 $ one, work, ldwork )
597 CALL
dgemm(
'N',
'N', k-l, n, m, one, v( kp, 1 ), ldv,
598 $ b, ldb, zero, work( kp, 1 ), ldwork )
602 work( i, j ) = work( i, j ) + a( i, j )
606 CALL
dtrmm(
'L',
'U', trans,
'N', k, n, one, t, ldt,
611 a( i, j ) = a( i, j ) - work( i, j )
615 CALL
dgemm(
'T',
'N', m-l, n, k, -one, v, ldv, work, ldwork,
617 CALL
dgemm(
'T',
'N', l, n, k-l, -one, v( kp, mp ), ldv,
618 $ work( kp, 1 ), ldwork, one, b( mp, 1 ), ldb )
619 CALL
dtrmm(
'L',
'L',
'T',
'N', l, n, one, v( 1, mp ), ldv,
623 b( m-l+i, j ) = b( m-l+i, j ) - work( i, j )
629 ELSE IF( row .AND. forward .AND. right )
THEN
649 work( i, j ) = b( i, n-l+j )
652 CALL
dtrmm(
'R',
'L',
'T',
'N', m, l, one, v( 1, np ), ldv,
654 CALL
dgemm(
'N',
'T', m, l, n-l, one, b, ldb, v, ldv,
655 $ one, work, ldwork )
656 CALL
dgemm(
'N',
'T', m, k-l, n, one, b, ldb,
657 $ v( kp, 1 ), ldv, zero, work( 1, kp ), ldwork )
661 work( i, j ) = work( i, j ) + a( i, j )
665 CALL
dtrmm(
'R',
'U', trans,
'N', m, k, one, t, ldt,
670 a( i, j ) = a( i, j ) - work( i, j )
674 CALL
dgemm(
'N',
'N', m, n-l, k, -one, work, ldwork,
675 $ v, ldv, one, b, ldb )
676 CALL
dgemm(
'N',
'N', m, l, k-l, -one, work( 1, kp ), ldwork,
677 $ v( kp, np ), ldv, one, b( 1, np ), ldb )
678 CALL
dtrmm(
'R',
'L',
'N',
'N', m, l, one, v( 1, np ), ldv,
682 b( i, n-l+j ) = b( i, n-l+j ) - work( i, j )
688 ELSE IF( row .AND. backward .AND. left )
THEN
709 work( k-l+i, j ) = b( i, j )
712 CALL
dtrmm(
'L',
'U',
'N',
'N', l, n, one, v( kp, 1 ), ldv,
713 $ work( kp, 1 ), ldwork )
714 CALL
dgemm(
'N',
'N', l, n, m-l, one, v( kp, mp ), ldv,
715 $ b( mp, 1 ), ldb, one, work( kp, 1 ), ldwork )
716 CALL
dgemm(
'N',
'N', k-l, n, m, one, v, ldv, b, ldb,
717 $ zero, work, ldwork )
721 work( i, j ) = work( i, j ) + a( i, j )
725 CALL
dtrmm(
'L',
'L ', trans,
'N', k, n, one, t, ldt,
730 a( i, j ) = a( i, j ) - work( i, j )
734 CALL
dgemm(
'T',
'N', m-l, n, k, -one, v( 1, mp ), ldv,
735 $ work, ldwork, one, b( mp, 1 ), ldb )
736 CALL
dgemm(
'T',
'N', l, n, k-l, -one, v, ldv,
737 $ work, ldwork, one, b, ldb )
738 CALL
dtrmm(
'L',
'U',
'T',
'N', l, n, one, v( kp, 1 ), ldv,
739 $ work( kp, 1 ), ldwork )
742 b( i, j ) = b( i, j ) - work( k-l+i, j )
748 ELSE IF( row .AND. backward .AND. right )
THEN
768 work( i, k-l+j ) = b( i, j )
771 CALL
dtrmm(
'R',
'U',
'T',
'N', m, l, one, v( kp, 1 ), ldv,
772 $ work( 1, kp ), ldwork )
773 CALL
dgemm(
'N',
'T', m, l, n-l, one, b( 1, np ), ldb,
774 $ v( kp, np ), ldv, one, work( 1, kp ), ldwork )
775 CALL
dgemm(
'N',
'T', m, k-l, n, one, b, ldb, v, ldv,
776 $ zero, work, ldwork )
780 work( i, j ) = work( i, j ) + a( i, j )
784 CALL
dtrmm(
'R',
'L', trans,
'N', m, k, one, t, ldt,
789 a( i, j ) = a( i, j ) - work( i, j )
793 CALL
dgemm(
'N',
'N', m, n-l, k, -one, work, ldwork,
794 $ v( 1, np ), ldv, one, b( 1, np ), ldb )
795 CALL
dgemm(
'N',
'N', m, l, k-l , -one, work, ldwork,
796 $ v, ldv, one, b, ldb )
797 CALL
dtrmm(
'R',
'U',
'N',
'N', m, l, one, v( kp, 1 ), ldv,
798 $ work( 1, kp ), ldwork )
801 b( i, j ) = b( i, j ) - work( i, k-l+j )
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dtprfb(SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK)
DTPRFB applies a real or complex "triangular-pentagonal" blocked reflector to a real or complex matri...
subroutine dtrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRMM