program CN1

c   Kreider 5/23/03  CN on 2d square
c   Gauss Elimination version

      implicit none
      double precision A(3000,3000),u(3000),b(3000),x(100),y(100)
      double precision L(3000),S(3000)
      double precision dx,dy,dt,alpha,rx,ry,d
      integer i,j,k,N,M,row

c     this is just ONE time step
c     assume initial condition is u(x,y,0) = 0 so righthand side is
c     source term only

c     assume homogeneous Dirichlet boundary conditions as well

      N = 61
      M = 41
      dx = 1.d0/dfloat(N-1)
      dy = 1.d0/dfloat(M-1)

      do i=1,N
         x(i) = dfloat(i-1)*dx
      end do
      do j=1,M
         y(j) = dfloat(j-1)*dy
      end do

      dt = 1.d-3
      alpha = 1.d0
      rx = alpha*dt/dx**2
      ry = alpha*dt/dy**2
      d = 1.d0 + rx + ry

      do i=1,N
         do j=1,M
            A(i,j) = 0.d0
         end do
      end do
      do i=1,N*M
         b(i) = 0.d0
      end do

c     build matrix A
      j=1
      do i=1,N
         row = (j-1)*N+i
         A(row,row) = 1.d0
      end do
      do j=2,M-1   
         i=1
         row = (j-1)*N+i
         A(row,row) = 1.d0
         do i=2,N-1
            row = (j-1)*N+i
            A(row,row) = d
            A(row,row-1) = -rx/2.d0
            A(row,row+1) = -rx/2.d0
            A(row,row-N) = -ry/2.d0
            A(row,row+N) = -ry/2.d0
         end do
         i=N
         row = (j-1)*N+i
         A(row,row) = 1.d0
      end do
      j=M
      do i=1,N
         row = (j-1)*N+i
         A(row,row) = 1.d0
      end do

c     build righthand side b
      do j=2,M-1
         do i=2,N-1
            row = (j-1)*N+i
            b(row) = 2.d0*exp(x(i)*y(j))
         end do
      end do

      call gauss(N*M,A,3000,L,S)
      call solve(N*M,A,3000,L,b,u)

      open(1,file='CN1.dat')
      do i=1,N
         do j=1,M
            write(1,*) x(i),y(j),u((j-1)*N+i)
         end do
         write(1,*)
      end do

      stop
      end

      SUBROUTINE GAUSS(N,A,IA,L,S)
      double precision A(IA,N),L(N),S(N),smax,rmax,r,LK,xmult
      DO 3 I = 1,N
        L(I) = I
        SMAX = 0.0
        DO 2 J = 1,N
          SMAX = DMAX1(SMAX,DABS(A(I,J)))
   2    CONTINUE
        S(I) = SMAX
   3  CONTINUE
      DO 7 K = 1,N-1
        RMAX = 0.0
        DO 4 I = K,N
          R = ABS(A(L(I),K))/S(L(I))
          IF(R .LE. RMAX)  GO TO 4
          J = I
          RMAX = R
   4    CONTINUE
        LK = L(J)
        L(J) = L(K)
        L(K) = LK
        DO 6 I = K+1,N
          XMULT = A(L(I),K)/A(LK,K)
          DO 5 J = K+1,N
            A(L(I),J) = A(L(I),J) - XMULT*A(LK,J)
   5      CONTINUE
          A(L(I),K) = XMULT
   6    CONTINUE
   7  CONTINUE
      RETURN
      END

      SUBROUTINE SOLVE(N,A,IA,L,B,X)
      double precision A(IA,N),L(N),B(N),X(N)
      double precision sum
      DO 3 K = 1,N-1
        DO 2 I = K+1,N
          B(L(I)) = B(L(I)) - A(L(I),K)*B(L(K))
   2    CONTINUE
   3  CONTINUE
      X(N) = B(L(N))/A(L(N),N)
      DO 5 I = N-1,1,-1
        SUM = B(L(I))
        DO 4 J = I+1,N
          SUM = SUM - A(L(I),J)*X(J)
   4    CONTINUE
        X(I) = SUM/A(L(I),I)
   5  CONTINUE
      RETURN
      END