subroutine read ( n, matrix, usol, rhsd, nx)
c ==================================================================
c
c     programmer       Kees Vuik
c
c     version   : 1.0
c     purpose   : reading of the matrix
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 read the matrix for the
c     Laplace symphony.
c
c ******************************************************************
c
c                       KEYWORDS
c
c
c ******************************************************************
c
c                       INPUT / OUTPUT   PARAMETERS
c
c matrix o    contents of the coefficient matrix
c n      i    product of (nx+1) and (ny+1)
c rhsd   o    contents of the right-hand-side vector
c usol   o    contents of the solution vector
c
c
c ******************************************************************
c
c                       COMMON BLOCKS
c
c ******************************************************************
c
c                       LOCAL PARAMETERS
c
c     i        local counting variable
c     j        local counting variable
c
c ******************************************************************
c
c                       SUBROUTINES CALLED
c
c ******************************************************************
c
c                       I/O
c
c ******************************************************************
c
c                       ERROR MESSAGES
c
c
c ******************************************************************
c
c                       PSEUDO CODE
c
c ==================================================================
c
c     Parameters:
c
      integer          n, nx
      double precision matrix(1:n,1:9), usol(n), rhsd(n)

      if ( nx .eq. 16 ) open(16,file = 'matrix16.dat')
      if ( nx .eq. 32 ) open(16,file = 'matrix32.dat')
      if ( nx .eq. 64 ) open(16,file = 'matrix64.dat')
      do j = 1, 9
         do i = 1, n
            read(16,*) matrix(i,j)
         enddo
      enddo
      do i = 1, n
         read(16,*) usol(i)
      enddo
      do i = 1, n
         read(16,*) rhsd(i)
      enddo
      close(16)
      end