! pdenu23_eq.f90 3D non equal h boundary value PDE  non iterative    
!  u(x,y,z) = x^3  + y^3  +z^3  + 2x^2y + 3xy^2 + 4z^2x + 5y + 6  
!     du/dx = 3x^2 + 0    + 0   + 4xy   + 3y^2  + 4z^2  + 0  + 0  
! d^2u/dx^2 = 6x   + 0    + 0   + 4y    + 0     + 0     + 0  + 0  
!     du/dy = 0    + 3y^2 + 0   + 2x^2  + 6xy   + 0     + 5  + 0  
! d^2u/dy^2 = 0    + 6y   + 0   + 0     + 6x    + 0     + 0  + 0  
!     du/dz = 0    + 0    +3z^2 + 0     + 0     + 8zx   + 0  + 0  
! d^2u/dz^2 = 0    + 0    +6z   + 0     + 0     + 8x    + 0  + 0  
!                                                                 
! solve d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 = c(x,y,z)              
! by second order method on cube -1,-1,-1 to 1.2,1,1              

! second order numerical difference                              
! d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 =                            
! 1/hx^2(u(x-hx,y,z)+u(x+hx,y,z)-2u(x,y,z)) +                    
! 1/hy^2(u(x,y-hy,z)+u(x,y+hy,z)-2u(x,y,z)) +                    
! 1/hz^2(u(x,y,z-hz)+u(x,y,z+hz)-2u(x,y,z)) + O(h^2)             
!  = c(x,y,z) = 12x + 10y - 8z                                   

! solve system of linear equations for                           
! u(i-1,j,k)/hx^2 + u(i+1,j,k)/hx^2 +                            
! u(i,j-1,k)/hy^2 + u(i,j+1,k)/hy^2 +                            
! u(i,j,k-1)/hz^2 + u(i,j,k+1)/hz^2 -                            
! (2/hx^2 + 2/hy^2 + 2/hz^2)*u(i,j,k) = c(x,y,z)                 
! for i=1,nx-2 j=1,ny-2 k=1,nz-2 while substituting              
! u(i,j,k) boundary values i=0 or nx-1, j=0 or ny-1, k=0 or nz-1 
! x=xmin+i*hx, y=ymin+j*hz, z=zmin+k*hz                          

! create two dimensional array with                              
! (nx-2)*(ny-2)*(nz-2) columns, one for each equation            
! and (nx-2)*(ny-2)*(nz-2)+1 rows for variables plus constant    
! Initialize the matrix, see init_matrix=the source code      
! Then solve the system of linear equations to get the solution. 

module global  
  implicit none
  integer :: nx = 10 ! one less than 'C' when Fortran subscript starts at zero
  integer :: ny = 10
  integer :: nz = 8
  double precision, dimension(0:566) :: us ! solution being iterated        
  double precision, dimension(0:566,0:567) :: ut ! temporary for partial results 
                 ! nx-1 by ny-1 by nz-1 internal, non boundary cells 
  double precision :: xmin = -1.0
  double precision :: ymin = -1.0
  double precision :: zmin = -1.0
  double precision :: xmax =  1.2
  double precision :: ymax =  1.0
  double precision :: zmax =  1.0
  double precision :: hx   ! x step = (xmax-xmin)/nx 
  double precision :: hy   ! y step = (ymax-ymin)/ny 
  double precision :: hz   ! z step = (zmax-zmin)/nz
  double precision :: error         ! sum of absolute exact-computed         
  double precision :: avgerror      ! average error                          
  double precision :: maxerror      ! maximum cell error                     
  interface
    function u(x, y, z) result(value) ! solution for boundary and check 
      implicit none
      double precision, intent(in) :: x, y, z
      double precision :: value
    end function u
  end interface
  interface
    function c(x, y, z) result(term) ! from solve
      implicit none
      double precision, intent(in) :: x, y, z
      double precision :: term
    end function c
  end interface
  interface
    subroutine simeq(n, B, X)
      implicit none
      integer, intent(in) :: n
      double precision, dimension(0:566,0:567), intent(in) :: b
      double precision, dimension(0:566), intent(out) :: X
     end subroutine simeq
  end interface
end module global

program pdenu23_eq
  use global
  implicit none
  integer :: i, j
  
  print *, "pdenu23_eq.f90 running"
  hx = (xmax-xmin)/nx ! step size=x 
  hy = (ymax-ymin)/ny ! step size=y 
  hz = (zmax-zmin)/nz ! step size=z 
  print *, "xmin=",xmin, "xmax=",xmax, "nx=",nx, "hx=",hx
  print *, "ymin=",ymin, "ymax=",ymax, "ny=",ny, "hy=",hy
  print *, "zmin=",zmin, "zmax=",zmax, "nz=",nz, "hz=",hz
  
  call init_matrix()
  print *, "matrix initialized"
  print "(/a)", "initial matrix, left side, upper"
  do i=0,8
    print "(9f8.3)", (ut(i,j),j=0,8)
  end do
  print "(/a)", "initial matrix, right side, upper"
  do i=0,8 
    print "(9f8.3)", (ut(i,j),j=559,567)
  end do
  print "(/a)", "initial matrix, right side, lower"
  do i=558,566
    print "(9f8.3)", (ut(i,j),j=559,567)
  end do
  print *, " "


  call simeq(566, ut, us)
  call printexact()
  call printus()
  call printdiff()
  call check()
  avgerror = error/((nx-1)*(ny-1)*(nz-1))
  print "('error=',e15.8,' avg=',e15.8,' max=',e15.8)", &
           error, avgerror, maxerror
