! pdenu_eq.f90  2D non uniform grid boundary value PDE  non iterative 
!                known solution to test method                         
!                u(x,y) = exp(x*y)*sin(x+y)
!                                                                      
!                differential equation to solve                       
! d^2u/dx^2 + 2*d^2u/dy^2 - 3*d^2u/dxdy + y*du/dx - x*du/dy = c(x,y)
! exp(x*y)*sin(x+y)*(x*x + 2*y*y -3*x*y -3)
! non uniform grid on rectangle -1.0,-1.1 to 1.2,1.3                
!                                                                   
! set up and solve system of linear equations                       
!                                                                   
! create two dimensional array with                                 
! (nx-1)*(ny-1) rows, one for each equation [nx in fortran is last subscript]
!
! uses nuderiv.f90
! uses simeq.f90


module global  
  implicit none
  integer :: nx = 10 ! one less than 'C' when Fortran subscript starts at zero
  integer :: ny = 9  ! includes boundary
  double precision, dimension(0:71) :: us ! solution being iterated        
  double precision, dimension(0:71,0:72) :: ut ! temporary for partial results 
  double precision, dimension(0:10) :: xg
  data xg /-1.0, -0.9, -0.8, -0.6, -0.3, 0.0, &
            0.35, 0.65, 0.85, 1.0,  1.2/
  double precision, dimension(0:9) :: yg
  data yg /-1.1, -1.0, -0.9, -0.6, -0.2, &
            0.3,  0.8,  1.0,  1.2,  1.3/

  double precision, dimension(0:10) :: cxx    ! nuderiv coefficients 
  double precision, dimension(0:9)  :: cyy
  double precision, dimension(0:10,0:9) :: cxy
  double precision, dimension(0:10) :: cx
  double precision, dimension(0:9)  :: cy
  double precision :: coefdxx, coefdyy, coefdxdy, coefdx, coefdy, coefu;
  double precision :: error         ! sum of absolute exact-computed         
  double precision :: avgerror      ! average error                          
  double precision :: maxerror      ! maximum cell error                     
  interface
    function u(x, y) result(value) ! solution for boundary and check 
      implicit none
      double precision, intent(in) :: x, y
      double precision :: value
    end function u
  end interface
  interface
    function c(x, y) result(term) ! from solve
      implicit none
      double precision, intent(in) :: x, y
      double precision :: term
    end function c
  end interface
  interface
    subroutine simeq(n, B, X)
      implicit none
      integer, intent(in) :: n
      double precision, dimension(0:71,0:72), intent(in) :: b
      double precision, dimension(0:71), intent(out) :: X
     end subroutine simeq
  end interface
  interface
     subroutine nuderiv(order, npoint, point, X, C)
     implicit none
     integer, intent(in) :: order
     integer, intent(in) :: npoint
     integer, intent(in) :: point
     double precision, dimension(npoint), intent(in) :: X
     double precision, dimension(npoint), intent(inout) :: C
     end subroutine nuderiv
  end interface
end module global

program pdenu_eq
  use global
  implicit none
  integer :: i, j
  
  print *, "pdenu_eq.f90 running"
  print *, "differential equation to solve"                       
  print *, "d^2u/dx^2 + 2*d^2u/dy^2 + 3*d^2u/dxdy + 4*du/dx + 5*du/dy +6*u = c(x,y)"
  print *, "c(x,y) = 6x^3+12y^3+18x^2y+24xy^2+57x^2+82y^2+64xy+122x+114y+110" 
  print *, "non uniform grid on rectangle -1.0,-1.1 to 1.2,1.3" 
  print *, "known solution to test method"
  print *, "u(x,y) = x^3+2y^3+3x^2y+4xy^2+5x^2+6y^2+7x+8"
  print *, "xmin=",xg(0), "xmax=",xg(nx), "nx=",nx
  print *, "ymin=",yg(0), "ymax=",yg(ny), "ny=",ny
  print *, "u(xg(0),yg(0))=",u(xg(0),yg(0))
  print *, "c(xg(0),yg(0))=",c(xg(0),yg(0))
  
  call init_matrix()
  print *, "matrix initialized"

  print "(/a)", "initial matrix, left side, upper"
  do i=0,7
    print "(8f9.3)", (ut(i,j),j=0,7)
  end do
  print "(/a)", "initial matrix, right side, upper"
  do i=0,7 
    print "(8f9.3)", (ut(i,j),j=65,72)
  end do
  print "(/a)", "initial matrix, right side, lower"
  do i=64,71
    print "(8f9.3)", (ut(i,j),j=65,72)
  end do
  print *, " "

  call simeq(71, ut, us)
  call printexact()

  call check()
  avgerror = error/((nx-1)*(ny-1))
  print "('error=',e15.8,' avg=',e15.8,' max=',e15.8)", &
           error, avgerror, maxerror
