! THIS IS BRATU 3D

module bmproblem

  use iso_c_binding, only: c_ptr,c_f_pointer

  implicit none

  type :: pdata_type
     integer :: np
     integer :: counters(2)
     real(kind=8) :: h,theta
     real(kind=8), allocatable :: phi(:,:,:),uopt(:,:,:),x0(:)
  end type pdata_type

  public :: pdata_type,iniprob,endprob,evalr,evalj

contains

  ! ******************************************************************
  ! ******************************************************************

  subroutine iniprob(n,jnnz,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(out) :: jnnz,n
    type(c_ptr), optional, intent(in) :: pdataptr

    ! LOCAL SCALARS
    integer :: allocstat,i,j,k
    real(kind=8) :: a,b,c,tmp
    type(pdata_type), pointer :: pdata

    call c_f_pointer(pdataptr,pdata)

    pdata%counters(1:2) = 0
    
    write(*,*) 'Number of points (np): '
    read(*,*) pdata%np

    n = ( pdata%np - 2 ) ** 3

    pdata%h = 1.0d0 / ( pdata%np - 1 )

    jnnz = 0
    do i = 2,pdata%np - 1
       do j = 2,pdata%np - 1
          do k = 2,pdata%np - 1
             if ( k .ne. pdata%np - 1 ) then
                jnnz = jnnz + 1
             end if
             if ( k .ne. 2 ) then
                jnnz = jnnz + 1
             end if
             if ( j .ne. pdata%np - 1 ) then
                jnnz = jnnz + 1
             end if
             if ( j .ne. 2 ) then
                jnnz = jnnz + 1
             end if
             if ( i .ne. pdata%np - 1 ) then
                jnnz = jnnz + 1
             end if
             if ( i .ne. 2 ) then
                jnnz = jnnz + 1
             end if
             jnnz = jnnz + 1
          end do
       end do
    end do

    write(*,*) 'Theta: '
    read(*,*) pdata%theta

    allocate(pdata%uopt(pdata%np,pdata%np,pdata%np), &
         pdata%phi(2:pdata%np-1,2:pdata%np-1,2:pdata%np-1),pdata%x0(n),stat=allocstat)
    if ( allocstat .ne. 0 ) then
       write(*,*) 'Error in allocate statement!'
       stop
    end if

    ! Initial guess
    
    pdata%x0(1:n) = 0.0d0
    
    ! "Known solution" used to generate the problem data

    do k = 1,pdata%np
       do j = 1,pdata%np
          do i = 1,pdata%np
             a = real( pdata%np - k ) / ( pdata%np - 1 )
             b = real( pdata%np - j ) / ( pdata%np - 1 )
             c = real( pdata%np - i ) / ( pdata%np - 1 )
             pdata%uopt(i,j,k) = 10.0d0 * a * b * c &
                  * (1.0d0 - a) * (1.0d0 - b) * (1.0d0 - c) * exp( a ** 4.5d0 )
          end do
       end do
    end do

    ! Right hand side of the system

    tmp = 1.0d0 / pdata%h ** 2
    
    do k = 2,pdata%np - 1
       do j = 2,pdata%np - 1
          do i = 2,pdata%np - 1
             pdata%phi(i,j,k) = - tmp * ( pdata%uopt(i,j,k+1) + pdata%uopt(i,j,k-1) + pdata%uopt(i,j+1,k) + &
                  pdata%uopt(i,j-1,k) + pdata%uopt(i+1,j,k) + pdata%uopt(i-1,j,k) - 6.0d0 * pdata%uopt(i,j,k) ) &
                  + pdata%theta * exp( pdata%uopt(i,j,k) )
          end do
       end do
    end do
    
  end subroutine iniprob

  ! ******************************************************************
  ! ******************************************************************

  subroutine endprob(pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    type(c_ptr), optional, intent(in) :: pdataptr

    ! LOCAL SCALARS
    integer :: allocstat
    type(pdata_type), pointer :: pdata

    call c_f_pointer(pdataptr,pdata)

    deallocate(pdata%uopt,pdata%phi,pdata%x0,stat=allocstat)
    if ( allocstat .ne. 0 ) then
       write(*,*) 'Error in allocate statement!'
       stop
    end if
 
  end subroutine endprob
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine evalr(n,x,r,ierr,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n
    integer, intent(inout) :: ierr
    type(c_ptr), optional, intent(in) :: pdataptr

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n)
    real(kind=8), intent(out) :: r(n)

    ! LOCAL SCALARS
    integer :: i,j,k,ell
    real(kind=8) :: tmp,xi0,xi1,xi2,xi3,xi4,xi5,xi6
    type(pdata_type), pointer :: pdata
    
    call c_f_pointer(pdataptr,pdata)

    pdata%counters(1) = pdata%counters(1) + 1

    if ( any( isnan( x(1:n) ) ) ) then
       write(*,*) 'EVALR: Given point x contains NaN!'
       ierr = - 93
       return
    end if
    
    tmp = 1.0d0 / pdata%h ** 2
    
    ell = 0
    do i = 2,pdata%np - 1
       do j = 2,pdata%np - 1
          do k = 2,pdata%np - 1
             ell = ell + 1
                
             xi0 = x(xind(i,j,k,pdata%np))

             if ( k .eq. pdata%np - 1 ) then
                xi1 = pdata%uopt(i,j,k+1)
             else
                xi1 = x(xind(i,j,k+1,pdata%np))
             end if

             if ( k .eq. 2 ) then
                xi2 = pdata%uopt(i,j,k-1)
             else
                xi2 = x(xind(i,j,k-1,pdata%np))
             end if

             if ( j .eq. pdata%np - 1 ) then
                xi3 = pdata%uopt(i,j+1,k)
             else
                xi3 = x(xind(i,j+1,k,pdata%np))
             end if

             if ( j .eq. 2 ) then
                xi4 = pdata%uopt(i,j-1,k)
             else
                xi4 = x(xind(i,j-1,k,pdata%np))
             end if

             if ( i .eq. pdata%np - 1 ) then
                xi5 = pdata%uopt(i+1,j,k)
             else
                xi5 = x(xind(i+1,j,k,pdata%np))
             end if

             if ( i .eq. 2 ) then
                xi6 = pdata%uopt(i-1,j,k)
             else
                xi6 = x(xind(i-1,j,k,pdata%np))
             end if
             
             r(ell) = - tmp * ( xi1 + xi2 + xi3 + xi4 + xi5 + xi6 - 6.0d0 * xi0 ) &
                  + pdata%theta * exp( xi0 ) - pdata%phi(i,j,k)
          end do
       end do
    end do
    
    if ( any( isnan( r(1:n) ) ) ) then
       write(*,*) 'EVALR: Computed r contains NaN!'
       ierr = - 93
       return
    end if
    
  end subroutine evalr

  ! ******************************************************************
  ! ******************************************************************

  subroutine evalj(n,x,jsta,jlen,lim,jvar,jval,ierr,pdataptr)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: lim,n
    integer, intent(inout) :: ierr
    type(c_ptr), optional, intent(in) :: pdataptr

    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n)
    integer, intent(out) :: jsta(n),jlen(n),jvar(lim)
    real(kind=8), intent(out) :: jval(lim)

    ! LOCAL SCALARS
    integer :: i,i0,iini,iend,j,jlentmp,jnnz,k,ell
    real(kind=8) :: tmp
    type(pdata_type), pointer :: pdata

    ! LOCAL ARRAYS
    integer :: jvartmp(7)
    real(kind=8) :: jvaltmp(7)
    
    call c_f_pointer(pdataptr,pdata)

    pdata%counters(2) = pdata%counters(2) + 1

    if ( any( isnan( x(1:n) ) ) ) then
       write(*,*) 'EVALJ: Given point x contains NaN!'
       ierr = - 94
       return
    end if
    
    tmp = 1.0d0 / pdata%h ** 2
    
    jnnz = 0
    
    ell = 0
    do i = 2,pdata%np - 1
       do j = 2,pdata%np - 1
          do k = 2,pdata%np - 1
             ell = ell + 1
                
             jvartmp(1:6) = - 1

             if ( k .ne. pdata%np - 1 ) then
                jvartmp(1) = xind(i,j,k+1,pdata%np) ! i1
             end if

             if ( k .ne. 2 ) then
                jvartmp(2) = xind(i,j,k-1,pdata%np) ! i2
             end if

             if ( j .ne. pdata%np - 1 ) then
                jvartmp(3) = xind(i,j+1,k,pdata%np) ! i3
             end if

             if ( j .ne. 2 ) then
                jvartmp(4) = xind(i,j-1,k,pdata%np) ! i4
             end if

             if ( i .ne. pdata%np - 1 ) then
                jvartmp(5) = xind(i+1,j,k,pdata%np) ! i5
             end if

             if ( i .ne. 2 ) then
                jvartmp(6) = xind(i-1,j,k,pdata%np) ! i6
             end if

             i0 = xind(i,j,k,pdata%np)
             jvartmp(7) = i0

             jvaltmp(1:7) = (/ -tmp, -tmp, -tmp, -tmp, -tmp, -tmp, 6.0d0 * tmp + pdata%theta * exp( x(i0) ) /)

             jlentmp = count( jvartmp(1:7) .ne. -1 )
                
             jsta(ell) = jnnz + 1

             jlen(ell) = jlentmp
                
             jnnz = jnnz + jlen(ell)

             if ( jnnz .gt. lim ) then
                ierr = - 94
                return
             end if
                
             iini = jsta(ell)
             iend = jsta(ell) + jlen(ell) - 1
                
             jvar(iini:iend) = pack( jvartmp(1:7), jvartmp(1:7) .ne. - 1 )
             jval(iini:iend) = pack( jvaltmp(1:7), jvartmp(1:7) .ne. - 1 )
          end do
       end do
    end do
    
    if ( any( isnan( jval(1:jnnz) ) ) ) then
       write(*,*) 'EVALJ: Computed Jacobian contains NaN!'
       ierr = - 93
       return
    end if
    
  end subroutine evalj

  ! ******************************************************************
  ! ******************************************************************

  integer function xind(i,j,k,np)

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: i,j,k,np

    xind = ( np - 2 ) ** 2 * ( i - 2 ) + ( np - 2 ) * ( j - 2 ) + k - 1
    
  end function xind
  
end module bmproblem
