SUBROUTINE MVMODE( N, M, X, LDX, Y, LDY )
*     .. 
*     .. Scalar Arguments ..
      INTEGER            LDY, LDX, M, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   Y( LDY, * ), X( LDX, * )
*     ..
*
*  Purpose
*  =======
*
*  Compute
*
*               Y(:,1:M) = A*X(:,1:M)
*
*  where op(A) is A or A' (the transpose of A). The matrix A is from
*  the discretization of an ODE boundary value problem.
*
*  NOTE: N should be smaller than 1001. If one wants to test a larger
*  size of N, the following parameter LDU should be changed to the 
*  desirable size. 
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix A.
*
*  M       (input) INTEGERS
*          The number of columns of X to multiply.
*
*  X       (input) REAL array, dimension ( LDX, M )
*          X contains the matrix X.
* 
*  LDX     (input) INTEGER
*          The leading dimension of the array X, LDX >= max( 1,N )
*
*  Y       (output) REAL array, dimension (LDY, M )
*          contains the product of the matrix A with Y.
*
*  LDY     (input) INTEGER
*          The leading dimension of the array Y, LDY >= max( 1,N )
*
*  Workspace: there are need 3*N workspace for computing the matrix-
*  vector product. 
*
*  ===================================================================