program optvorma use iso_c_binding, only: c_ptr,c_loc,c_f_pointer use optvor, only: pdata_type,evalfg,proj,approxg,drawsol,InA,compscale implicit none ! PARAMETERS real(kind=8), parameter :: PI = dacos( - 1.0d0 ), & oct_circ_rad = 1.0d0 / ( 2.0d0 * dsin(PI / 8.0d0) ), & dist_orig_to_bisec_segs = sqrt( ( 0.5d0 * ( oct_circ_rad * dcos( PI / 4.0d0 ) + & oct_circ_rad * dcos( PI / 2.0d0 ) ) ) ** 2.0d0 + & ( 0.5d0 * ( oct_circ_rad * dsin( PI / 4.0d0 ) + & oct_circ_rad * dsin( PI / 2.0d0 ) ) ) ** 2.0d0 ), & dist_orig_to_tipstar = dist_orig_to_bisec_segs + 2.0d0 * oct_circ_rad ! LOCAL SCALARS integer :: allocerr,fcnt,i,inform,iprint,iter,itrial,j,maxfc,maxit,maxnv,n, & nsites,ntrials,spginfo real(kind=8) :: area,epsopt,f,fbest,finish,ftarget,gpsupn,marg,scale,seed,start type(pdata_type), target :: pdata ! LOCAL ARRAYS real(kind=8) :: p(2) real(kind=8), allocatable :: g(:),x(:),x0(:),vol(:) character(len=80) :: filename ! FUNCTIONS real(kind=8) :: drand ! Problem data ! Set sites write(*,*) 'Enter nsites>0: ' read(*,*) nsites ! Set disjoint convex polygons Aj whose union define A pdata%npols = 9 allocate(pdata%nvpols(pdata%npols),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'Allocation error.' stop end if pdata%nvpols(1:pdata%npols) = (/ 9, 4, 4, 4, 4, 4, 4, 4, 4 /) 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 pdata%xpol(1, 1:pdata%nvpols(1)) = (/ oct_circ_rad * dcos(PI / 4.0d0), & oct_circ_rad * dcos(PI / 2.0d0), & oct_circ_rad * dcos(135.0d0 * PI / 180.0d0), & oct_circ_rad * dcos(PI), & oct_circ_rad * dcos(225.0d0 * PI / 180.0d0), & oct_circ_rad * dcos(3.0d0 * PI / 2.0d0), & oct_circ_rad * dcos(315.0d0 * PI / 180.0d0), & oct_circ_rad * dcos(2.0d0 * PI), & oct_circ_rad * dcos(PI / 4.0d0) /) pdata%ypol(1, 1:pdata%nvpols(1)) = (/ oct_circ_rad * dsin(PI / 4.0d0), & oct_circ_rad * dsin(PI / 2.0d0), & oct_circ_rad * dsin(135.0d0 * PI / 180.0d0), & oct_circ_rad * dsin(PI), & oct_circ_rad * dsin(225.0d0 * PI / 180.0d0), & oct_circ_rad * dsin(3.0d0 * PI / 2.0d0), & oct_circ_rad * dsin(315.0d0 * PI / 180.0d0), & oct_circ_rad * dsin(2.0d0 * PI), & oct_circ_rad * dsin(PI / 4.0d0) /) pdata%xpol(2, 1:pdata%nvpols(2)) = (/ oct_circ_rad * dcos(PI / 4.0d0), & oct_circ_rad * dcos(2.0d0 * PI), & dist_orig_to_tipstar * dcos(22.5d0 * PI / 180.0d0), & oct_circ_rad * dcos(PI / 4.0d0) /) pdata%ypol(2, 1:pdata%nvpols(2)) = (/ oct_circ_rad * dsin(PI / 4.0d0), & oct_circ_rad * dsin(2.0d0 * PI), & dist_orig_to_tipstar * dsin(22.5d0 * PI / 180.0d0), & oct_circ_rad * dsin(PI / 4.0d0) /) pdata%xpol(3, 1:pdata%nvpols(3)) = (/ oct_circ_rad * dcos(PI / 2.0d0), & oct_circ_rad * dcos(PI / 4.0d0), & dist_orig_to_tipstar * dcos(67.5d0 * PI / 180.0d0), & oct_circ_rad * dcos(PI / 2.0d0) /) pdata%ypol(3, 1:pdata%nvpols(3)) = (/ oct_circ_rad * dsin(PI / 2.0d0), & oct_circ_rad * dsin(PI / 4.0d0), & dist_orig_to_tipstar * dsin(67.5d0 * PI / 180.0d0), & oct_circ_rad * dsin(PI / 2.0d0) /) pdata%xpol(4, 1:pdata%nvpols(4)) = (/ oct_circ_rad * dcos(135.0d0 * PI / 180.0d0), & oct_circ_rad * dcos(PI / 2.0d0), & dist_orig_to_tipstar * dcos(112.5d0 * PI / 180.0d0), & oct_circ_rad * dcos(135.0d0 * PI / 180.0d0) /) pdata%ypol(4, 1:pdata%nvpols(4)) = (/ oct_circ_rad * dsin(135.0d0 * PI / 180.0d0), & oct_circ_rad * dsin(PI / 2.0d0), & dist_orig_to_tipstar * dsin(112.5d0 * PI / 180.0d0), & oct_circ_rad * dsin(135.0d0 * PI / 180.0d0) /) pdata%xpol(5, 1:pdata%nvpols(5)) = (/ oct_circ_rad * dcos(PI), & oct_circ_rad * dcos(135.0d0 * PI / 180.0d0), & dist_orig_to_tipstar * dcos(157.5d0 * PI / 180.0d0), & oct_circ_rad * dcos(PI) /) pdata%ypol(5, 1:pdata%nvpols(5)) = (/ oct_circ_rad * dsin(PI), & oct_circ_rad * dsin(135.0d0 * PI / 180.0d0), & dist_orig_to_tipstar * dsin(157.5d0 * PI / 180.0d0), & oct_circ_rad * dsin(PI) /) pdata%xpol(6, 1:pdata%nvpols(6)) = (/ oct_circ_rad * dcos(225.0d0 * PI / 180.0d0), & oct_circ_rad * dcos(PI), & dist_orig_to_tipstar * dcos(202.5d0 * PI / 180.0d0), & oct_circ_rad * dcos(225.0d0 * PI / 180.0d0) /) pdata%ypol(6, 1:pdata%nvpols(6)) = (/ oct_circ_rad * dsin(225.0d0 * PI / 180.0d0), & oct_circ_rad * dsin(PI), & dist_orig_to_tipstar * dsin(202.5d0 * PI / 180.0d0), & oct_circ_rad * dsin(225.0d0 * PI / 180.0d0) /) pdata%xpol(7, 1:pdata%nvpols(7)) = (/ oct_circ_rad * dcos(3.0d0 * PI / 2.0d0), & oct_circ_rad * dcos(225.0d0 * PI / 180.0d0), & dist_orig_to_tipstar * dcos(247.5d0 * PI / 180.0d0), & oct_circ_rad * dcos(3.0d0 * PI / 2.0d0) /) pdata%ypol(7, 1:pdata%nvpols(7)) = (/ oct_circ_rad * dsin(3.0d0 * PI / 2.0d0), & oct_circ_rad * dsin(225.0d0 * PI / 180.0d0), & dist_orig_to_tipstar * dsin(247.5d0 * PI / 180.0d0), & oct_circ_rad * dsin(3.0d0 * PI / 2.0d0) /) pdata%xpol(8, 1:pdata%nvpols(8)) = (/ oct_circ_rad * dcos(315.0d0 * PI / 180.0d0), & oct_circ_rad * dcos(3.0d0 * PI / 2.0d0), & dist_orig_to_tipstar * dcos(292.5d0 * PI / 180.0d0), & oct_circ_rad * dcos(315.0d0 * PI / 180.0d0) /) pdata%ypol(8, 1:pdata%nvpols(8)) = (/ oct_circ_rad * dsin(315.0d0 * PI / 180.0d0), & oct_circ_rad * dsin(3.0d0 * PI / 2.0d0), & dist_orig_to_tipstar * dsin(292.5d0 * PI / 180.0d0), & oct_circ_rad * dsin(315.0d0 * PI / 180.0d0) /) pdata%xpol(9, 1:pdata%nvpols(9)) = (/ oct_circ_rad * dcos(2.0d0 * PI), & oct_circ_rad * dcos(315.0d0 * PI / 180.0d0), & dist_orig_to_tipstar * dcos(337.5d0 * PI / 180.0d0), & oct_circ_rad * dcos(2.0d0 * PI) /) pdata%ypol(9, 1:pdata%nvpols(9)) = (/ oct_circ_rad * dsin(2.0d0 * PI), & oct_circ_rad * dsin(315.0d0 * PI / 180.0d0), & dist_orig_to_tipstar * dsin(337.5d0 * PI / 180.0d0), & oct_circ_rad * dsin(2.0d0 * PI) /) pdata%poledges(1,1:pdata%nvpols(1)-1) = (/ 1, 1, 1, 1, 1, 1, 1, 1 /) pdata%poledges(2,1:pdata%nvpols(2)-1) = (/ 1, 0, 0 /) pdata%poledges(3,1:pdata%nvpols(3)-1) = (/ 1, 0, 0 /) pdata%poledges(4,1:pdata%nvpols(4)-1) = (/ 1, 0, 0 /) pdata%poledges(5,1:pdata%nvpols(5)-1) = (/ 1, 0, 0 /) pdata%poledges(6,1:pdata%nvpols(6)-1) = (/ 1, 0, 0 /) pdata%poledges(7,1:pdata%nvpols(7)-1) = (/ 1, 0, 0 /) pdata%poledges(8,1:pdata%nvpols(8)-1) = (/ 1, 0, 0 /) pdata%poledges(9,1:pdata%nvpols(9)-1) = (/ 1, 0, 0 /) pdata%xmin = - dist_orig_to_tipstar - 1.0d0 pdata%xmax = dist_orig_to_tipstar + 1.0d0 pdata%ymin = - dist_orig_to_tipstar - 1.0d0 pdata%ymax = dist_orig_to_tipstar + 1.0d0 pdata%c(1:2) = 0.5d0 * (/ pdata%xmax + pdata%xmin, pdata%ymax + pdata%ymin /) pdata%rmax = norm2(pdata%c - (/ pdata%xmin - 0.5d0, pdata%ymin - 0.5d0 /)) pdata%drawscale = 0.8d0 ! Scale A so cells have ideal area near unity call compscale(nsites,scale,area,c_loc(pdata)) pdata%volA = area do j = 1,pdata%npols pdata%xpol(j,1:pdata%nvpols(j)) = scale * pdata%xpol(j,1:pdata%nvpols(j)) pdata%ypol(j,1:pdata%nvpols(j)) = scale * pdata%ypol(j,1:pdata%nvpols(j)) end do 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)) ) do j = 2,pdata%npols pdata%xmin = min( pdata%xmin, minval( pdata%xpol(j,1:pdata%nvpols(j)) ) ) pdata%xmax = max( pdata%xmax, maxval( pdata%xpol(j,1:pdata%nvpols(j)) ) ) pdata%ymin = min( pdata%ymin, minval( pdata%ypol(j,1:pdata%nvpols(j)) ) ) pdata%ymax = max( pdata%ymax, maxval( pdata%ypol(j,1:pdata%nvpols(j)) ) ) end do 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 = 0.75d0 / scale ! Goal function data pdata%Jweight(0:5) = 0.0d0 pdata%Jweight(0) = 10.0d0 pdata%Jweight(1) = 1.0d0 pdata%J0delta = 0.1d0 ! Optimize write(*,*) 'Enter ntrials>0: ' read(*,*) ntrials n = 2 * nsites allocate(x(n),x0(n),g(n),vol(nsites),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'Allocation error.' stop end if 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) epsopt = 1.0d-08 ftarget = 1.0d-08 maxit = 10000 maxfc = 10000 iprint = 0 call cpu_time(start) call spg(n,x,epsopt,ftarget,maxit,maxfc,iprint,f,gpsupn,iter,fcnt,spginfo,inform, & evalfg,proj,c_loc(pdata)) call cpu_time(finish) write(*,*) 'Statistics: f = ',f,' gpsupn = ',gpsupn,' it = ',iter,' fcnt = ',fcnt, & ' spginfo = ',spginfo,' time = ',finish-start if ( f .lt. fbest ) then fbest = f write(*,*) 'Solution was improved. New best f = ',fbest write(filename,"(A13,I0,A3)") 'fig-star-ini-',nsites,'.mp' call drawsol(n,x0,filename,c_loc(pdata)) write(filename,"(A13,I0,A3)") 'fig-star-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,gpsupn,iter,fcnt,spginfo,finish-start close(10) if ( fbest .le. ftarget ) then write(*,*) 'A global minimizer has been found!' exit end if end if end do deallocate(x,x0,g,vol,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'Deallocation error.' stop end if stop 1000 format(1P,E12.5,1X,1P,E12.5,1X,I3,1X,I6,1X,I6,1X,1P,E12.5,1X,1P,E7.1,2(1X,I6),1X,I2,0P,F14.6) end program optvorma