SUBROUTINE MVMRWK( 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)
*
*  The matrix is from a random walk example. The order of the random 
*  walk test matrix is N = (K+1)*(K+2)/2,  where K is the size of 
*  triangular random walk grid.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix A.
*          N should be equal to (K+1)*(K+2)/2 for an integer K > 0.
*
*  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 X.
*
*  LDY     (input) INTEGER
*          The leading dimension of the array y, LDY >= max(1,N)
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, HALF
      PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J, K, KK, LK, NN
      DOUBLE PRECISION   UPPER
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          REAL, SQRT
*     ..
*     .. Statement Functions ..
      INTEGER            IINDEX
*     ..
*     .. Statement Function definitions ..
      IINDEX( K, I, J ) = I*K - ( ( I-1 )*( I-2 ) ) / 2 + 1 + J + 1
*     ..
*     .. Executable Statements ..
*
      NN = ( -3+SQRT( REAL( 8*N+1 ) ) ) / 2
*
*
*        Compute Y( :,1:M ) = A*X( :,1:M ) 
*
         DO 30 KK = 1, M
*
*           Compute the vector Y(1:N,KK) = A*X(1:N,KK)
*
            DO 20 I = 0, NN
               DO 10 J = 0, NN - I
*
                  K = IINDEX( NN, I, J )
                  UPPER = REAL( I+J ) / REAL( 2*NN )
                  Y( K, KK ) = ZERO
*
                  IF( I.EQ.0 .AND. J.EQ.0 ) THEN
                     LK = IINDEX( NN, I, J+1 )
                     Y( K, KK ) = Y( K, KK ) + 
     $                             ( HALF-UPPER )*X( LK, KK )
                     LK = IINDEX( NN, I+1, J )
                     Y( K, KK ) = Y( K, KK ) + 
     $                             ( HALF-UPPER )*X( LK, KK )
*
                  ELSE IF( I.EQ.0 .AND. J.LE.NN-1 ) THEN
                     LK = IINDEX( NN, I, J-1 )
                     Y( K, KK ) = Y( K, KK ) + 2*UPPER*X( LK, KK )
                     LK = IINDEX( NN, I+1, J )
                     Y( K, KK ) = Y( K, KK ) + 
     $                             ( HALF-UPPER )*X( LK, KK )
                     LK = IINDEX( NN, I, J+1 )
                     Y( K, KK ) = Y( K, KK ) + 
     $                             ( HALF-UPPER )*X( LK, KK )
*
                  ELSE IF( I.EQ.0 .AND. J.EQ.NN ) THEN
                     LK = IINDEX( NN, I, J-1 )
                     Y( K, KK ) = Y( K, KK ) + 2*UPPER*X( LK, KK )
*
                  ELSE IF( I.LE.NN-1 .AND. J.EQ.0 ) THEN
                     LK = IINDEX( NN, I-1, J )
                     Y( K, KK ) = Y( K, KK ) + 2*UPPER*X( LK, KK )
                     LK = IINDEX( NN, I+1, J )
                     Y( K, KK ) = Y( K, KK ) + 
     $                             ( HALF-UPPER )*X( LK, KK )
                     LK = IINDEX( NN, I, J+1 )
                     Y( K, KK ) = Y( K, KK ) + 
     $                             ( HALF-UPPER )*X( LK, KK )
*
                  ELSE IF( I.EQ.NN .AND. J.EQ.0 ) THEN
                     LK = IINDEX( NN, I-1, J )
                     Y( K, KK ) = Y( K, KK ) + 2*UPPER*X( LK, KK )
*
                  ELSE IF( ( I+J ).EQ.NN ) THEN
                     LK = IINDEX( NN, I-1, J )
                     Y( K, KK ) = Y( K, KK ) + UPPER*X( LK, KK )
                     LK = IINDEX( NN, I, J-1 )
                     Y( K, KK ) = Y( K, KK ) + UPPER*X( LK, KK )
*
                  ELSE
                     LK = IINDEX( NN, I-1, J )
                     Y( K, KK ) = Y( K, KK ) + UPPER*X( LK, KK )
                     LK = IINDEX( NN, I, J-1 )
                     Y( K, KK ) = Y( K, KK ) + UPPER*X( LK, KK )
                     LK = IINDEX( NN, I, J+1 )
                     Y( K, KK ) = Y( K, KK ) + 
     $                             ( HALF-UPPER )*X( LK, KK )
                     LK = IINDEX( NN, I+1, J )
                     Y( K, KK ) = Y( K, KK ) + 
     $                             ( HALF-UPPER )*X( LK, KK )
*
                  END IF
   10          CONTINUE
   20       CONTINUE
*
   30    CONTINUE
*
      RETURN
*
*     End of MVMRWK
*
      END