Chimie Théorique » OpenPath

      SUBROUTINE DLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,

     $                   LDV, T, LDT, C, LDC, WORK, LDWORK )

*

*  -- LAPACK routine (version 3.2) --

*  -- LAPACK is a software package provided by Univ. of Tennessee,    --

*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--

*     November 2006

*

*     .. Scalar Arguments ..

      CHARACTER          DIRECT, SIDE, STOREV, TRANS

      INTEGER            K, L, LDC, LDT, LDV, LDWORK, M, N

*     ..

*     .. Array Arguments ..

      DOUBLE PRECISION   C( LDC, * ), T( LDT, * ), V( LDV, * ),

     $                   WORK( LDWORK, * )

*     ..

*

*  Purpose

*  =======

*

*  DLARZB applies a real block reflector H or its transpose H**T to

*  a real distributed M-by-N  C from the left or the right.

*

*  Currently, only STOREV = 'R' and DIRECT = 'B' are supported.

*

*  Arguments

*  =========

*

*  SIDE    (input) CHARACTER*1

*          = 'L': apply H or H' from the Left

*          = 'R': apply H or H' from the Right

*

*  TRANS   (input) CHARACTER*1

*          = 'N': apply H (No transpose)

*          = 'C': apply H' (Transpose)

*

*  DIRECT  (input) CHARACTER*1

*          Indicates how H is formed from a product of elementary

*          reflectors

*          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)

*          = 'B': H = H(k) . . . H(2) H(1) (Backward)

*

*  STOREV  (input) CHARACTER*1

*          Indicates how the vectors which define the elementary

*          reflectors are stored:

*          = 'C': Columnwise                        (not supported yet)

*          = 'R': Rowwise

*

*  M       (input) INTEGER

*          The number of rows of the matrix C.

*

*  N       (input) INTEGER

*          The number of columns of the matrix C.

*

*  K       (input) INTEGER

*          The order of the matrix T (= the number of elementary

*          reflectors whose product defines the block reflector).

*

*  L       (input) INTEGER

*          The number of columns of the matrix V containing the

*          meaningful part of the Householder reflectors.

*          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.

*

*  V       (input) DOUBLE PRECISION array, dimension (LDV,NV).

*          If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.

*

*  LDV     (input) INTEGER

*          The leading dimension of the array V.

*          If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.

*

*  T       (input) DOUBLE PRECISION array, dimension (LDT,K)

*          The triangular K-by-K matrix T in the representation of the

*          block reflector.

*

*  LDT     (input) INTEGER

*          The leading dimension of the array T. LDT >= K.

*

*  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)

*          On entry, the M-by-N matrix C.

*          On exit, C is overwritten by H*C or H'*C or C*H or C*H'.

*

*  LDC     (input) INTEGER

*          The leading dimension of the array C. LDC >= max(1,M).

*

*  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)

*

*  LDWORK  (input) INTEGER

*          The leading dimension of the array WORK.

*          If SIDE = 'L', LDWORK >= max(1,N);

*          if SIDE = 'R', LDWORK >= max(1,M).

*

*  Further Details

*  ===============

*

*  Based on contributions by

*    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA

*

*  =====================================================================

*

*     .. Parameters ..

      DOUBLE PRECISION   ONE

      PARAMETER          ( ONE = 1.0D+0 )

*     ..

*     .. Local Scalars ..

      CHARACTER          TRANST

      INTEGER            I, INFO, J

*     ..

*     .. External Functions ..

      LOGICAL            LSAME

      EXTERNAL           LSAME

*     ..

*     .. External Subroutines ..

      EXTERNAL           DCOPY, DGEMM, DTRMM, XERBLA

*     ..

*     .. Executable Statements ..

*

*     Quick return if possible

*

      IF( M.LE.0 .OR. N.LE.0 )

     $   RETURN

*

*     Check for currently supported options

*

      INFO = 0

      IF( .NOT.LSAME( DIRECT, 'B' ) ) THEN

         INFO = -3

      ELSE IF( .NOT.LSAME( STOREV, 'R' ) ) THEN

         INFO = -4

      END IF

      IF( INFO.NE.0 ) THEN

         CALL XERBLA( 'DLARZB', -INFO )

         RETURN

      END IF

*

      IF( LSAME( TRANS, 'N' ) ) THEN

         TRANST = 'T'

      ELSE

         TRANST = 'N'

      END IF

*

      IF( LSAME( SIDE, 'L' ) ) THEN

*

*        Form  H * C  or  H' * C

*

*        W( 1:n, 1:k ) = C( 1:k, 1:n )'

*

         DO 10 J = 1, K

            CALL DCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 )

   10    CONTINUE

*

*        W( 1:n, 1:k ) = W( 1:n, 1:k ) + ...

*                        C( m-l+1:m, 1:n )' * V( 1:k, 1:l )'

*

         IF( L.GT.0 )

     $      CALL DGEMM( 'Transpose', 'Transpose', N, K, L, ONE,

     $                  C( M-L+1, 1 ), LDC, V, LDV, ONE, WORK, LDWORK )

*

*        W( 1:n, 1:k ) = W( 1:n, 1:k ) * T'  or  W( 1:m, 1:k ) * T

*

         CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T,

     $               LDT, WORK, LDWORK )

*

*        C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )'

*

         DO 30 J = 1, N

            DO 20 I = 1, K

               C( I, J ) = C( I, J ) - WORK( J, I )

   20       CONTINUE

   30    CONTINUE

*

*        C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...

*                            V( 1:k, 1:l )' * W( 1:n, 1:k )'

*

         IF( L.GT.0 )

     $      CALL DGEMM( 'Transpose', 'Transpose', L, N, K, -ONE, V, LDV,

     $                  WORK, LDWORK, ONE, C( M-L+1, 1 ), LDC )

*

      ELSE IF( LSAME( SIDE, 'R' ) ) THEN

*

*        Form  C * H  or  C * H'

*

*        W( 1:m, 1:k ) = C( 1:m, 1:k )

*

         DO 40 J = 1, K

            CALL DCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 )

   40    CONTINUE

*

*        W( 1:m, 1:k ) = W( 1:m, 1:k ) + ...

*                        C( 1:m, n-l+1:n ) * V( 1:k, 1:l )'

*

         IF( L.GT.0 )

     $      CALL DGEMM( 'No transpose', 'Transpose', M, K, L, ONE,

     $                  C( 1, N-L+1 ), LDC, V, LDV, ONE, WORK, LDWORK )

*

*        W( 1:m, 1:k ) = W( 1:m, 1:k ) * T  or  W( 1:m, 1:k ) * T'

*

         CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, ONE, T,

     $               LDT, WORK, LDWORK )

*

*        C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k )

*

         DO 60 J = 1, K

            DO 50 I = 1, M

               C( I, J ) = C( I, J ) - WORK( I, J )

   50       CONTINUE

   60    CONTINUE

*

*        C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...

*                            W( 1:m, 1:k ) * V( 1:k, 1:l )

*

         IF( L.GT.0 )

     $      CALL DGEMM( 'No transpose', 'No transpose', M, L, K, -ONE,

     $                  WORK, LDWORK, V, LDV, ONE, C( 1, N-L+1 ), LDC )

*

      END IF

*

      RETURN

*

*     End of DLARZB

*

      END