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

program algencanma

  use iso_c_binding, only: c_ptr,c_loc,c_f_pointer
  use optvor, only: pdata_type,compJ,drawsol,InA,compscale
  
  implicit none

  ! PARAMETERS
  integer, parameter :: nvert = 20
  real(kind=8), parameter :: PI = acos( - 1.0d0 ) 

  ! LOCAL SCALARS
  logical :: checkder
  integer :: allocerr,hnnzmax,i,inform,itrial,jcnnzmax,m,maxnv, &
       n,nsites,ntrials,nvparam
  real(kind=8) :: ang,area,cnorm,efacc,efstain,eoacc,eostain,epsfeas,epsopt, &
       f,fbest,finish,marg,nlpsupn,scale,seed,snorm,start
  character(len=80) :: filename,specfnm,outputfnm
  type(pdata_type), target :: pdata

  ! LOCAL ARRAYS
  logical :: coded(11)
  real(kind=8) :: p(2)
  logical, pointer :: equatn(:),linear(:)
  real(kind=8), allocatable :: l(:),lambda(:),u(:),x(:),x0(:)
  character(len=80) :: vparam(10)

  ! FUNCTIONS
  real(kind=8) :: drand
  
  ! EXTERNAL SUBROUTINES
  !external :: myevalf,myevalg,myevalh,myevalc,myevaljac,myevalhc, &
  !            myevalfc,myevalgjac,myevalgjacp,myevalhl,myevalhlp

  ! Set problem data

  ! Number of sites

  write(*,*) 'Enter nsites>0: '
  read(*,*) nsites

  ! Set disjoint convex polygons Aj whose union define A (a single Aj in this case)
  
  pdata%npols = 1
  
  allocate(pdata%nvpols(pdata%npols),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error.'
     stop
  end if

  pdata%nvpols(1:pdata%npols) = (/ nvert + 1 /)
  
  maxnv = maxval( pdata%nvpols(1:pdata%npols) )
  
  allocate(pdata%poledges(pdata%npols,maxnv),pdata%xpol(pdata%npols,maxnv), &
       pdata%ypol(pdata%npols,maxnv),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error.'
     stop
  end if

  do i = 1,nvert + 1
     ang = 2.0d0 * PI * ( i - 1 ) / nvert
     pdata%xpol(1,i) = cos( ang )
     pdata%ypol(1,i) = sin( ang )
  end do
  
  pdata%poledges(1,1:nvert) = 0

  ! Scale A so cells have ideal area near unity
  
  call compscale(nsites,scale,area,c_loc(pdata))

  pdata%volA = area 
  pdata%xpol(1,1:pdata%nvpols(1)) = scale * pdata%xpol(1,1:pdata%nvpols(1))
  pdata%ypol(1,1:pdata%nvpols(1)) = scale * pdata%ypol(1,1:pdata%nvpols(1))
  
  pdata%xmin = minval( pdata%xpol(1,1:pdata%nvpols(1)) )
  pdata%xmax = maxval( pdata%xpol(1,1:pdata%nvpols(1)) )
  pdata%ymin = minval( pdata%ypol(1,1:pdata%nvpols(1)) )
  pdata%ymax = maxval( pdata%ypol(1,1:pdata%nvpols(1)) )
  
  marg = 0.1d0 * max( 1.0d0, abs( pdata%xmin ), abs( pdata%xmax ), abs( pdata%ymin), abs( pdata%ymax ) )

  pdata%xmin = pdata%xmin - marg
  pdata%xmax = pdata%xmax + marg
  pdata%ymin = pdata%ymin - marg
  pdata%ymax = pdata%ymax + marg
  
  pdata%c(1:2) = 0.5d0 * (/ pdata%xmax + pdata%xmin, pdata%ymax + pdata%ymin /)
  pdata%rmax = norm2(pdata%c - (/ pdata%xmin, pdata%ymin /))
  pdata%drawscale = 2.5d0 / scale

  if ( nsites .ge. 10000 ) then
     pdata%drawscale = 5.0d0 / scale
  end if

  ! Coded subroutines

  coded(1:11) = .false.
  coded(7:8) = .true. ! fcsub,gjacsub

  ! Upper bounds on the number of sparse-matrices non-null elements

  jcnnzmax = huge( 1 )
  hnnzmax = huge( 1 )

  ! Checking derivatives?

  checkder = .true.

  ! Parameters setting

  epsfeas = 1.0d-08
  epsopt = 1.0d-08

  efstain = sqrt( epsfeas )
  eostain = epsopt ** 1.5d0

  efacc = sqrt( epsfeas )
  eoacc = sqrt( epsopt )

  outputfnm = ''
  specfnm = ''

  nvparam = 0

  ! Set objective function and constraints scheme

  n = 2 * nsites
  
  allocate(x0(n),x(n),l(n),u(n),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error in main program.'
     stop
  end if

  l(1:n) = - 1.0d+20
  u(1:n) =   1.0d+20

  l(1:n) = pdata%c - pdata%rmax
  u(1:n) = pdata%c + pdata%rmax

  pdata%Jweight(0)   =   0.0d0
  pdata%Jweight(1)   =   1.0d0
  pdata%Jweight(2)   =   0.0d0
  pdata%Jweight(3:4) = - 1.0d0
  pdata%Jweight(5)   =   0.0d0
  
  m = count( pdata%Jweight(3:5) .lt. 0.0d0 ) 

  allocate(equatn(m),linear(m),lambda(m),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error in main program'
     stop
  end if

  equatn(1:m) = .false.
  linear(1:m) = .false.

  lambda(1:m) = 0.0d0

  ! Optimize
  
  write(*,*) 'Enter ntrials>0: '
  read(*,*) ntrials

  fbest = huge( 1.0d0 )
  
  do itrial = 1,ntrials
     seed = 123456.0d0 * itrial

     write(*,*) 'Processing trial itrial = ',itrial
     
     do i = 1,nsites
        p(1) = pdata%xmin + ( pdata%xmax - pdata%xmin ) * drand(seed)
        p(2) = pdata%ymin + ( pdata%ymax - pdata%ymin ) * drand(seed)

        do while ( .not. InA(p,c_loc(pdata)) )
           p(1) = pdata%xmin + ( pdata%xmax - pdata%xmin ) * drand(seed)
           p(2) = pdata%ymin + ( pdata%ymax - pdata%ymin ) * drand(seed)
        end do
       
        x(2*i-1) = p(1)
        x(2*i)   = p(2)
     end do

     x0(1:n) = x(1:n)

     call cpu_time(start)
  
     call algencan(myevalf,myevalg,myevalh,myevalc,myevaljac,myevalhc, &
          myevalfc,myevalgjac,myevalgjacp,myevalhl,myevalhlp,jcnnzmax, &
          hnnzmax,epsfeas,epsopt,efstain,eostain,efacc,eoacc,outputfnm, &
          specfnm,nvparam,vparam,n,x,l,u,m,lambda,equatn,linear,coded, &
          checkder,f,cnorm,snorm,nlpsupn,inform)
     
     call cpu_time(finish)

     write(*,*) 'Statistics: f = ',f,' cnorm = ',cnorm,' snorm = ',snorm,' nlpsupn = ',nlpsupn, &
          ' inform = ',inform,' time = ',finish-start
     
     if ( cnorm .le. epsfeas .and. f .lt. fbest ) then
        fbest = f

        write(*,*) 'Solution was improved. New best f = ',fbest

        write(filename,"(A15,I0,A3)") 'fig-regpol-ini-',nsites,'.mp'
        call drawsol(n,x0,filename,c_loc(pdata))
  
        write(filename,"(A15,I0,A3)") 'fig-regpol-end-',nsites,'.mp'
        call drawsol(n,x,filename,c_loc(pdata))

        open(unit=10,file='tabline.txt')   
        write(10,1000) scale,pdata%volA,pdata%npols,nsites,itrial,f,cnorm,snorm,nlpsupn,inform,finish-start
        close(10)
     end if
  end do
  
  deallocate(x0,x,l,u,lambda,equatn,linear,stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Deallocation error in main program.'
     stop
  end if

  stop

1000 format(1P,E12.5,1X,1P,E12.5,1X,I3,1X,I6,1X,I6,1X,1P,E12.5,3(1X,1P,E7.1),1X,I2,0P,F14.6)

contains
  
! ******************************************************************
! ******************************************************************

subroutine myevalf(n,x,f,flag)

  implicit none

  ! SCALAR ARGUMENTS
  integer,      intent(in)  :: n
  integer,      intent(out) :: flag
  real(kind=8), intent(out) :: f

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

  flag = - 1

end subroutine myevalf

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

subroutine myevalg(n,x,g,flag)

  implicit none

  ! SCALAR ARGUMENTS
  integer,      intent(in)  :: n
  integer,      intent(out) :: flag

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

  flag = - 1

end subroutine myevalg

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

subroutine myevalh(n,x,hrow,hcol,hval,hnnz,lim,lmem,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical,      intent(out) :: lmem
  integer,      intent(in)  :: lim,n
  integer,      intent(out) :: flag,hnnz

  ! ARRAY ARGUMENTS
  integer,      intent(out) :: hcol(lim),hrow(lim)
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: hval(lim)

  flag = - 1

end subroutine myevalh

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

subroutine myevalc(n,x,ind,c,flag)

  implicit none

  ! SCALAR ARGUMENTS
  integer,      intent(in)  :: ind,n
  integer,      intent(out) :: flag
  real(kind=8), intent(out) :: c

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

  flag = - 1

end subroutine myevalc

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

subroutine myevaljac(n,x,ind,jcvar,jcval,jcnnz,lim,lmem,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical, intent(out) :: lmem
  integer, intent(in)  :: ind,lim,n
  integer, intent(out) :: flag,jcnnz

  ! ARRAY ARGUMENTS
  integer,      intent(out) :: jcvar(lim)
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: jcval(lim)

  flag = - 1

end subroutine myevaljac

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

subroutine myevalhc(n,x,ind,hcrow,hccol,hcval,hcnnz,lim,lmem,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical, intent(out) :: lmem
  integer, intent(in)  :: ind,lim,n
  integer, intent(out) :: flag,hcnnz

  ! ARRAY ARGUMENTS
  integer,      intent(out) :: hccol(lim),hcrow(lim)
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: hcval(lim)

  flag = - 1

end subroutine myevalhc

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

subroutine myevalfc(n,x,f,m,c,flag)

  implicit none

  ! SCALAR ARGUMENTS
  integer,      intent(in)  :: m,n
  integer,      intent(out) :: flag
  real(kind=8), intent(out) :: f

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

  ! LOCAL SCALARS
  integer :: ind,k
  
  ! LOCAL ARRAYS
  real(kind=8) :: farray(0:5),garray(n,0:5)
  
  flag = 0

  call compJ(n,x,farray,garray,flag,c_loc(pdata))

  ind = 0
  f = 0.0d0
  
  do k = 0,5
     if ( pdata%Jweight(k) .gt. 0.0d0 ) then
        f = f + pdata%Jweight(k) * farray(k)
     else if ( pdata%Jweight(k) .lt. 0.0d0 ) then
        ind = ind + 1
        c(ind) = farray(k) - 0.01d0
     end if
  end do
  
end subroutine myevalfc

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

subroutine myevalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical,      intent(out) :: lmem
  integer,      intent(in)  :: lim,m,n
  integer,      intent(out) :: flag,jcnnz

  ! ARRAY ARGUMENTS
  integer,      intent(out) :: jcfun(lim),jcvar(lim)
  real(kind=8), intent(in)  :: x(n)
  real(kind=8), intent(out) :: g(n),jcval(lim)

  ! LOCAL SCALARS
  integer :: i,ind,k
  
  ! LOCAL ARRAYS
  real(kind=8) :: farray(0:5),garray(n,0:5)
  
  flag = 0
  lmem = .false.

  call compJ(n,x,farray,garray,flag,c_loc(pdata))

  ind = 0
  jcnnz = 0
  g(1:n) = 0.0d0
  
  do k = 0,5
     if ( pdata%Jweight(k) .gt. 0.0d0 ) then
        g(1:n) = g(1:n) + pdata%Jweight(k) * garray(1:n,k)
     else if ( pdata%Jweight(k) .lt. 0.0d0 ) then
        if ( jcnnz + n .gt. lim ) then
           lmem = .true.
           return
        end if
        
        ind = ind + 1
        jcfun(jcnnz+1:jcnnz+n) = ind
        jcvar(jcnnz+1:jcnnz+n) = (/ (i,i=1,n) /)
        jcval(jcnnz+1:jcnnz+n) = garray(1:n,k)
        jcnnz = jcnnz + n
     end if
  end do
  
end subroutine myevalgjac

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

subroutine myevalgjacp(n,x,g,m,p,q,work,gotj,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical,   intent(inout) :: gotj
  integer,   intent(in)    :: m,n
  integer,   intent(out)   :: flag
  character, intent(in)    :: work

  ! ARRAY ARGUMENTS
  real(kind=8), intent(in)    :: x(n)
  real(kind=8), intent(inout) :: p(m),q(n)
  real(kind=8), intent(out)   :: g(n)

  flag = - 1

end subroutine myevalgjacp

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

subroutine myevalhl(n,x,m,lambda,sf,sc,hlrow,hlcol,hlval,hlnnz,lim,lmem,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical,      intent(out) :: lmem
  integer,      intent(in)  :: lim,m,n
  integer,      intent(out) :: flag,hlnnz
  real(kind=8), intent(in)  :: sf

  ! ARRAY ARGUMENTS
  integer,      intent(out) :: hlcol(lim),hlrow(lim)
  real(kind=8), intent(in)  :: lambda(m),sc(m),x(n)
  real(kind=8), intent(out) :: hlval(lim)

  flag = - 1

end subroutine myevalhl

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

subroutine myevalhlp(n,x,m,lambda,sf,sc,p,hp,goth,flag)

  implicit none

  ! SCALAR ARGUMENTS
  logical,      intent(inout) :: goth
  integer,      intent(in)    :: m,n
  integer,      intent(out)   :: flag
  real(kind=8), intent(in)    :: sf

  ! ARRAY ARGUMENTS
  real(kind=8), intent(in)  :: lambda(m),p(n),sc(m),x(n)
  real(kind=8), intent(out) :: hp(n)

  flag = - 1

end subroutine myevalhlp

end program algencanma

