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

program spgma

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

  ! LOCAL SCALARS
  integer :: allocstat,fcnt,flag,indj,inform,iprint,iter,j,maxfc,maxit,n,spginfo
  real(kind=8) :: epsopt,f,finish,ftarget,gpsupn,start,timecd,timeiniguess,xerror

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

  ! COMMON BLOCKS
  integer :: counters(2)
  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='spg-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

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

  x(1:n) = x0(1:n)
  
  counters(1:2) = 0

  iprint = 0
  
  maxit = huge( 1 )
  maxfc = huge( 1 )
  
  ftarget = 1.0d-10
  epsopt  = 1.0d-08
  
  call cpu_time(start)

  call spg(n,x,ftarget,epsopt,maxit,maxfc,iprint,f,gpsupn,iter,fcnt,spginfo,inform)
  
  call cpu_time(finish)

  timecd = 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='spg-tabline.txt')
     write(10,1010) probname,'procrustes',flag
     close(10)
     stop
  end if

  open(unit=10,file=trim(probname)//'-computedsol-spg.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               = ',counters(1)
  write(*,*) 'Number of evaluations of gradient        = ',counters(2)
  write(*,*) 'CPU time of SPG                          = ',timecd
  write(*,*) 'spginfo                                  = ',spginfo
  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='spg-tabline.txt')
  write(10,1000) probname,np,tau,100.0d0*tau/(np**2-np),timeiniguess,iter,counters(1),counters(2), &
       timecd,spginfo,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,1X,0PF5.2,1X,0PF8.2,3(1X,I10),1X,0PF8.2,1X,I10,1X,1PD17.10,2(1X,1PD8.2))
1010 format(A20,1X,A20,1X,I10)
  
end program spgma

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

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(2)
  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(2)
  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 proj(n,x,flag)

  implicit none

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

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

  flag = 0

end subroutine proj