end program pdenu_eq 

function u(x, y) result(value) ! solution for boundary and check 
  implicit none
  double precision, intent(in) :: x, y
  double precision :: value
  value = exp(x*y)*sin(x+y)
end function u

function c(x, y) result(term) ! from solve
  implicit none
  double precision, intent(in) :: x, y
  double precision :: term
  term = exp(x*y)*sin(x+y)*(x*x + 2.0*y*y -3.0*x*y -3.0)
end function c

subroutine check() ! typically not known, instrumentation 
  use global
  implicit none
  integer :: i, j
  double precision :: uexact, err
  
  error = 0.0
  maxerror = 0.0
  do i=1,nx-1
    do j=1,ny-1
      uexact = u(xg(i),yg(j))
      err = abs(us((i-1)+(nx-1)*(j-1))-uexact)
      error = error + err
      if(err>maxerror) maxerror=err
    end do
  end do
end subroutine check 

subroutine printexact()
  use global
  implicit none
  integer :: i, j
  
  print "(/a)", "exact solution u, computed us, error"
  do i=1,nx-1
    do j=1,ny-1
      print "('xg=',f6.3,', yg=',f6.3,', u=',f6.3,', us=',f9.3,'err=',e15.7)", &
      xg(i), yg(j), u(xg(i),yg(j)), us((i-1)+(nx-1)*(j-1)), &
      u(xg(i),yg(j))-us((i-1)+(nx-1)*(j-1))
    end do
  end do
end subroutine printexact 

subroutine init_matrix()
  use global
  implicit none
  integer :: i, j, k, ii, jj, kk, cs 
  integer :: ik, kj, kx, ky, kxky
  
  cs = (nx-1)*(ny-1) ! constant column subscript 
  
  ! zero internal cells 
  do i=1,nx-1
    do j=1,ny-1
      do ii=1,nx-1
        do jj=1,ny-1
          ut((i-1)+(nx-1)*(j-1), (ii-1)+(nx-1)*(jj-1)) = 0.0
        end do ! jj
      end do ! ii
    end do ! j
  end do ! i

  print *, "internal cells zeroed "          
  do i=1,nx-1  ! inside boundary 
    do j=1,ny-1
      ii = (i-1)+(nx-1)*(j-1) ! column of matrix 
      ! for each term 

      ! coefdxx * d^2 u(x,y)/dx^2 
      coefdxx = 1.0
      call nuderiv(2,nx+1,i+1,xg,cxx)
      do k=0,nx
        cxx(k) = cxx(k) * coefdxx
        kj = (k-1)+(nx-1)*(j-1)
	if(k.eq.0 .or. k.eq.nx) then
          ut(ii,cs) = ut(ii,cs) - cxx(k)*u(xg(k),yg(j))
	else
          ut(ii,kj) = ut(ii,kj) + cxx(k)
          !  ut(ii, kj=(i-1)+1)+(nx-1)*(j-1) ) = 1.0/hx2
        end if
      end do

      ! coefdyy * d^2 u(x,y)/dy^2 
      coefdyy = 2.0
      call nuderiv(2,ny+1,j+1,yg,cyy)
      do k=0,ny
        cyy(k) = cyy(k) * coefdyy
        ik = (i-1)+(nx-1)*(k-1)
	if(k.eq.0 .or. k.eq.ny) then
          ut(ii,cs) = ut(ii,cs) - cyy(k)*u(xg(i),yg(k))
	else
          ut(ii,ik) = ut(ii,ik) + cyy(k)
        end if
      end do

      ! coefdxdy * d^2 u(x,y)/dxdy 
      coefdxdy = -3.0
      call nuderiv(1,nx+1,i+1,xg,cx)
      call nuderiv(1,ny+1,j+1,yg,cy)
      do kx=0,nx
        do ky=0,ny
          cxy(kx,ky) = cx(kx) * cy(ky) * coefdxdy
          kxky = (kx-1)+(nx-1)*(ky-1)
	  if(kx.eq.0 .or. kx.eq.nx .or. ky.eq.0 .or. ky.eq.ny) then
            ut(ii,cs)   = ut(ii,cs)   - cxy(kx,ky)*u(xg(kx),yg(ky))
	  else
            ut(ii,kxky) = ut(ii,kxky) + cxy(kx,ky)
          end if
        end do
      end do

      ! coefdx * d u(x,y)/dx 
      coefdx = yg(j)
      call nuderiv(1,nx+1,i+1,xg,cx)
      do k=0,nx
        cx(k) = cx(k) * coefdx
        kj = (k-1)+(nx-1)*(j-1)
        if(k.eq.0 .or. k.eq.nx) then
          ut(ii,cs) = ut(ii,cs) - cx(k)*u(xg(k),yg(j))
	else
          ut(ii,kj) = ut(ii,kj) + cx(k)
        end if
      end do

      ! coefdy * d u(x,y)/dy 
      coefdy = -xg(i)
      call nuderiv(1,ny+1,j+1,yg,cy)
      do k=0,ny
        cy(k) = cy(k) * coefdy
        ik = (i-1)+(nx-1)*(k-1)
	if(k.eq.0 .or. k.eq.ny) then
          ut(ii,cs) = ut(ii,cs) - cy(k)*u(xg(i),yg(k))
	else
          ut(ii,ik) = ut(ii,ik) + cy(k)
        end if
      end do

      ! coefu * u(x,y) 
      coefu = 0.0
      ut(ii,ii) = ut(ii,ii) + 1.0*coefu

      ! RHS constant 
      ut(ii,cs) = ut(ii,cs) + c(xg(i),yg(j))
      ! end terms 
    end do
  end do
