! pdeCnu22_eq1.f90  2D non uniform grid boundary value PDE  non iterative 
!                known solution to test method  double complex
!                 U(x,y) = x^3+2y^3+3x^2y+4xy^2+5x^2+6y^2+7x+8
!
!                differential equation to solve                       
! d^2U/dx^2 + 2*d^2U/dy^2 + 3*d^2U/dxdy + 4*dU/dx + 5*dU/dy +6*U = c(x,y)
! c(x,y) = 6x^3+12y^3+18x^2y+24xy^2+57x^2+82y^2+64xy+122x+114y+110 
! non uniform grid on rectangle (-1.0,-0.9) (-1.1,-0.8) to (1.2,1.4) (1.3,1.5)
!                                                                   
! 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 nuderivC.f90  uses inverseC.f90
! uses simeqC1.f90

module global  
  implicit none
  integer, parameter :: nx = 11 ! Fortran subscript starts at 1
  integer, parameter :: ny = 10 ! includes boundary (11-2)*(10-2) = 72
  integer, parameter :: neqn = (nx-2)*(ny-2)
  integer, parameter :: cs = neqn+1
  double complex, dimension(neqn) :: us ! solution being iterated        
  double complex, dimension(neqn,cs) :: ut ! temporary for partial results 
  double complex, dimension(nx) :: xg
  data xg / (-1.0,-0.9), (-0.9,-0.8), (-0.8,-0.7), (-0.6,-0.4), &
            (-0.3,-0.2), ( 0.0, 0.0), ( 0.35,0.4), ( 0.65,0.7), &
            ( 0.85,0.9), ( 1.0, 1.2), ( 1.2, 1.4) /

  double complex, dimension(ny) :: yg
  data yg / (-1.1,-0.8), (-1.0,-0.7), (-0.9,-0.6), &
            (-0.6,-0.4), (-0.2,-0.1), ( 0.3,0.5),  &
            ( 0.8, 0.7), ( 1.0, 1.0), ( 1.2,1.4), (1.3,1.5) /

  double complex, dimension(nx) :: cxx    ! nuderiv coefficients 
  double complex, dimension(ny)  :: cyy
  double complex, dimension(nx,ny) :: cxy
  double complex, dimension(nx) :: cx
  double complex, dimension(ny)  :: cy
  double complex :: 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 complex, intent(in) :: x, y
      double complex :: value
    end function u
  end interface
  interface
    function c(x, y) result(term) ! from solve
      implicit none
      double complex, intent(in) :: x, y
      double complex :: term
    end function c
  end interface
  interface
    subroutine simeqC1(n, B, X)
      implicit none
      integer, intent(in) :: n
      double complex, dimension(n,n+1), intent(in) :: b
      double complex, dimension(n), intent(out) :: X
     end subroutine simeqC1
  end interface
  interface
     subroutine nuderivC(order, npoint, point, X, C)
     implicit none
     integer, intent(in) :: order
     integer, intent(in) :: npoint
     integer, intent(in) :: point
     double complex, dimension(npoint), intent(in) :: X
     double complex, dimension(npoint), intent(inout) :: C
     end subroutine nuderivC
  end interface
end module global

module util
  interface
    function S(i,j) result(k)
    implicit none
    integer, intent(in) :: i
    integer, intent(in) :: j
    integer :: k
    end function S
  end interface
end module util

module support
  interface
     subroutine check() ! typically not known, instrumentation 
     end subroutine check
  end interface

  interface
     subroutine printexact()
     end subroutine printexact
  end interface

  interface
     subroutine init_matrix()
     end subroutine init_matrix
  end interface

end module support

program pdeCnu22_eq1
  use global
  use util
  use support
  implicit none
  integer :: i, j
  
  print *, "pdeCnu22_eq1.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 complex grid on rectangle " 
  print *, "  (-1.0,-0.9) (-1.1,-0.8) to (1.2,1.4) (1.3,1.5) " 
  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(1)
  print *, "xmax=",xg(nx), "nx=",nx
  print *, "ymin=",yg(1)
  print *, "ymax=",yg(ny), "ny=",ny
  print *, "U(xg(1),yg(1))=",u(xg(1),yg(1))
  print *, "c(xg(1),yg(1))=",c(xg(1),yg(1))
  
  call init_matrix()
  print *, "matrix initialized"

  print "(/a)", "initial matrix, left side, upper"
  do i=1,4
    print *, (ut(i,j),j=1,4)
  end do
  print "(/a)", "initial matrix, right side, upper"
  do i=1,4 
    print *, (ut(i,j),j=neqn-3,cs)
  end do
  print "(/a)", "initial matrix, right side, lower"
  do i=neqn-3,cs
    print *, (ut(i,j),j=neqn-3,cs)
  end do
  print *, " "

  call simeqC1(neqn, ut, us)
  call printexact()

  call check()
  avgerror = error/neqn
  print "('error=',e15.8,' avg=',e15.8,' max=',e15.8)", &
           error, avgerror, maxerror
end program pdeCnu22_eq1 

function u(x, y) result(value) ! solution for boundary and check 
  implicit none
  double complex, intent(in) :: x, y
  double complex :: value
  value = x*x*x + 2.0*y*y*y + 3.0*x*x*y + 4.0*x*y*y + 5.0*x*x + &
	  6.0*y*y + 7.0*x + 8.0;
end function u

function c(x, y) result(term) ! from solve
  implicit none
  double complex, intent(in) :: x, y
  double complex :: term
  term = 6.0*x*x*x + 12.0*y*y*y + 18.0*x*x*y + 24.0*x*y*y + &
        57.0*x*x + 82.0*y*y + 64.0*x*y + 122.0*x + 114.0*y + 110.0
end function c

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

subroutine printexact()
  use global
  use util
  implicit none
  integer :: i, j
  
  print "(/a)", "xg, yg,exact solution u, computed us, error"
  do i=2,nx-1
    do j=2,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(S(i,j))
    print *, xg(i)
    print *, yg(j)
    print *, u(xg(i),yg(j))
    print *, us(S(i,j))
    print "('err= ',e15.7)", abs(u(xg(i),yg(j))-us(S(i,j)))
    print "(/a)"
    end do
  end do
end subroutine printexact 

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

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

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

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

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

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

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

      ! coefu * u(x,y) 
      coefu = 6.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
  
function S(i,j) result(k)
  use global
  implicit none
  integer, intent(in) :: i
  integer, intent(in) :: j
  integer :: k
  k = (i-2)+(nx-2)*(j-2)+1  ! i,j in us solution  
  if(k.lt.1 .or. k.gt.neqn) then
    print *, 'S(i,j)=k ', i, j, k
    k = 1
  end if
end function S