c-----------------------------------------------------------------------
c example of using LINPACK to solve the wave equation using implicit 
c methods.
c-----------------------------------------------------------------------
c Driver for implicit method
		program driver

		integer	nmax, nx
		parameter (nmax = 200)

		integer	ipvt(nmax)

		real*8	unp(nmax), unm(nmax), un(nmax)
		real*8	mat(nmax, nmax)

		real*8	rhs(nmax)

c		
		real*8	dt, dx, xmax, xmin, cfl, vx, width, ampl

		nx = 100

		xmin = -1.0
		xmax = 1.0
		vx   = 0.0

		width = 0.05
		ampl = 0.05

		cfl = 1.0

		dx = (xmax - xmin) / (nx - 1)
		dt = cfl * dx

c		Initialize data
		call initialp(nx, dt, dx, width, ampl,
     .            vx,xmin, xmax, un)
		call copy(nx, un, unm)

c		Carry out one implicit step

		call impl_update(nx, dt, dx, unp, unm, un, mat,
     .           ipvt, rhs)

		call  dumpdata1d(nx, unp, dx,xmin, 'output.dat')

		stop
		end

c-----------------------------------------------------------------------
c The following subroutine takes in as its arguments u at n and n-1 and
c computes u at n+1. The routine requires a temporary array of size NxN
c where N is the size of the 1D grid on which we are solving the wave
c equation.
c To form the matrix we have to take into account the boundary conditions.
c Here we take zero boundary conditions at u_1=0 and u_N=0

c nx := grid spacing
c dt := time step
c dx := space step
c unp:= u at n+1
c unm:= u at n-1
c un := u at n
c mat:= Matrix to be formed
c 
c ipvt := Integer array of dimension N. Pivot array for LINPACK
c rhs  := array that will contain rhs of matrix equation and solution on return
c         from LINPACK. Could recycle unp and use as rhs and solution...

      subroutine impl_update(nx, dt, dx, unp, unm, un, mat,
     .           ipvt, rhs)
         implicit none

         integer   nx
      
         real*8   dt, dx

c         Grid function declarations
         real*8   unp(nx), unm(nx), un(nx)

c         Matrix and LINPACK arrays
         real*8   mat(nx,nx), ipvt(nx), rhs(nx)

c      Local variables

         integer   i, j, job

         real*8   m, m2, fact, rcond

         logical trace

         trace = .TRUE.

         m = dt / dx
         m2= m**2
         fact = 1.0 + 2.0 * m2

			write(0,*)m,m2,fact

c      First zero out matrix

         do i = 1, nx
           do j = 1, nx
             mat(i,j) = 0.0
           end do
         end do

c      Form the matrix taking into account the boundary conditions
c      Here we take fixed boundaries. The fixed boundaries yield that
c     unp(1) = 0.0 and unp(nx) = 0.0
c     We are left with a nx-2 by nx-2 matrix to solve for the interior points

c        First initialize row 1 and row nx
         mat(1,1) = 1.0

         mat(2,2) = fact
			mat(2,3) = -m2

         mat(nx-1,nx-1) = fact
			mat(nx-1,nx-2) = -m2

         mat(nx,nx) = 1.0

c         Set up rest of matrix

         do i = 3, nx-2
           mat(i,i-1) = -m2
           mat(i,i) = fact   
           mat(i,i+1) = -m2
         end do

			write(0,*)mat(1,1), mat(1,2)
			do i = 2, 5
				write(0,*)mat(i,i-1), mat(i,i), mat(i,i+1)
			end do
	

c      Form right hand side of matrix equation
c      Don't forget that we require the interior only

			rhs(1) = 0.0
         do i = 2, nx-1
           rhs(i) = 2.0 * un(i) - unm(i)
         end do
			rhs(nx) = 0.0

c      Now use DGECO to factorize matrix and get condition number
c     Could use DGEFA to get just factorization...
c      unp passed in as workspace
         call DGECO(mat,nx, nx, ipvt, rcond, unp)

         if(trace)then
           write(0,*)'Reciprocal condition number = ',rcond
         end if
         if(1.0 + rcond .eq. 1.0)then
            write(0,*)'Matrix formed is singular!!! Exiting...'
            stop
         end if

c      Use DGESL to solve the linear system
         job = 0
         call DGESL(mat, nx, nx, ipvt, rhs, job)

c      Copy the solution from rhs to unp

         unp(1) = 0.0
         do i = 2, nx-1
            unp(i) = rhs(i)
         end do
         unp(nx) = 0.0

         return
      end
c-----------------------------------------------------------------------
c  routine to copy arrays
      subroutine copy(nx, fromme, toyou)
         implicit none

         integer  nx
  
         real*8   fromme(nx), toyou(nx)

c     Local variables
         integer  i

         do i = 1, nx

            toyou(i) = fromme(i)

         end do
         return
      end
c-----------------------------------------------------------------------
c	Routine to set up initial data
c  put down Gaussian at center of grid
c	Pass h into program

		subroutine initialp(nx,  dt, h, width, ampl, 
     .				vx,xmin, xmax, phi0)
			implicit none

			integer	nx

			real*8	dt, h, dz, width, ampl, vx, 
     .				xmin, xmax
	
			real*8	phi0(nx)

c		Local variables

			integer	i, j, k

			real*8	x, t


c		Carry out simple minded initial data for dtphi0 for now and
c		Gaussian wave-packet for phi0

			t = 0.0

			do i = 1, nx

				 x = h * (i-1) + xmin
				 if( x .ge. -0.2 .and. x .le. 0.2)then
				   phi0(i) = ampl * exp(-((x - vx*t)**2)/width**2) 
     .				+	ampl * exp(-((x +  vx*t)**2)/width**2)
				 else
					phi0(i) = 0.0
				 end if

			end do

			return
		end	

c-----------------------------------------------------------------------
      subroutine dumpdata1d(nx, array, dx,xmin, fname)
         implicit none

         integer  nx

         real*8   array(nx), dx, xmin

         character*64 fname

c     Local

         integer  i,j,k

         open(unit=1,file=fname)

         do i = 1, nx
               write(1,'(2F20.15)')(i-1)*dx + xmin,array(i)
         end do
         close(unit=1)
         return
      end
c-----------------------------------------------------------------------