program bmbunchma

  use mdgp, only: np,r,tau,x0,xgiven,inimdgp,endmdgp,procrustes
  
  implicit none

  ! LOCAL SCALARS
  integer :: allocstat,fcnt,flag,hnnz,inform,indj,j,iter,maxfcnt,maxit,n
  real(kind=8) :: eps,f,finish,ftarget,gnorm,start,timebmbunch,timeiniguess,xerror
  
  ! LOCAL ARRAYS
  character(len=80) :: probname
  integer, allocatable :: hcol(:),hrow(:)
  real(kind=8), allocatable :: x(:),g(:)

  ! COMMON BLOCKS
  integer :: counters(3)
  common /cntdata/ counters

  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='bmbunch-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

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

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

  x(1:n) = x0(1:n)
  
  counters(1:3) = 0
  
  eps     = 1.0d-08
  ftarget = 1.0d-10
  maxit   = 10000
  maxfcnt = 10 * maxit
  
  call gethpat(hnnz,hrow,hcol)
  
  call cpu_time(start)

  call CRNewton(n,x,f,g,gnorm,eps,ftarget,maxit,maxfcnt,hnnz,hrow,hcol,iter,fcnt,inform)

  call cpu_time(finish)
        
  timebmbunch = 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='bmbunch-tabline.txt')
     write(10,1010) probname,'procrustes',flag
     close(10)
     stop
  end if

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

  do j = 1,np
     indj = ( j - 1 ) * r
     write(10,*) x(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               = ',counters(1)
  write(*,*) 'Number of evaluations of gradient        = ',counters(2)
  write(*,*) 'Number of evaluations of Hessian         = ',counters(3)
  write(*,*) 'CPU time of bmbunch                      = ',timebmbunch
  write(*,*) 'inform                                   = ',inform
  write(*,*) 'Final objective function value           = ',f
  write(*,*) 'Final sup-norm of the gradient           = ',gnorm
  write(*,*) 'Error with respect to the known solution = ',xerror
  
  open(unit=10, file='bmbunch-tabline.txt')
  write(10,1000) probname,np,tau,100.0d0*tau/(np**2-np),timeiniguess,iter,counters(1),counters(2), &
       counters(3),timebmbunch,inform,f,gnorm,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,1PD17.10,2(1X,1PD8.2))
1010 format(A20,1X,A20,1X,I10)

end program bmbunchma

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

subroutine gethpat(hnnz,hrow,hcol)

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

  implicit none

  !SCALAR ARGUMENTS
  integer, intent(in) :: hnnz

  ! ARRAY ARGUMENTS
  integer, intent(out) :: hrow(hnnz),hcol(hnnz)

 ! LOCAL SCALARS
 integer :: cont,i,indi,indj,j,s,t

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

    do s = 1,r
       do t = s,r
          if ( indj + t .eq. indj + s ) then
             cont = cont + 1
             hrow(cont) = indj + s
             hcol(cont) = indj + s
          else if ( indj + t .gt. indj + s ) then
             cont = cont + 1
             hrow(cont) = indj + t
             hcol(cont) = indj + s
          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
                cont = cont + 1
                hrow(cont) = indj + s
                hcol(cont) = indj + s
             else if ( indi + t .gt. indj + s ) then
                if ( drow(i) .eq. j ) then
                   cont = cont + 1
                   hrow(cont) = indi + t
                   hcol(cont) = indj + s
                else if ( s .eq. t ) then
                   cont = cont + 1
                   hrow(cont) = indi + t
                   hcol(cont) = indj + s
                else
                   cont = cont + 1
                   hrow(cont) = indi + t
                   hcol(cont) = indj + s
                end if
             end if
          end do
       end do
    end do
 end do

 if ( cont .ne. hnnz ) then
    write(*,*) 'There is something wrong in gethpat.'
    write(*,*) 'hnnz = ',hnnz,' cont = ',cont
    stop
 end if
 
end subroutine gethpat

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

