subroutine matvec ( matrix, x, b, n, nx )

c ==================================================================
c
c     programmer       Kees Vuik
c
c     version   : 1.0
c     purpose   : matrix vector product
c     date      : 07-06-2004
c     programmer: Kees Vuik
c     contact   : c.vuik@math.tudelft.nl
c                 TU Delft, The Netherlands
c
c
c ******************************************************************
c
c                       DESCRIPTION
c
c     The subroutine read is used to compute a matrix vector product
c     for the Laplace symphony.
c
c ******************************************************************

      integer          n, nx
      double precision matrix(1:n,1:9), x(n), b(n)
      integer          i,n1

      n1 = nx+1
      i = 1
         b(i) = matrix(i,1)*x(i)+
     +          matrix(i,6)*x(i+1)+
     +          matrix(i,7)*x(i+n1-1)+
     +          matrix(i,8)*x(i+n1)+
     +          matrix(i,9)*x(i+n1+1)
      do i = 2, n1-1
         b(i) = matrix(i,1)*x(i)+
     +          matrix(i,5)*x(i-1)+
     +          matrix(i,6)*x(i+1)+
     +          matrix(i,7)*x(i+n1-1)+
     +          matrix(i,8)*x(i+n1)+
     +          matrix(i,9)*x(i+n1+1)
      end do
      i = n1
         b(i) = matrix(i,1)*x(i)+
     +          matrix(i,4)*x(i-n1+1)+
     +          matrix(i,5)*x(i-1)+
     +          matrix(i,6)*x(i+1)+
     +          matrix(i,7)*x(i+n1-1)+
     +          matrix(i,8)*x(i+n1)+
     +          matrix(i,9)*x(i+n1+1)
      i = n1+1
         b(i) = matrix(i,1)*x(i)+
     +          matrix(i,3)*x(i-n1)+
     +          matrix(i,4)*x(i-n1+1)+
     +          matrix(i,5)*x(i-1)+
     +          matrix(i,6)*x(i+1)+
     +          matrix(i,7)*x(i+n1-1)+
     +          matrix(i,8)*x(i+n1)+
     +          matrix(i,9)*x(i+n1+1)
      do i = n1+2, n-n1-1
	 b(i) = matrix(i,1)*x(i)+
     +          matrix(i,2)*x(i-n1-1)+
     +          matrix(i,3)*x(i-n1)+
     +          matrix(i,4)*x(i-n1+1)+
     +          matrix(i,5)*x(i-1)+
     +          matrix(i,6)*x(i+1)+
     +          matrix(i,7)*x(i+n1-1)+
     +          matrix(i,8)*x(i+n1)+
     +          matrix(i,9)*x(i+n1+1)
      end do
      i = n-n1
         b(i) = matrix(i,1)*x(i)+
     +          matrix(i,2)*x(i-n1-1)+
     +          matrix(i,3)*x(i-n1)+
     +          matrix(i,4)*x(i-n1+1)+
     +          matrix(i,5)*x(i-1)+
     +          matrix(i,6)*x(i+1)+
     +          matrix(i,7)*x(i+n1-1)+
     +          matrix(i,8)*x(i+n1)
      i = n-n1+1
         b(i) = matrix(i,1)*x(i)+
     +          matrix(i,2)*x(i-n1-1)+
     +          matrix(i,3)*x(i-n1)+
     +          matrix(i,4)*x(i-n1+1)+
     +          matrix(i,5)*x(i-1)+
     +          matrix(i,6)*x(i+1)+
     +          matrix(i,7)*x(i+n1-1)
      do i = n-n1+2, n-1
         b(i) = matrix(i,1)*x(i)+
     +          matrix(i,2)*x(i-n1-1)+
     +          matrix(i,3)*x(i-n1)+
     +          matrix(i,4)*x(i-n1+1)+
     +          matrix(i,5)*x(i-1)+
     +          matrix(i,6)*x(i+1)
      end do
      i = n
         b(i) = matrix(i,1)*x(i)+
     +          matrix(i,2)*x(i-n1-1)+
     +          matrix(i,3)*x(i-n1)+
     +          matrix(i,4)*x(i-n1+1)+
     +          matrix(i,5)*x(i-1)
      return
      end