end program pdenu23_eq 

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

function c(x, y, z) result(term) ! from solve
  implicit none
  double precision, intent(in) :: x, y, z
  double precision :: term
  term = 20.0*x + 10.0*y + 6.0*z
end function c

subroutine check() ! typically not known, instrumentation 
  use global
  implicit none
  integer :: i, j, k
  double precision :: zexact, err
  
  error = 0.0
  maxerror = 0.0
  do i=1,nx-1
    do j=1,ny-1
      do k=1,nz-1
        zexact = u(xmin+i*hx,ymin+j*hy,zmin+k*hz)
        err = abs(us((i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1))-zexact)
        error = error + err
        if(err>maxerror) maxerror=err
      end do
    end do
  end do
end subroutine check 

subroutine printus()
  use global
  implicit none
  integer :: i, j, k
  
  print "(/a)", "computed solution"
  do k=1,nz-1
    print *, "u(x,y,",k,") plane"
    do i=1,nx-1
      print "(11f7.3)", (us((i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1)), j=1,min(ny,10))
    end do
  end do
end subroutine printus 

subroutine printexact()
  use global
  implicit none
  integer :: i, j, k
  
  print "(/a)", "exact solution"
  do k=0,nz
    print *, "u(x,y,",k,") plane"
    do i=0,nx
      print "(11f7.3)", (u(xmin+i*hx,ymin+j*hy,zmin+k*hz),j=0,min(ny,10))
    end do
  end do
end subroutine printexact 

subroutine printdiff()
  use global
  implicit none
  integer :: i, j, k
  
  print "(/a)", "difference exact-computed solution" 
  do k=1,nz-1
    print *, "diff(x,y,",k,") plane"
    do i=1,nx-1
      print "(10f8.4)", (u(xmin+i*hx,ymin+j*hy,zmin+k*hz)- &
            us((i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1)), j=1,min(ny-1,10))
    end do
  end do
end subroutine printdiff 

subroutine init_matrix()
  use global
  implicit none
  integer :: i, j, k, ii, jj, kk, cs 
  double precision ::  hx2, hy2, hz2
  
  hx2 = hx*hx
  hy2 = hy*hy
  hz2 = hz*hz
  cs = (nx-1)*(ny-1)*(nz-1) ! constant column subscript 
  
  ! zero internal cells 
  do i=1,nx-1
    do j=1,ny-1
      do k=1,nz-1
        do ii=1,nx-1
          do jj=1,ny-1
            do kk=1,nz-1 
              ut(((i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1)), &
                 ((ii-1)+(nx-1)*(jj-1)+(nx-1)*(ny-1)*(kk-1))) = 0.0
            end do
          end do
        end do
      end do
    end do
  end do

  print *, "internal cells zeroed "          
  do i=1,nx-1  ! inside boundary 
    do j=1,ny-1
      do k=1,nz-1
        ii = (i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1) ! column of matrix 
        ut(ii,ii) = -(2.0/hx2 + 2.0/hy2 + 2.0/hz2)
        ut(ii,cs) = c(xmin+i*hx, ymin+j*hy, zmin+k*hz)
        if((i+1)<(nx-1)) then
          ut(ii,((i-1)+1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1)) = 1.0/hx2
        else
          ut(ii,cs) = ut(ii,cs) - u(xmin+(i+1)*hx, ymin+j*hy, zmin+k*hz)/hx2
        end if
        if((i-1).ne.0) then
          ut(ii,((i-1)-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*(k-1)) = 1.0/hx2
        else
          ut(ii,cs) = ut(ii,cs) - u(xmin+(i-1)*hx, ymin+j*hy, zmin+k*hz)/hx2
        end if
        if((j+1)<(ny-1)) then
          ut(ii,(i-1)+(nx-1)*((j-1)+1)+(nx-1)*(ny-1)*(k-1)) = 1.0/hy2
        else
          ut(ii,cs) = ut(ii,cs) - u(xmin+i*hx, ymin+(j+1)*hy, zmin+k*hz)/hy2
        end if
        if((j-1).ne.0) then
          ut(ii,(i-1)+(nx-1)*((j-1)-1)+(nx-1)*(ny-1)*(k-1)) = 1.0/hy2
        else
          ut(ii,cs) = ut(ii,cs) - u(xmin+i*hx, ymin+(j-1)*hy, zmin+k*hz)/hy2
        end if
        if((k+1)<(nz-1)) then
          ut(ii,(i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*((k-1)+1)) = 1.0/hz2
        else
          ut(ii,cs) = ut(ii,cs) - u(xmin+i*hx, ymin+j*hy, zmin+(k+1)*hz)/hz2
        end if
        if((k-1).ne.0) then 
          ut(ii,(i-1)+(nx-1)*(j-1)+(nx-1)*(ny-1)*((k-1)-1)) = 1.0/hz2
        else
          ut(ii,cs) = ut(ii,cs) - u(xmin+i*hx, ymin+j*hy, zmin+(k-1)*hz)/hz2
        end if
      end do
    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:566,0:567), intent(inout) :: B
  double precision, dimension(0:566), 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