subroutine evalf(n,x,f,flag)

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

  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)

  ! COMMON BLOCKS
  integer :: counters(3)
  common /cntdata/ counters

  ! LOCAL SCALARS
  integer :: i,indi,indj,j

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

  counters(1) = counters(1) + 1

  flag = 0
  
  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,flag)

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

  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)

  ! COMMON BLOCKS
  integer :: counters(3)
  common /cntdata/ counters

  ! LOCAL SCALARS
  integer :: i,indi,indj,j
  
  ! LOCAL ARRAYS
  real(kind=8) :: tmp(r)

  counters(2) = counters(2) + 1

  flag = 0
  
  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,hnnz,hrow,hcol,hval,flag)

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

  implicit none

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

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

 ! COMMON BLOCKS
 integer :: counters(3)
 common /cntdata/ counters

 ! LOCAL SCALARS
 integer :: cont,i,indi,indj,indk,j,k,s,t

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

 counters(3) = counters(3) + 1

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

    do s = 1,r
       do t = s,r
          if ( indj + t .eq. indj + s ) then
             cont = cont + 1
             if ( hrow(cont) .ne. indj + s .or. hcol(cont) .ne. indj + s ) then
                flag = - 1
                return
             end if
             hval(cont) = 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(cont) = hval(cont) + 2.0d0 * ( x(indj+s) - x(indk+s) ) ** 2 + &
                     dot_product( tmp(1:r), tmp(1:r) ) - dval(k)
             end do
             hval(cont) = ( 8.0d0 / tau ) * hval (cont)
          else if ( indj + t .gt. indj + s ) then
             cont = cont + 1
             if ( hrow(cont) .ne. indj + t .or. hcol(cont) .ne. indj + s ) then
                flag = - 1
                return
             end if
             hval(cont) = 0.0d0 
             do k = dcolsta(j),dcolsta(j+1) - 1
                indk = ( drow(k) - 1 ) * r
                hval(cont) = hval(cont) + ( x(indj+t) - x(indk+t) ) * ( x(indj+s) - x(indk+s) )
             end do
             hval(cont) = ( 16.0d0 / tau ) * hval (cont)
          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
                cont = cont + 1
                if ( hrow(cont) .ne. indj + s .or. hcol(cont) .ne. indj + s ) then
                   flag = - 1
                   return
                end if
                hval(cont) = 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(cont) = hval(cont) + 2.0d0 * ( x(indj+s) - x(indk+s) ) ** 2 + &
                        dot_product( tmp(1:r), tmp(1:r) ) - dval(k)
                end do
                hval(cont) = ( 8.0d0 / tau ) * hval (cont)
                
             else if ( indi + t .gt. indj + s ) then
                if ( drow(i) .eq. j ) then
                   cont = cont + 1
                   if ( hrow(cont) .ne. indi + t .or. hcol(cont) .ne. indj + s ) then
                      flag = - 1
                      return
                   end if
                   hval(cont) = 0.0d0
                   do k = dcolsta(j),dcolsta(j+1) - 1
                      indk = ( drow(k) - 1 ) * r
                      hval(cont) = hval(cont) + ( x(indj+t) - x(indk+t) ) * ( x(indj+s) - x(indk+s) )
                   end do
                   hval(cont) = ( 16.0d0 / tau ) * hval (cont)
                else if ( s .eq. t ) then
                   cont = cont + 1
                   if ( hrow(cont) .ne. indi + t .or. hcol(cont) .ne. indj + s ) then
                      flag = - 1
                      return
                   end if
                   tmp(1:r) = x(indi+1:indi+r) - x(indj+1:indj+r)
                   hval(cont) = ( 2.0d0 / tau ) * ( - 8.0d0 * ( x(indj+s) - x(indi+s) ) ** 2 + &
                        dot_product( tmp(1:r), tmp(1:r) ) - dval(i) )
                else
                   cont = cont + 1
                   if ( hrow(cont) .ne. indi + t .or. hcol(cont) .ne. indj + s ) then
                      flag = - 1
                      return
                   end if
                   hval(cont) = ( 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
