      SUBROUTINE MVMBWM( TRANS, N, M, X, LDX, Y, LDY )
*     ..
*     .. Scalar Arguments ..
      INTEGER            LDY, LDX, M, N, TRANS
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   Y( LDY, * ), X( LDX, * )
*     ..
*
*  Purpose
*  =======
*
*  Compute
*
*               Y(:,1:M) = op(A)*X(:,1:M)
*
*  where op(A) is A or A' (the transpose of A). The matrix A is a 2 by 2
*  block matrix resulted from the finite difference discretization
*  of the Brusselator wave model. 
*
*  Arguments
*  =========
*
*  TRANS   (input) INTEGER
*          If TRANS = 0, compute Y(:,1:M) = A*Q(:,1:M) 
*          If TRANS = 1, compute Y(:,1:M) = A'*Q(:,1:M) 
*           
*  N       (input) INTEGER
*          The order of the matrix A. N has to be an even number.
*
*  M       (input) INTEGERS
*          The number of Q to multiply.
*
*  X       (input) DOUBLE PRECISION array, dimension ( LDX, M )
*          contains the vectors X.
*
*  LDX     (input) INTEGER
*          The leading dimension of the array X, LDX >= max( 1, N )
*
*  Y       (output) DOUBLE PRECISION array, dimension (LDY, M )
*          contains the product of the matrix op(A) with X.
*
*  LDY     (input) INTEGER
*          The leading dimension of the array Y, LDY >= max( 1, N )
*
*  ===================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, TWO
      PARAMETER          ( ONE = 1.0D+0, TWO = 2.0D+0 )
      DOUBLE PRECISION   ALPHA, BETA, DELT1, DELT2, L
      PARAMETER          ( ALPHA = 2.0D+0, BETA = 5.45D+0, 
     $                     DELT1 = 8.0D-3, DELT2 = 4.0D-3,
     $                     L = 0.51302D+0 )
*
*     .. Local Variables .. 
      INTEGER            I, K, NHALF
      DOUBLE PRECISION   HSTEP, TAU1, TAU2, DIAG1, DIAG2, ALPHA2 
*
*     set up middle parameters
*
      NHALF = N / 2 
      HSTEP = ONE / ( DBLE( NHALF ) + 1 ) 
      TAU1 = DELT1 / (HSTEP*L)**2  
      TAU2 = DELT2 / (HSTEP*L)**2  
      DIAG1 = -TWO*TAU1 + BETA - ONE 
      ALPHA2 = ALPHA**2 
      DIAG2 = -TWO*TAU2 - ALPHA2 
* 
      IF ( TRANS.EQ.0 ) THEN
*
*        Compute Y(:,1*M) = A*X(:,1:M)  
*
         DO 30 K = 1, M 
*
*           the 1st row 
*
            Y( 1, K ) = DIAG1*X( 1, K ) + TAU1*X( 2, K ) + 
     $                   ALPHA2*X( NHALF+1, K )    
*
*           the 2nd row to (NHALF-1)th row              
*
            DO 10 I = 2, NHALF - 1 
               Y( I, K ) = TAU1*X( I-1, K ) + DIAG1*X( I, K ) +
     $                      TAU1*X( I+1, K ) + ALPHA2*X( NHALF+I, K )   
   10       CONTINUE  
*
*           the NHALFth row
*
            Y( NHALF, K ) = TAU1*X( NHALF-1, K ) + DIAG1*X( NHALF, K ) 
     $                       + ALPHA2*X( N,K )  
*          
*           the (NHALF+1)th row
*
            Y( NHALF+1, K ) = -BETA*X( 1, K ) + DIAG2*X( NHALF+1, K ) 
     $                          + TAU2*X( NHALF+2, K )  
*
*           the (NHALF+2)th row to (N-1)th row
*
            DO 20 I = NHALF+2, N - 1 
               Y( I, K ) = -BETA*X( I-NHALF, K ) + TAU2*X( I-1, K ) + 
     $                       DIAG2*X( I, K ) + TAU2*X( I+1, K ) 
   20       CONTINUE     
*
*           the last (Nth) row  
*
            Y( N, K ) = -BETA*X( N-NHALF, K ) + TAU2*X( N-1,K ) + 
     $                   DIAG2*X( N, K ) 
*
   30    CONTINUE 
*
      ELSE IF ( TRANS.EQ. 1 ) THEN
* 
*        Compute Y(:,1:M) = A'*X(:,1:M)  
*
         DO 60 K = 1, M 
*
*           the 1st row 
*
            Y( 1, K ) = DIAG1*X( 1, K ) + TAU1*X( 2, K )  
     $                   - BETA*X( NHALF+1, K )    
*
*           the 2nd row to (NHALF-1)th row              
*
            DO 40 I = 2, NHALF - 1 
               Y( I, K ) = TAU1*X( I-1, K ) + DIAG1*X( I, K ) +
     $                      TAU1*X( I+1, K ) - BETA*X( NHALF+I, K )   
   40       CONTINUE  
*
*           the NHALFth row
*
            Y( NHALF, K ) = TAU1*X( NHALF-1, K ) + DIAG1*X( NHALF, K ) 
     $                       - BETA*X( N,K )  
*          
*           the (NHALF+1)th row
*
            Y( NHALF+1, K ) = ALPHA2*X( 1, K ) + DIAG2*X( NHALF+1, K ) 
     $                          + TAU2*X( NHALF+2, K )  
*
*           the (NHALF+2)th row to (N-1)th row
*
            DO 50 I = NHALF+2, N - 1 
               Y( I, K ) = ALPHA2*X( I-NHALF, K ) + TAU2*X( I-1, K ) + 
     $                       DIAG2*X( I, K ) + TAU2*X( I+1, K ) 
   50       CONTINUE     
*
*           the last (Nth) row  
*
            Y( N, K ) = ALPHA2*X( N-NHALF, K ) + TAU2*X( N-1,K ) + 
     $                   DIAG2*X( N, K ) 
*
   60    CONTINUE 
*
      END IF 
*
      RETURN
*
*     End of MVMBWM
*
      END
