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 )
*
*  ===================================================================