end subroutine init_matrix
  
subroutine simeq(n, B, X)
  implicit none ! tailored version for this application
  integer, intent(in) :: n
  double precision, dimension(0:71,0:72), intent(inout) :: B
  double precision, dimension(0:71), intent(out) :: X
  integer, dimension(0:n) :: ROW           ! ROW INTERCHANGE INDICES 
  integer :: HOLD , I_PIVOT                ! PIVOT INDICES 
  double precision :: PIVOT                ! PIVOT ELEMENT VALUE 
  double precision :: ABS_PIVOT
  integer :: i, j, k, m

  m = n+1

  ! SET UP ROW  INTERCHANGE VECTORS 
  do k=0,n
    ROW(k) = k
  end do ! k

  ! BEGIN MAIN REDUCTION LOOP 
  do k=0,n

    ! FIND LARGEST ELEMENT FOR PIVOT 
    PIVOT = B(ROW(k),k)
    ABS_PIVOT = abs(PIVOT)
    I_PIVOT = k
    do i=k,n
      if( abs(B(ROW(i),k)) > ABS_PIVOT) then
        I_PIVOT = i
        PIVOT = B(ROW(i),k)
        ABS_PIVOT = abs ( PIVOT )
      end if
    end do ! i

    ! HAVE PIVOT, INTERCHANGE ROW POINTERS 
    HOLD = ROW(k)
    ROW(k) = ROW(I_PIVOT)
    ROW(I_PIVOT) = HOLD

    ! CHECK FOR NEAR SINGULAR 
    if( ABS_PIVOT < 1.0E-10 ) then
      do j=k+1,n+1
        B(ROW(k),j) = 0.0
      end do ! j
      print *, 'redundant row (singular) ', ROW(k)
    else

      ! REDUCE ABOUT PIVOT 
      do j=k+1,n+1
        B(ROW(k),j) = B(ROW(k),j) / B(ROW(k),k)
      end do ! j

      ! INNER REDUCTION LOOP 
      do i=0,n
        if( i .ne. k) then
          do j=k+1,n+1
            B(ROW(i),j) = B(ROW(i),j) - B(ROW(i),k) * B(ROW(k),j)
          end do ! j
        end if
      end do ! i
    end if 
    ! FINISHED INNER REDUCTION 
  end do ! k

  ! END OF MAIN REDUCTION LOOP 
  ! BUILD  X  FOR RETURN, UNSCRAMBLING ROWS 
  do i=0,n
    X(i) = B(ROW(i),n+1)
  end do ! i
end subroutine simeq