program ivmmahess C IVM Test driver. This master file uses IVM (Inexact Variable C Metric Method) to solve some the 3D-lOcation problem. C C R. Andreani, E. G. Birgin, J. M. Martínez and J. Yuan, "Spectral C projected gradient and variable metric methods for optimization C with linear inequalities", IMA Journal of Numerical Analysis 25, C pp. 221-252, 2005. C C This version 25 SET 2003 by E.G.Birgin, J.M.Martinez and J.Yuan. C Final revision ... by E.G.Birgin, J.M.Martinez and J. Yuan. C PARAMETERS integer mmax parameter (mmax = 3 650 000) integer nmax parameter (nmax = 4 000) integer anzemax parameter (anzemax = 3 * mmax) C COMMON SCALARS integer anze,n C COMMON ARRAYS integer arow(anzemax),acol(anzemax) double precision aval(anzemax) C LOCAL SCALARS logical extrap,hess integer fcnt,flag,gcnt,hcnt,i,innit,j,kprint,m,maxfc,maxinnit, + maxoutit,modh,outit,output,p,ntpmin,ntpmax,actconstr, + inform double precision amax,dsupn,deucn2,epsdsn,epsden,epsmac,epsnlsn, + eta,f,gamma,infty,lmin,lmax,steabs,sterel,tau,tau1,tau2, + theta,tol,prob,rmin,rmax,bmin,fsol,diff character * 10 pname real time,dtime C LOCAL ARRAYS integer ngp(3) double precision b(mmax),l(nmax),u(nmax),x(nmax),y(mmax),step(3), + sol(nmax) real dum(2) C EXTERNAL SUBROUTINES external ivmhess C INTRINSIC FUNCTIONS intrinsic dtime C COMMON BLOCKS common /ldata1/ arow,acol,aval,anze,n C SET PROBLEM DATA pname = 'LOCATI3D' ngp(1) = 11 ngp(2) = 11 ngp(3) = 11 step(1) = 1.0d0 step(2) = 1.0d0 step(3) = 1.0d0 prob = 0.01d0 rmin = 1.0d0 rmax = 1.0d0 ntpmin = 200 ntpmax = 200 read(*,*) prob,ntpmin,ntpmax c prob = 0.0015d0 c ntpmin = 100 c ntpmax = 100 if ( prob .eq. 0.0d0 ) then stop end if call genpro(ngp,step,prob,rmin,rmax,ntpmin,ntpmax,m,n,arow,acol, + aval,anze,b,l,u,x,sol) call evalf(n,sol,fsol,inform) write(*,*) 'There is a feasible point where f = ',fsol C INCLUDE THE BOUND CONSTRAINTS IN THE FORMAT A x <= b if ( m + 2 * n .gt. mmax ) then write(*,*) 'Increase mmax and run again' stop end if do i = 1,n m = m + 1 anze = anze + 1 arow(anze) = m acol(anze) = i aval(anze) = - 1.0d0 b(m) = - l(i) m = m + 1 anze = anze + 1 arow(anze) = m acol(anze) = i aval(anze) = 1.0d0 b(m) = u(i) end do if ( anze .gt. anzemax ) then write(*,*) 'Increase anzemax and run again' stop end if write(*,*) 'PROBLEM DATA' write(*,*) 'Name : 3D Location Problem' write(*,*) 'Number of variables : ',n write(*,*) 'Number of constraints of type A x <= b : ',m write(*,*) 'Number of non-zero elements in matrix : ',anze C SET THE INITIAL ESTIMATION OF THE LAGRANGE MULTIPLIERS do i = 1,m y(i) = 0.0d0 end do C SETUP OPTIMIZER PARAMETERS call iniivm(epsdsn,epsden,epsnlsn,maxoutit,maxinnit,maxfc,p, + extrap,hess,output,kprint,eta,theta,tau,tol, + lmin,lmax,gamma,sterel,steabs,epsmac,infty,tau1, + tau2) C CALL THE OPTIMIZER time= dtime(dum) call ivmhess(n,m,x,y,arow,acol,aval,anze,b,epsdsn,epsden,epsnlsn, + maxoutit,maxinnit,maxfc,p,extrap,hess,output,kprint, + eta,theta,tau,tol,lmin,lmax,gamma,sterel,steabs, + epsmac,infty,tau1,tau2,f,dsupn,deucn2,amax,outit, + innit,fcnt,gcnt,hcnt,modh,flag) time= dtime(dum) time= dum(1) C WRITE STATISTICS write(*,*) write(*,*) 'Solution (f = ',f,' Known feasible point (f =', + fsol,')' diff = 0.0d0 do i = 1,n/3 write(*,fmt=9010) (x(3*(i-1)+j),j=1,3),(sol(3*(i-1)+j),j=1,3) do j = 1,3 diff = max( diff, abs(x(3*(i-1)+j)-sol(3*(i-1)+j)) ) end do end do write(*,*) 'Difference = ',diff C TEST FEASIBILITY ( b - Ax >= 0 ) do i = 1,anze b(arow(i)) = b(arow(i)) - aval(i) * x(acol(i)) end do actconstr = 0 bmin = 1.0d+99 do i = 1,m bmin = min( bmin, b(i) ) if ( b(i) .le. 1.0d-6 ) then actconstr = actconstr + 1 end if end do write(*,*) write(*,*) 'If this is positive, the final point is feasible',bmin write(*,*) 'Number of active constraints = ',actconstr open(30, file='ivmline.txt') write(30,fmt=9000) pname,n,m,outit,fcnt,gcnt,hcnt,modh,innit,time, + f,dsupn,sqrt(deucn2),amax,flag,actconstr close(30) stop C NON EXECUTABLE STATEMENTS 9000 format(A8,' &',I7,' &',I7,' &',I6,' &',I6,' &',I6,' &',I6,' &', + I6,' &',I6,' &',F9.2,' &',1PD14.6,' &',1PD9.1,' &',1PD9.1, + ' &',1PD9.1,' &',I1,' &',I7) 9010 format(6(1X,1P,D12.4)) end C ****************************************************** C ****************************************************** subroutine evalf(n,x,f,inform) C This subroutine computes the objective function. C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C point at which the function will be evaluated. C C On Return C C f double precision, C function value at x, C C inform integer, C termination parameter: C 0= the function was successfully evaluated, C 1= some error occurs in the function evaluation. C SCALAR ARGUMENTS double precision f integer n,inform C ARRAY ARGUMENTS double precision x(n) C LOCAL SCALARS double precision dist integer i,j,npol C INTRINSIC FUNCTIONS intrinsic sqrt inform = 0 npol = n / 3 - 1 f = 0.0d0 do i = 1,npol dist = 0.0d0 do j = 1,3 dist = dist + (x(3 * npol + j) - x(3 * (i - 1) + j)) ** 2 end do dist = sqrt( dist ) if ( dist .le. 1.0d-4 ) then write (*,fmt=*) + 'ERROR IN PROBLEM DEFINITION (DIST TOO SMALL)' inform = 1 end if f = f + dist end do f = f / npol end C ****************************************************** C ****************************************************** subroutine evalg(n,x,g,inform) C This subroutine computes the gradient of the C objective function. C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C point at which the gradient will be evaluated. C C On Return C C g double precision g(n), C gradient vector at x, C C inform integer, C termination parameter: C 0= the gradient was successfully evaluated, C 1= some error occurs in the gradient evaluation. C SCALAR ARGUMENTS integer n,inform C ARRAY ARGUMENTS double precision g(n),x(n) C LOCAL SCALARS double precision dist integer i,j,npol C LOCAL ARRAYS double precision diff(3) C INTRINSIC FUNCTIONS intrinsic sqrt inform = 0 do j = 1,n g(j) = 0.0d0 end do npol = n / 3 - 1 do i = 1,npol dist = 0.0d0 do j = 1,3 diff(j) = x(3 * npol + j) - x(3 * (i - 1) + j) dist = dist + diff(j) ** 2 end do dist = sqrt( dist ) do j = 1,3 g(3 * (i - 1) + j) = - ( diff(j) / dist ) / npol g(3 * npol + j) = g(3 * npol + j) - g(3 * (i - 1) + j) end do end do end C ****************************************************** C ****************************************************** subroutine evalH(n,x,ldh,H,inform) C This subroutine computes the hessian of the C objective function. C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C point at which the Hessian will be evaluated, C C lda integer, C leading dimension of matrix H. C C On Return C C H double precision H(nmax,n), C Hessian matrix at x, C C inform integer, C termination parameter: C 0= the Hessian was successfully evaluated, C 1= some error occurs in the Hessian evaluation. C SCALAR ARGUMENTS integer n,ldh,inform C ARRAY ARGUMENTS double precision H(ldh,n),x(n) C LOCAL SCALARS integer i,j,k,npol,ind1,ind2,ind3,ind4 double precision dist,dist2 C LOCAL ARRAYS double precision diff(3) inform = 0 do i = 1,n do j = 1,n H(i,j) = 0.0d0 end do end do npol = n / 3 - 1 do i = 1,npol dist2 = 0.0d0 do j = 1,3 diff(j) = x(3 * npol + j) - x(3 * (i - 1) + j) dist2 = dist2 + diff(j) ** 2 end do dist = sqrt( dist2 ) * npol do j = 1,3 ind1 = 3 * (i - 1) + j do k = 1,3 ind2 = 3 * (i - 1) + k H(ind1,ind2) = - diff(j) * diff(k) / dist / dist2 if ( j .eq. k ) then H(ind1,ind2) = H(ind1,ind2) + 1.0d0 / dist end if ind3 = 3 * npol + j ind4 = 3 * npol + k H(ind1,ind4) = - H(ind1,ind2) H(ind3,ind2) = - H(ind1,ind2) H(ind3,ind4) = H(ind3,ind4) + H(ind1,ind2) end do end do end do end C ****************************************************** C ****************************************************** subroutine genpro(ngp,step,prob,rmin,rmax,ntpmin,ntpmax,m,n, + arow,acol,aval,anze,b,l,u,x,sol) C This subroutine generates a 3D-location problem. First, a regular C grid with nx, ny and nz points at edges x, y and z, respectively, C in the positive orthant is considered. The points of the grid start C at the origin with an horizontal distances of xstep, ystep and zstep, C respectively. This grid will be the space at which the cities C (represented by polygons) will be distributed. Before building the C cities, an area (rectangular polygon) of preservation where C almost nothing can be done, is defined. This area of preservation C will receive, after the construction of the cities, an hydraulic C plant of energy generation (to supply the energy to the cities). Then, C in the rest of the space, the cities are built. At each point of the C grid (out of the central region) a city (represented by a polygon) C will be built with probability prob. The definition of the polygon C uses variables ntpmin, ntpmax, rmin and rmax in a way described in C the genpol (generate polygon) subroutine. To transmit the energy C from the plant to the cities, a tower inside each city and a tower C inside the central region must be built. The objective of the C problem is to determine the location of this towers in order C to minimize the sum of the distances from each city tower to the C central one. C C On Entry: C C ngp integer ngp(3), C number grid points at each one of the 3 dimensions, C C step double precision step(3), C distance between points of the grid at each dimension, C C prob double precision, C probability of defining a city at a point of the grid C (0 <= prob <= 1), C C ntpmin integer, C minimum number of half spaces for each polygon, C C ntpmax integer, C maximum number of half spaces for each polygon, C C rmin double precision, C minimum ratio of the sphere contained inside each polygon, C C rmax double precision, C maximum ratio of the sphere contained inside each polygon. C C On output: C C m integer, C number of constraints, C C n integer, C number of variables, C C arow integer arow(mmax*nmax), C acol integer acol(mmax*nmax), C aval double precision aval(mmax*nmax), C coefficients matrix of the constraints in sparse format, C C anze integer, C number of nonzero elements of the coefficients matrix, C C b double precision b(mmax), C right hand side of the linear inequality constraints. C C x double precision x(nmax), C an interior feasible point. C PARAMETERS integer mmax parameter (mmax = 3 650 000) integer nmax parameter (nmax = 4 000) integer anzemax parameter (anzemax = 3 * mmax) C SCALAR ARGUMENTS double precision prob,rmax,rmin integer m,n,ntpmin,ntpmax,anze C ARRAY ARGUMENTS integer ngp(3),arow(anzemax),acol(anzemax) double precision step(3),aval(anzemax),b(mmax),l(nmax),u(nmax), + x(nmax),sol(nmax) C LOCAL SCALARS double precision r,seed integer i,j,k,nac,q C LOCAL ARRAYS double precision c(3),cbox(3),cpl(3),cpu(3) C EXTERNAL FUNCTIONS double precision drand external drand C EXTERNAL SUBROUTINES external genpol C INITIALIZE NUMBER OF VARIABLES AND CONSTRAINTS m = 0 n = 0 C DEFINE CENTRAL POLYGON. THE SIX CONSTRAINTS WHICH DEFINE THE C CENTRAL POLYGON ARE, IN FACT, BOUND CONSTRAINTS. do i = 1,3 C The lower bound for component i cpl(i) = 0.4d0 * ( ngp(i) - 1 ) * step(i) C The upper bound for component i cpu(i) = 0.6d0 * ( ngp(i) - 1 ) * step(i) C The center of the central polygon cbox(i) = ( cpl(i) + cpu(i) ) / 2.0d0 end do C SEED FOR THE RANDOM GENERATION seed = 760013 C DEFINE CITY-POLYGONS CENTERED AT THE GRID POINTS AND 'SAFELY' C OUTSIDE THE CENTRAL REGION. SAFELY MEANS rmax FAR FROM THE C CENTRAL REGION. anze = 0 do i = 0,ngp(1) - 1 c(1) = i * step(1) do j = 0,ngp(2) - 1 c(2) = j * step(2) do k = 0,ngp(3) - 1 c(3) = k * step(3) if ( c(1) .le. cpl(1) - 2.0d0 * rmax .or. + c(1) .ge. cpu(1) + 2.0d0 * rmax .or. + c(2) .le. cpl(2) - 2.0d0 * rmax .or. + c(2) .ge. cpu(2) + 2.0d0 * rmax .or. + c(3) .le. cpl(3) - 2.0d0 * rmax .or. + c(3) .ge. cpu(3) + 2.0d0 * rmax ) then C GENERATE A NEW POLYGON if ( drand( seed ) .le. prob ) then c if ( c + c + ( c(1).eq.2.0d0.and.c(2).eq.5.0d0.and.c(3).eq.5.0d0 ) .or. c + ( c(1).eq.8.0d0.and.c(2).eq.5.0d0.and.c(3).eq.5.0d0 ) .or. c + ( c(1).eq.5.0d0.and.c(2).eq.2.0d0.and.c(3).eq.5.0d0 ) .or. c + ( c(1).eq.5.0d0.and.c(2).eq.8.0d0.and.c(3).eq.5.0d0 ) .or. c + ( c(1).eq.5.0d0.and.c(2).eq.5.0d0.and.c(3).eq.2.0d0 ) .or. c + ( c(1).eq.5.0d0.and.c(2).eq.5.0d0.and.c(3).eq.8.0d0 ) c + c + ) then C SET THE COMPONENTS OF THE FEASIBLE INTERIOR POINT if ( n + 3 .gt. nmax ) then write(*,*) 'Increase nmax and run again' stop end if do q = 1,3 x(n+q) = c(q) end do C GENERATE THE RATIO OF THE SPHERE r = rmin + ( rmax - rmin ) * drand( seed ) C GENERATE THE SIX CONSTRAINTS WHICH DETERMINE C THE SMALLEST BOX WHICH CONTAINS THE SPHERE OF C RATIO r CENTERED AT c. THIS SIX CONSTRAINTS C ARE, IN FACT, BOUND CONSTRAINTS do q = 1,3 C The lower bound for component i l(n+q) = c(q) - r C The upper bound for component i u(n+q) = c(q) + r end do C GENERATE THE OTHER CONSTRAINTS WHICH DEFINE C THE POLYGON call genpol(cbox,c,r,ntpmin,ntpmax,seed,m,n, + arow,acol,aval,anze,b,nac,sol) C INCREASE THE NUMBER OF VARIABLES AND CONSTRAINTS m = m + nac n = n + 3 end if end if end do end do end do C DEFINE CENTRAL POLYGON. THE SIX CONSTRAINTS WHICH DEFINE THE C CENTRAL POLYGON ARE, IN FACT, BOUND CONSTRAINTS. do i = 1,3 C The lower bound for component i l(n+i) = cpl(i) C The upper bound for component i u(n+i) = cpu(i) x(n+i) = ( l(n+i) + u(n+i) ) / 2.0d0 sol(n+i) = x(n+i) end do C INCREASE THE NUMBER OF VARIABLES n = n + 3 end C ****************************************************** C ****************************************************** subroutine genpol(cbox,c,r,ntpmin,ntpmax,seed,m,n,arow,acol,aval, + anze,b,nac,sol) C This subroutine generates a polygon in R^3 as the intersection C of half spaces determined by planes tangent to the sphere with C center c and ratio r. The number of tangent planes is randomly C generated satisfying ntpmin <= number of tangent planes <= ntpmax. C C On Entry: C C c double precision c(3), C center of the sphere, C C r double precision, C ratio of the sphere, C C ntpmin integer, C minimum number of tangent planes, C C ntpmax integer, C maximum number of tangent planes, C C seed double precision, C seed for the random generation, C C m integer, C current number of constraints, C C n integer, C current number of variables. C C arow integer arow(mmax*nmax), C acol integer acol(mmax*nmax), C aval double precision aval(mmax*nmax), C coefficients matrix of the current constraints in C sparse format, C C anze integer, C current number of non zero elements of matrix A, C C b double precision b(mmax), C right hand side of the current linear inequality C constraints A x <= b. C C On Return: C C arow integer arow(mmax*nmax), C acol integer acol(mmax*nmax), C aval double precision aval(mmax*nmax), C coefficients matrix of the current constraints in C sparse format, C C anze integer, C current number of non zero elements of matrix A, C C b double precision b(mmax), C right hand side of the current linear inequality C constraints A x <= b, C C nac integer, C number of added constraints. C PARAMETERS integer mmax parameter (mmax = 3 650 000) integer nmax parameter (nmax = 4 000) integer anzemax parameter (anzemax = 3 * mmax) C SCALAR ARGUMENTS integer ntpmin,ntpmax,m,n,anze,nac double precision r,seed C ARRAY ARGUMENTS integer arow(anzemax),acol(anzemax) double precision cbox(3),c(3),aval(anzemax),b(mmax),sol(nmax) C LOCAL SCALARS double precision penorm,rhs integer i,j C LOCAL ARRAYS double precision p(3) C EXTERNAL FUNCTIONS double precision drand external drand C INTRINSIC FUNCTIONS intrinsic int,sqrt C GENERATE THE NUMBER OF TANGENT PLANES nac = ntpmin + int( ( ntpmax - ntpmin + 1 ) * drand( seed ) ) if ( m + nac .gt. mmax ) then write(*,*) 'Increase mmax and run again' stop end if C GENERATE THE TANGENT PLANES C GENERATE A SPECIFIC PLANE if ( nac .ge. 1 ) then penorm = 0.0d0 do j = 1,3 p(j) = cbox(j) - c(j) penorm = penorm + p(j) ** 2 end do penorm = sqrt( penorm ) rhs = r ** 2 do j = 1,3 p(j) = p(j) / penorm * r rhs = rhs + p(j) * c(j) end do do j = 1,3 sol(n+j) = c(j) + p(j) end do do j = 1,3 anze = anze + 1 if ( anze .gt. anzemax ) then write(*,*) 'Increase anzemax and run again' stop end if arow(anze) = m + 1 acol(anze) = n + j aval(anze) = p(j) end do end if C SET THE RIGHT HAND SIDE VECTOR b b(m+1) = rhs do i = 2,nac C GENERATE A POINT p IN THE BORDER OF THE SPHERE OF RATIO r C CENTERED AT THE ORIGIN. THE PLANE TANGENT TO THAT POINT, C TRANSLATED TO THE SPHERE OF CENTER c IS DEFINED BY C p^t x = p^t ( p + c ) = r^2 + p^t c. THE HALF SPACE WHICH C CONTAINS c IS p^t x <= r^2 + p^t c. penorm = 0.0d0 do j = 1,3 p(j) = drand(seed) penorm = penorm + p(j) ** 2 end do penorm = sqrt( penorm ) rhs = r ** 2 do j = 1,3 p(j) = p(j) / penorm * r rhs = rhs + p(j) * c(j) end do C ADD THE COEFFICIENT TO THE SPARSE MATRIX A do j = 1,3 anze = anze + 1 if ( anze .gt. anzemax ) then write(*,*) 'Increase anzemax and run again' stop end if arow(anze) = m + i acol(anze) = n + j aval(anze) = p(j) end do C SET THE RIGHT HAND SIDE VECTOR b b(m+i) = rhs end do end c **************************************************** c **************************************************** double precision function drand(ix) double precision ix double precision a,p,b15,b16,xhi,xalo,leftlo,fhi,k data a/16807.d0/,b15/32768.d0/,b16/65536.d0/,p/2147483647.d0/ xhi= ix/b16 xhi= xhi - dmod(xhi,1.d0) xalo= (ix-xhi*b16)*a leftlo= xalo/b16 leftlo= leftlo - dmod(leftlo,1.d0) fhi= xhi*a + leftlo k= fhi/b15 k= k - dmod(k,1.d0) ix= (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k if (ix.lt.0) ix= ix + p drand= ix*4.656612875d-10 return end c ****************************************************** c ****************************************************** subroutine iniivm(epsdsn,epsden,epsnlsn,maxoutit,maxinnit,maxfc,p, + extrap,hess,output,kprint,eta,theta,tau,tol, + lmin,lmax,gamma,sterel,steabs,epsmac,infty,tau1, + tau2) C This subroutine set the IVM parameters with its default values C SCALAR ARGUMENTS logical hess,extrap integer kprint,maxfc,maxinnit,maxoutit,output,p double precision epsdsn,epsden,epsnlsn,eta,theta,tau,tol,lmin, + lmax,gamma,tau1,tau2,sterel,steabs,epsmac,infty epsdsn = 1.0d-06 epsden = 1.0d-06 epsnlsn = 1.0d-08 maxoutit = 100 maxinnit = 100000 maxfc = 2000 p = 10 extrap = .true. hess = .true. output = 3 kprint = 5 c eta = 0.68d0 c theta = 0.85d0 eta = 0.90d0 theta = 0.95d0 tau = 1.0d-01 tol = 0.99d0 lmin = 1.0d-03 lmax = 1.0d+03 gamma = 1.0d-04 sterel = 1.0d-07 steabs = 1.0d-10 infty = 1.0d+20 call mcheps(epsmac) tau1 = epsmac ** ( 1.0d0 / 3.0d0 ) tau2 = epsmac ** ( 1.0d0 / 3.0d0 ) end c ****************************************************** c ****************************************************** subroutine mcheps(epsmac) C This subroutine computes an approximation of machine epsilon. C SCALAR ARGUMENTS double precision epsmac epsmac = 1.0d0 20 if ( ( 1.0d0 + epsmac ) .ne. 1.0d0 ) then epsmac = epsmac / 2.0d0 go to 20 end if epsmac = 2.0d0 * epsmac end