SUBROUTINE MVMMCD( 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 block
*  tridiagonal matrix resulted from the finite difference discretization       
*  of a 2-D model convection diffusion operator,
*
*           L[u] = - u_xx - u_yy + 2p_1 u_x + 2p_2 u_y - p_3 u.
*
*  The constant parameters p1, p2 and p3 in the convection-diffusion 
*  operator may be changed.  
*
*  The convection-diffusion operator is discretized on a K x K square 
*  grid, and resulted in the matrix A of order N = K^2. 
*
*  Arguments
*  =========
*
*  TRANS   (input) INTEGER
*          If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M) 
*          If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M) 
*           
*  N       (input) INTEGER
*          The order of the matrix A. N has to be K*K for an integer K.
*
*  M       (input) INTEGERS
*          The number of columns of X to multiply.
*
*  X       (input) DOUBLE PRECISION array, dimension ( LDX, M )
*          X contains the matrix (vectors) X.
*
*  LDX     (input) INTEGER
*          The leading dimension of array X, LDX >= max( 1, N )
*
*  Y       (output) DOUBLE PRECISION array, dimension (LDY, M )
*          contains the product of the matrix and vectors. 
*
*  LDY     (input) INTEGER
*          The leading dimension of array Y, LDY >= max( 1, N )
*
*  ===================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, FOUR
      PARAMETER          ( ONE = 1.0D+0, FOUR = 4.0D+0 )
      DOUBLE PRECISION   P1, P2, P3
      PARAMETER          ( P1 = 1.0D+0, P2 = 1.0D+0, P3 = 1.0D+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, II, K, NS
      DOUBLE PRECISION   BETA, GAMMA, HSTEP, SIGMA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          REAL, SQRT
*     ..
*     .. Executable Statements ..
*
      NS = SQRT( REAL( N ) )
      HSTEP = ONE / ( REAL( NS )+ONE )
      SIGMA = P3*HSTEP*HSTEP
      GAMMA = P2*HSTEP
      BETA = P1*HSTEP
*
      IF ( TRANS.EQ.0 ) THEN 
         DO 70 K = 1, M
*
*           Compute the vector Y(1:N,K) = A*X(1:N,K)
*
            DO 60 I = 1, NS
*
*              diagonal elements
*
               DO 10 II = ( I-1 )*NS + 1, I*NS
                  Y( II, K ) = ( FOUR-SIGMA )*X( II, K )
   10          CONTINUE
*
*              sub diagonal elements
*
               DO 20 II = ( I-1 )*NS + 1 + 1, I*NS
                  Y( II, K ) = Y( II, K ) + 
     $                                  ( -GAMMA-ONE )*X( II-1, K )
   20          CONTINUE
*
*              super diagonal elements 
*
               DO 30 II = ( I-1 )*NS + 1, I*NS - 1
                  Y( II, K ) = Y( II, K ) + 
     $                                  ( GAMMA-ONE )*X( II+1, K ) 
   30          CONTINUE
*
*              (NS+1)th sub diagonal 
*
               IF( I.GT.1 ) THEN
                  DO 40 II = ( I-1 )*NS + 1, I*NS
                     Y( II, K ) = Y( II, K ) + 
     $                                  ( BETA-ONE )*X( II-NS, K )
   40             CONTINUE
               END IF
*
*              (NS+1)th super diagonal 
*
               IF( I.LT.NS ) THEN
                  DO 50 II = ( I-1 )*NS + 1, I*NS
                     Y( II, K ) = Y( II, K ) + 
     $                                  ( -BETA-ONE )*X( II+NS, K )
   50             CONTINUE
               END IF
*
   60       CONTINUE
*
   70    CONTINUE
*
      ELSE IF ( TRANS.EQ.1 ) THEN 
*
         DO 140 K = 1, M
*
*           Compute the vector Y(1:N,K) = trans(A)*X(1:N,K)
*
            DO 130 I = 1, NS
*
               DO 80 II = ( I-1 )*NS + 1, I*NS
                  Y( II, K ) = ( FOUR-SIGMA )*X( II, K )
   80          CONTINUE
*
               DO 90 II = ( I-1 )*NS + 1, I*NS - 1
                  Y( II, K ) = Y( II, K ) + 
     $                                  ( -GAMMA-ONE )*X( II+1, K )       
   90          CONTINUE
*
               DO 100 II = ( I-1 )*NS + 1 + 1, I*NS
                  Y( II, K ) = Y( II, K ) + 
     $                                  ( GAMMA-ONE )*X( II-1, K )
  100          CONTINUE
*
*              (NS+1)th sub diagonal
*
               IF( I.GT.1 ) THEN
                  DO 110 II = ( I-1 )*NS + 1, I*NS
                     Y( II, K ) = Y( II, K ) + 
     $                                  (-BETA-ONE )*X( II-NS, K )
  110             CONTINUE
               END IF
*
*              (NS+1)th super diagonal
*
               IF( I.LT.NS ) THEN
                  DO 120 II = ( I-1 )*NS + 1, I*NS
                     Y( II, K ) = Y( II, K ) + 
     $                                  ( BETA-ONE )*X( II+NS, K )
  120             CONTINUE
               END IF
*
  130       CONTINUE
*
  140    CONTINUE
*
      END IF 
*
      RETURN
*
*     End of MVMMCD
      END