program gencanma

  use bmgencan, only: gencan, genunc
  use iso_c_binding, only: c_ptr, c_loc, c_f_pointer
  use mdgp, only: np,r,tau,x0,xgiven,inimdgp,endmdgp,procrustes

  implicit none

  type pdata_type
     integer :: counters(3) = 0
  end type pdata_type

  ! LOCAL SCALARS
  logical :: extallowed,hfixstr
  integer :: allocstat,flag,hnnzmax,ierr,indj,istop,j,iter,maxit,n
  real(kind=8) :: eps,f,finish,ftarget,gpsupn,start,timegencan,timeiniguess,xerror
  type(pdata_type), target :: pdata

  ! LOCAL ARRAYS
  character(len=80) :: probname
  real(kind=8), allocatable :: x(:),g(:)

  write(*,*) 'Enter problem name: '
  read(*,*) probname

  write(*,*) 'Initializing problem ...'

  call cpu_time(start)

  call inimdgp(probname,flag)
  if ( flag .ne. 0 ) then
     open(unit=10, file='gencan-tabline.txt')
     write(10,1010) probname,'inimdgp',flag
     close(10)
     stop
  end if

  call cpu_time(finish)

  timeiniguess = finish - start

  write(*,*) 'Solving optimization problem ...'

  n = np * r

  hnnzmax = 9 * ( tau / 2 ) + 6 * np

  allocate(x(n),g(n),stat=allocstat )
  if ( allocstat .ne. 0 ) then
     write(*,*) 'Error in allocate!'
     stop
  end if

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

  ftarget     =  1.0d-10
  eps         =  1.0d-08
  maxit       =    50000
  hfixstr     =   .true.
  extallowed  =   .true.

  call cpu_time(start)

  call genunc(evalf,evalg,evalh,hnnzmax,hfixstr,n,x,f,g,gpsupn, &
       ftarget,eps,maxit,extallowed,iter,ierr,istop,pdata=c_loc(pdata))

  call cpu_time(finish)

  timegencan = finish - start

  write(*,*) 'Saving solution and comparing with known solution (via Procrustes) ...'

  call procrustes(n,x,flag)
  if ( flag .ne. 0 ) then
     open(unit=10, file='gencan-tabline.txt')
     write(10,1010) probname,'procrustes',flag
     close(10)
     stop
  end if

  open(unit=10,file=trim(probname)//'-computedsol-gencan.txt')

  do j = 1,np
     indj = ( j - 1 ) * r
     write(10,*) x(indj+1:indj+r)
  end do

  close(10)

  open(unit=10,file=trim(probname)//'-givensol-transl.txt')

  do j = 1,np
     indj = ( j - 1 ) * r
     write(10,*) xgiven(indj+1:indj+r)
  end do

  close(10)
  
  xerror = 0.0d0
  do j = 1,np
     indj = ( j - 1 ) * r
     xerror = max( xerror, maxval( abs( x(indj+1:indj+r) - xgiven(indj+1:indj+r) ) ) / &
          max( 1.0d0, maxval( abs( xgiven(indj+1:indj+r) ) ) ) )
  end do

  write(*,*) 'Problem                                  = ',probname
  write(*,*) 'np                                       = ',np
  write(*,*) 'tau                                      = ',tau
  write(*,*) 'Density                                  = ',100.0d0 * tau / ( np ** 2 - np )
  write(*,*) 'CPU time for computing initial guess     = ',timeiniguess
  write(*,*) 'Number of iterations                     = ',iter
  write(*,*) 'Number of evaluations of f               = ',pdata%counters(1)
  write(*,*) 'Number of evaluations of gradient        = ',pdata%counters(2)
  write(*,*) 'Number of evaluations of Hessian         = ',pdata%counters(3)
  write(*,*) 'CPU time of gencan                       = ',timegencan
  write(*,*) 'ierr                                     = ',ierr
  write(*,*) 'istop                                    = ',istop
  write(*,*) 'Final objective function value           = ',f
  write(*,*) 'Final sup-norm of the gradient           = ',gpsupn
  write(*,*) 'Error with respect to the known solution = ',xerror

  open(unit=10, file='gencan-tabline.txt')
  write(10,1000) probname,np,tau,100.0d0*tau/(np**2-np),timeiniguess,iter,pdata%counters(1), &
       pdata%counters(2),pdata%counters(3),timegencan,ierr,istop,f,gpsupn,xerror
  close(10)

  deallocate(x,stat=allocstat)
  if ( allocstat .ne. 0 ) then
     write(*,*) 'Deallocation error in cdma.'
     stop
  end if

  call endmdgp(flag)

1000 format(A20,1X,I10,1X,I10,0P,F5.2,1X,0PF8.2,4(1X,I10),0PF8.2,1X,I10,1X,I10,1PD17.10,2(1X,1PD8.2))
1010 format(A20,1X,A20,1X,I10)

contains

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

  subroutine evalf(n,x,f,inform,pdataptr)

    use mdgp, only: np,r,tau,dcolsta,drow,dval

    implicit none

    ! SCALAR ARGUMENTS
    integer, intent(in) :: n
    integer, intent(inout) :: inform
    real(kind=8), intent(out) :: f
    type(c_ptr), optional, intent(in) :: pdataptr

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

    ! LOCAL SCALARS
    integer :: i,indi,indj,j
    type(pdata_type), pointer :: pdata

    ! LOCAL ARRAYS
    real(kind=8) :: tmp(r)

    call c_f_pointer(pdataptr,pdata)
    pdata%counters(1) = pdata%counters(1) + 1

    f = 0.0d0

    do j = 1,np
       indj = ( j - 1 ) * r

       do i = dcolsta(j),dcolsta(j+1) - 1
          if ( drow(i) .gt. j ) then
             indi = ( drow(i) - 1 ) * r
             tmp(1:r) = x(indi+1:indi+r) - x(indj+1:indj+r)
             f = f + ( dot_product( tmp(1:r), tmp(1:r) ) - dval(i) ) ** 2
          end if
       end do
    end do

    f = ( 2.0d0 / tau ) * f

  end subroutine evalf

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

  subroutine evalg(n,x,g,inform,pdataptr)

    use mdgp, only: np,r,tau,dcolsta,drow,dval

    implicit none

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

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

    ! LOCAL SCALARS
    integer :: i,indi,indj,j
    type(pdata_type), pointer :: pdata

    ! LOCAL ARRAYS
    real(kind=8) :: tmp(r)

    call c_f_pointer(pdataptr,pdata)
    pdata%counters(2) = pdata%counters(2) + 1

    g(1:n) = 0.0d0

    do j = 1,np
       indj = ( j - 1 ) * r

       do i = dcolsta(j),dcolsta(j+1) - 1
          indi = ( drow(i) - 1 ) * r
          tmp(1:r) = x(indj+1:indj+r) - x(indi+1:indi+r)
          g(indj+1:indj+r) = g(indj+1:indj+r) + ( dot_product( tmp(1:r), tmp(1:r) ) - dval(i) ) * tmp(1:r)
       end do
    end do

    g(1:n) = ( 8.0d0 / tau ) * g(1:n)

  end subroutine evalg

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

  subroutine evalh(n,x,lim,hnnz,hrow,hcol,hval,inform,pdataptr)

    use mdgp, only: np,r,tau,dcolsta,drow,dval

    implicit none

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

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

    ! LOCAL SCALARS
    integer :: i,indi,indj,indk,j,k,s,t
    type(pdata_type), pointer :: pdata

    ! LOCAL ARRAYS
    real(kind=8) :: tmp(r)

    call c_f_pointer(pdataptr,pdata)
    pdata%counters(3) = pdata%counters(3) + 1

    hnnz = 0

    do j = 1,np
       indj = ( j - 1 ) * r

       do s = 1,r
          do t = 1,s
             if ( indj + t .eq. indj + s ) then
                hnnz = hnnz + 1
                hrow(hnnz) = indj + s
                hcol(hnnz) = indj + s
                hval(hnnz) = 0.0d0
                do k = dcolsta(j),dcolsta(j+1) - 1
                   indk = ( drow(k) - 1 ) * r
                   tmp(1:r) = x(indk+1:indk+r) - x(indj+1:indj+r)
                   hval(hnnz) = hval(hnnz) + 2.0d0 * ( x(indj+s) - x(indk+s) ) ** 2 + &
                        dot_product( tmp(1:r), tmp(1:r) ) - dval(k)
                end do
                hval(hnnz) = ( 8.0d0 / tau ) * hval (hnnz)
             else if ( indj + t .lt. indj + s ) then
                hnnz = hnnz + 1
                hrow(hnnz) = indj + t
                hcol(hnnz) = indj + s
                hval(hnnz) = 0.0d0
                do k = dcolsta(j),dcolsta(j+1) - 1
                   indk = ( drow(k) - 1 ) * r
                   hval(hnnz) = hval(hnnz) + ( x(indj+t) - x(indk+t) ) * ( x(indj+s) - x(indk+s) )
                end do
                hval(hnnz) = ( 16.0d0 / tau ) * hval (hnnz)
             end if
          end do

          do i = dcolsta(j),dcolsta(j+1) - 1
             indi = ( drow(i) - 1 ) * r
             do t = 1,r
                if ( indi + t .eq. indj + s ) then
                   hnnz = hnnz + 1
                   hrow(hnnz) = indj + s
                   hcol(hnnz) = indj + s
                   hval(hnnz) = 0.0d0
                   do k = dcolsta(j),dcolsta(j+1) - 1
                      indk = ( drow(k) - 1 ) * r
                      tmp(1:r) = x(indk+1:indk+r) - x(indj+1:indj+r)
                      hval(hnnz) = hval(hnnz) + 2.0d0 * ( x(indj+s) - x(indk+s) ) ** 2 + &
                           dot_product( tmp(1:r), tmp(1:r) ) - dval(k)
                   end do
                   hval(hnnz) = ( 8.0d0 / tau ) * hval (hnnz)

                else if ( indi + t .lt. indj + s ) then
                   if ( drow(i) .eq. j ) then
                      hnnz = hnnz + 1
                      hrow(hnnz) = indi + t
                      hcol(hnnz) = indj + s
                      hval(hnnz) = 0.0d0
                      do k = dcolsta(j),dcolsta(j+1) - 1
                         indk = ( drow(k) - 1 ) * r
                         hval(hnnz) = hval(hnnz) + ( x(indj+t) - x(indk+t) ) * ( x(indj+s) - x(indk+s) )
                      end do
                      hval(hnnz) = ( 16.0d0 / tau ) * hval (hnnz)
                   else if ( s .eq. t ) then
                      hnnz = hnnz + 1
                      hrow(hnnz) = indi + t
                      hcol(hnnz) = indj + s
                      tmp(1:r) = x(indi+1:indi+r) - x(indj+1:indj+r)
                      hval(hnnz) = ( 2.0d0 / tau ) * ( - 8.0d0 * ( x(indj+s) - x(indi+s) ) ** 2 + &
                           dot_product( tmp(1:r), tmp(1:r) ) - dval(i) )
                   else
                      hnnz = hnnz + 1
                      hrow(hnnz) = indi + t
                      hcol(hnnz) = indj + s
                      hval(hnnz) = ( 2.0d0 / tau ) * ( 8.0d0 * ( x(indi+t) - x(indj+t) ) * ( x(indj+s) - x(indi+s) ) )
                   end if
                end if
             end do
          end do
       end do
    end do

  end subroutine evalh

end program gencanma
