C      ALGORITHM 813, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 27,NO. 3, September, 2001, P.  340--349.
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	Doc/
#	Doc/Makefile.Dp
#	Doc/Readme
#	Fortran77/
#	Fortran77/Dp/
#	Fortran77/Dp/Drivers/
#	Fortran77/Dp/Drivers/data2
#	Fortran77/Dp/Drivers/driver1.f
#	Fortran77/Dp/Drivers/driver2.f
#	Fortran77/Dp/Drivers/res1
#	Fortran77/Dp/Drivers/res2
#	Fortran77/Dp/Src/
#	Fortran77/Dp/Src/src.f
# This archive created: Thu Jan 31 16:39:31 2002
export PATH; PATH=/bin:$PATH
if test ! -d 'Doc'
then
	mkdir 'Doc'
fi
cd 'Doc'
if test -f 'Makefile.Dp'
then
	echo shar: will not over-write existing file "'Makefile.Dp'"
else
cat << "SHAR_EOF" > 'Makefile.Dp'
include makefile.inc
all: res1 res2
Libs1= $(PORT)
Objs1= driver1.o src.o
Libs2= $(PORT)
Objs2= driver2.o src.o

driver1: $(Objs1)
	$(F77LINK) $(F77LINKOPTS) -o driver1 $(Objs1)  $(Libs1)
res1: driver1
	driver1 >res1

driver2: $(Objs2)
	$(F77LINK) $(F77LINKOPTS) -o driver2 $(Objs2)  $(Libs2)
res2: driver2 data2
	driver2 <data2 >res2
SHAR_EOF
fi # end of overwriting check
if test -f 'Readme'
then
	echo shar: will not over-write existing file "'Readme'"
else
cat << "SHAR_EOF" > 'Readme'
   SPG: NONMONOTONE SPECTRAL PROJECTED GRADIENT METHOD
   ---------------------------------------------------

This file briefly describes Fortran 77 software that implements the
nonmonotone spectral projected gradient (SPG) algorithm. The SPG method
applies to constrained optimization problems of the form

                 min f(x) subject to x in Omega,

where Omega is a closed convex set in R^n. It is assumed that f is
defined and has continuous partial derivatives on an open set that
contains Omega. Users of the software must supply subroutines to compute
the function f(x), its gradient g(x), and projections of an arbitrary
point x onto Omega. Information about the Hessian matrix is not
required and the storage requirements are minimal. Hence the algorithm
is appropriate for large-scale convex-constrained optimization problems
with affordable projections onto the feasible set. Note that the
algorithm is also suitable for unconstrained optimization problems.
(In such cases, Omega = R^n and the projection routine does nothing.)

The algorithm is fully described in

E. G. Birgin, J. M. Martinez, and M. Raydan, "Nonmonotone spectral
projected gradient methods on convex sets", SIAM Journal on Optimization
10, pp. 1196-1211, 2000.

and

E. G. Birgin, J. M. Martinez, and M. Raydan, "SPG: software for
convex-constrained optimization", ACM Transactions on Mathematical
Software, 2001 (to appear).

It combines the projected gradient method with two new features in
optimization. First, it extends the typical globalization strategies
associated with these methods to the nonmonotone line search scheme
developed by Grippo, Lampariello and Lucidi. Second, it uses the
spectral steplength introduced by Barzilai and Borwein and analyzed by
Raydan. This choice of steplength requires little computational work and
greatly speeds up the convergence of gradient methods for unconstrained
problems.


USAGE OF THE SPG PACKAGE
------------------------

The distribution of the software has 3 source files called spg.f,
spgma1.f and spgma2.f. The first one contains the most important
subroutine of the package, which implements the SPG method itself. The
other two files are main programs to test the SPG method in examples 1
and 2 described below.


Example 1:
----------

This uses SPG to solve an easy (very easy) bound-constrained problem. The
aim of this example is not to show the capabilities of the method but to
be useful as an easy master file that could be modified to solve your
own problem.

The compilation process uses the Fortran compiler f77 or any other
compatible compiler. To compile this example on a Unix or Linux system,
type
   f77 spgma1.f spg.f -o spgma1
to generate an executable file called spgma1. Then type
   spgma1   (or  ./spgma1)
to see the results, which will appear on the screen.

To test SPG on your own problem you will need to modify subroutines
EvalF, EvalG and Proj in file spgma1.f. These subroutines evaluate the
objective function and its gradient and project an arbitrary point
onto the feasible region, respectively.


Example 2:
----------

This uses SPG to solve a location problem described in the ACM TOMS paper.
The master file of this example (spgma2.f) generates the problems,
calls the optimizer to solve them, and writes some reports (tables).

Brief description of the problem:

A regular grid is considered. This grid is the space at which a set
of cities (represented by polygons) will be built. Beforehand, a reserved
area (rectangle) is defined, where almost nothing can be done.
This reserved area will contain a hydro generation plant to supply energy
to the cities.  In the remaining space, the cities are built. At each point
of the grid (excluding the central region) a city (represented by a polygon)
will be built with a previously defined probability. To transmit energy from
the plant to the cities, a tower inside each city and a tower inside the
central region must be constructed. The objective of the problem is to
determine the location of these towers in order to minimize the sum of the
distances from each city tower to the central one.

To compile this example on a Unix or Linux system, type

  f77 spgma2.f spg.f -o spgma2

to generate an executable file called spgma2. Then type spgma2 and
follow the instructions on the screen. You will be asked for the number
of horizontal and vertical points in the grid, the probability of having
a polygon (city) at a grid point, and the minimum and maximum number of
vertices of the polygons. If you enter the values

   1  100  100  0.01  3  5

you will be generating and solving the first problem given in the paper.
These numbers mean that you are generating problem number 1, with 100
horizontal and 100 vertical points in the grid, that the probability of
having a polygon at a grid point is 0.01, and that the polygons will
have between 3 and 5 vertices or sides.

The data for the whole set of 45 problems is given below. Each column
means: identifier of the problem, number of horizontal and vertical
points in the grid, the probability of having a polygon (city) at a grid
point, and the minimum and maximum number of vertices of the polygons,
respectively. The characteristics of each of these problems (dimension
and number of constraints) as well as some information on the solutions
are given in the paper. For each problem, spgma2.f also generates an
output file with the optimal location of the city towers, the central
tower, and the distances from each city tower to the central one.

  1  100  100 0.01 3   5
  2  100  100 0.01 3  13
  3  100  100 0.01 3  21
  4  100  100 0.02 3   5
  5  100  100 0.02 3  13
  6  100  100 0.02 3  21
  7  100  100 0.03 3   5
  8  100  100 0.03 3  13
  9  100  100 0.03 3  21
 10  100  100 0.04 3   5
 11  100  100 0.04 3  13
 12  100  100 0.04 3  21
 13  100  100 0.05 3   5
 14  100  100 0.05 3  13
 15  100  100 0.05 3  21
 16  100 1000 0.01 3   5
 17  100 1000 0.01 3  13
 18  100 1000 0.01 3  21
 19  100 1000 0.02 3   5
 20  100 1000 0.02 3  13
 21  100 1000 0.02 3  21
 22  100 1000 0.03 3   5
 23  100 1000 0.03 3  13
 24  100 1000 0.03 3  21
 25  100 1000 0.04 3   5
 26  100 1000 0.04 3  13
 27  100 1000 0.04 3  21
 28  100 1000 0.05 3   5
 29  100 1000 0.05 3  13
 30  100 1000 0.05 3  21
 31 1000 1000 0.01 3   5
 32 1000 1000 0.01 3  13
 33 1000 1000 0.01 3  21
 34 1000 1000 0.02 3   5
 35 1000 1000 0.02 3  13
 36 1000 1000 0.02 3  21
 37 1000 1000 0.03 3   5
 38 1000 1000 0.03 3  13
 39 1000 1000 0.03 3  21
 40 1000 1000 0.04 3   5
 41 1000 1000 0.04 3  13
 42 1000 1000 0.04 3  21
 43 1000 1000 0.05 3   5
 44 1000 1000 0.05 3  13
 45 1000 1000 0.05 3  21

Observation:
------------

As the aim of this README file is not to reproduce the information in
the above publications, please read the numerical results section of the
papers for details. If, after reading it, you continue having problems
to run the programs, please do not hesitate in contact any of the authors.


AUTHORS INFORMATION
-------------------

Ernesto G. Birgin (egbirgin@ime.usp.br)
Department of Computer Science, 
Institute of Mathematics and Statistics, 
University of Sao Paulo,
Rua do Matao 1010, Cidade Universitaria, 
05508-900 Sao Paulo SP, Brazil

Jose Mario Martinez (martinez@ime.unicamp.br)
Department of Applied Mathematics, 
IMECC-UNICAMP, CP 6065, 
13081-970 Campinas SP, Brazil

Marcos Raydan (mraydan@reacciun.ve)
Departamento de Computacion, 
Facultad de Ciencias, 
Universidad Central de Venezuela, 
Ap. 47002, Caracas 1041-A,
Venezuela 

SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Fortran77'
then
	mkdir 'Fortran77'
fi
cd 'Fortran77'
if test ! -d 'Dp'
then
	mkdir 'Dp'
fi
cd 'Dp'
if test ! -d 'Drivers'
then
	mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'data2'
then
	echo shar: will not over-write existing file "'data2'"
else
cat << "SHAR_EOF" > 'data2'
  1  100  100 0.01 3   5
  2  100  100 0.01 3  13
  3  100  100 0.01 3  21
  4  100  100 0.02 3   5
  5  100  100 0.02 3  13
  6  100  100 0.02 3  21
  7  100  100 0.03 3   5
  8  100  100 0.03 3  13
  9  100  100 0.03 3  21
 10  100  100 0.04 3   5
 11  100  100 0.04 3  13
 12  100  100 0.04 3  21
 13  100  100 0.05 3   5
 14  100  100 0.05 3  13
 15  100  100 0.05 3  21
 16  100 1000 0.01 3   5
 17  100 1000 0.01 3  13
 18  100 1000 0.01 3  21
 19  100 1000 0.02 3   5
 20  100 1000 0.02 3  13
 21  100 1000 0.02 3  21
 22  100 1000 0.03 3   5
 23  100 1000 0.03 3  13
 24  100 1000 0.03 3  21
 25  100 1000 0.04 3   5
 26  100 1000 0.04 3  13
 27  100 1000 0.04 3  21
 28  100 1000 0.05 3   5
 29  100 1000 0.05 3  13
 30  100 1000 0.05 3  21
 31 1000 1000 0.01 3   5
 32 1000 1000 0.01 3  13
 33 1000 1000 0.01 3  21
 34 1000 1000 0.02 3   5
 35 1000 1000 0.02 3  13
 36 1000 1000 0.02 3  21
 37 1000 1000 0.03 3   5
 38 1000 1000 0.03 3  13
 39 1000 1000 0.03 3  21
 40 1000 1000 0.04 3   5
 41 1000 1000 0.04 3  13
 42 1000 1000 0.04 3  21
 43 1000 1000 0.05 3   5
 44 1000 1000 0.05 3  13
 45 1000 1000 0.05 3  21
SHAR_EOF
fi # end of overwriting check
if test -f 'driver1.f'
then
	echo shar: will not over-write existing file "'driver1.f'"
else
cat << "SHAR_EOF" > 'driver1.f'
      program spgma1

C     SPG Test driver. This master file uses SPG to solve
C     an easy (very easy) bound constrained problem.
C
C     This version 02 FEB 2001 by E.G.Birgin, J.M.Martinez and M.Raydan.
C     Final revision 03 JUL 2001 by E.G.Birgin, J.M.Martinez and M.Raydan.

C     PARAMETERS
      integer nmax
      parameter (nmax=100)

C     COMMON ARRAYS
      double precision l(nmax),u(nmax)

C     LOCAL SCALARS
      double precision pginfn,pgtwon,eps,eps2,f
      integer fcnt,flag,gcnt,i,iter,m,maxfc,maxit,n
      logical output

C     LOCAL ARRAYS
      double precision x(nmax)

C     EXTERNAL SUBROUTINES
      external spg

C     INTRINSIC FUNCTIONS
      intrinsic sqrt

C     COMMON BLOCKS
      common /bounds/l,u

C     DEFINE DIMENSION OF THE PROBLEM

      n = 10

C     DEFINE BOUNDS

      do i = 1,n
          l(i)= -100.0d0
          u(i)=   50.0d0
      end do

C     DEFINE INITIAL POINT

      do i = 1,n
          x(i) = 60.0d0
      end do

C     SET UP OPTIMIZER PARAMETERS 

      output= .true.
      maxit = 1000
      maxfc = 2000
      eps = 0.0d0
      eps2 = 1.0d-6
      m = 10

C     CALL THE OPTIMIZER

      call spg(n,x,m,eps,eps2,maxit,maxfc,output,f,pginfn,pgtwon,iter,
     +         fcnt,gcnt,flag)

C     WRITE STATISTICS

      write (*,fmt=9010) f,pginfn,sqrt(pgtwon),flag
      write (*,fmt=9020) iter,fcnt,gcnt

      stop

 9010 format (/' F= ',8X,1P,D17.10,/' PGINFNORM= ',1X,1P,D16.10,
     +       /' PGTWONORM= ',1X,1P,D16.10,/' FLAG= ',6X,I1)
 9020 format (/' ITER= ',1X,I10,/' FCNT= ',1X,I10,/' GCNT= ',1X,I10)

      end


      subroutine evalf(n,x,f,inform)

C     This subroutine computes the objective function.
C
C     On Entry:
C
C     n     integer,
C           size of the problem,
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
      integer i

      inform = 0

      f = 0.0d0
      do i = 1,n
          f = f + x(i)**2
      end do

      end


      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           size of the problem,
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
      integer i

      inform = 0

      do i = 1,n
          g(i) = 2.0d0 * x(i)
      end do

      end


      subroutine proj(n,x,inform)

C     This subroutine computes the projection of an arbitrary 
C     point onto the feasible set. Since the feasible set can be 
C     described in many ways, its information is in a common 
C     block.
C
C     On Entry:
C
C     n     integer,
C           size of the problem,
C
C     x     double precision x(n),
C           point that will be projected.
C
C     On Return
C
C     x     double precision x(n),
C           projected point,
C
C     inform integer,
C           termination parameter:
C           0= the projection was successfully done,
C           1= some error occurs in the projection.

C     PARAMETERS
      integer nmax
      parameter (nmax=100)

C     COMMON ARRAYS
      double precision l(nmax),u(nmax)

C     SCALAR ARGUMENTS
      integer n,inform

C     ARRAY ARGUMENTS
      double precision x(n)

C     LOCAL SCALARS
      integer i

C     COMMON BLOCKS
      common /bounds/l,u

      inform = 0

      do i = 1,n
          x(i) = max(l(i),min(x(i),u(i)))
      end do

      end
SHAR_EOF
fi # end of overwriting check
if test -f 'driver2.f'
then
	echo shar: will not over-write existing file "'driver2.f'"
else
cat << "SHAR_EOF" > 'driver2.f'
      program spgma2

C     SPG Test driver. This master file uses SPG to solve
C     the location problem described in Test Problems
C     section. It generates the problems, calls the optimizer
C     to solve them and writes the reports (tables).
C
C     This version 17 JAN 2000 by E.G.Birgin, J.M.Martinez and M.Raydan.
C     Reformatted 03 OCT 2000 by Tim Hopkins.
C     Final revision 03 JUL 2001 by E.G.Birgin, J.M.Martinez and M.Raydan.

C     PARAMETERS
      integer npmax
      parameter (npmax=50000)
      integer nvsmax
      parameter (nvsmax=30)
      integer nmax
      parameter (nmax=npmax*2)

C     COMMON SCALARS
      integer np,totnvs

C     COMMON ARRAYS
      double precision edges(npmax*nvsmax*3),vert(npmax*nvsmax*2)
      integer nvs(npmax)

C     LOCAL SCALARS
      double precision pginfn,pgtwon,eps,eps2,f,prob,rmax,rmin,xstep,
     +       ystep
      real time
      integer fcnt,flag,gcnt,i,iter,m,maxfc,maxit,n,nvsvma,nvsvmi,nx,ny,
     +        pnum
      logical output

C     LOCAL ARRAYS
      double precision x(nmax)
      real dum(2)

C     EXTERNAL SUBROUTINES
      external dtime,genpro,spg

C     INTRINSIC FUNCTIONS
      intrinsic sqrt

C     COMMON BLOCKS
      common /polyg/nvs,vert,edges,np,totnvs

C     DEFINE THE PROBLEM

1     continue

      write(*,fmt=9000)
      read (*,fmt=*,end=2) pnum,nx,ny,prob,nvsvmi,nvsvma

      if (nvsvmi.lt.3) then
          write (*,fmt=*) 'NVSVMI MUST BE GREATER THAN OR EQUAL TO 3'
          read (*,fmt=*)
      end if

      if (nvsvma.gt.nvsmax) then
          write (*,fmt=*) 'NVSVMA MUST BE LESS THAN OR EQUAL TO',nvsmax
          write (*,fmt=*) 'REDUCE NVSVMA OR INCREASE NVSMAX'
          read (*,fmt=*)
      end if

      if (nx*ny*prob.gt.npmax) then
          write (*,fmt=*) 'NX*NY*PROB MUST BE LESS THAN ',npmax
          write (*,fmt=*) 'REDUCE NX, NY OR PROB OR INCREASE NPMAX'
          read (*,fmt=*)
      end if

      xstep = 5.0d0
      ystep = 5.0d0

      rmin = 1.0d0
      rmax = 2.0d0

C     GENERATE THE PROBLEM

      call genpro(nx,ny,xstep,ystep,prob,nvsvmi,nvsvma,rmin,rmax)

      write (*,fmt=*) 'NUMBER OF POLYGONS: ',np
      write (*,fmt=*) 'NUMBER OF VERTICES (AND EDGES): ',totnvs

C     DEFINE INITIAL POINT

      n = 2*np
      do i = 1,n
          x(i) = 0.0d0
      end do

C     SET UP THE INPUT DATA OF THE OPTIMIZATION ALGORITHM

      output= .false.
      maxit = 1000
      maxfc = 2000
      eps = 0.0d0
      eps2 = 1.0d-06
      m = 10

C     CALL THE OPTIMIZER

      call dtime(dum)
      call spg(n,x,m,eps,eps2,maxit,maxfc,output,f,pginfn,pgtwon,iter,
     +         fcnt,gcnt,flag)
      call dtime(dum)

      time = dum(1) + dum(2)

C     WRITE STATISTICS

      write (*,fmt=9010) f,pginfn,sqrt(pgtwon),flag
      write (*,fmt=9020) iter,fcnt,gcnt
      write (*,fmt=9030) time

C     WRITE SOLUTION

      goto 1
2     continue

      stop

 9000 format (//'---------------------------------------'//)
 9010 format (/' F= ',8X,1P,D17.10,/' PGINFNORM= ',1X,1P,D16.10,
     +       /' PGTWONORM= ',1X,1P,D16.10,/' FLAG= ',6X,I1)
 9020 format (/' ITER= ',1X,I10,/' FCNT= ',1X,I10,/' GCNT= ',1X,I10)
 9030 format (/' TIME= ',F12.2,' SECONDS')
      end


      subroutine evalf(n,x,f,inform)

C     This subroutine computes the objective function.
C
C     On Entry:
C
C     n     integer,
C           size of the problem,
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 diff1,diff2,dist
      integer i,ndist

C     INTRINSIC FUNCTIONS
      intrinsic sqrt

      inform = 0

      f = 0.0d0
      ndist= n/2-1
      do i = 1,ndist
          diff1 = x(1) - x(2*i+1)
          diff2 = x(2) - x(2*i+2)
          dist = sqrt(diff1**2+diff2**2)
          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 / ndist

      end


      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           size of the problem,
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 diff1,diff2,dist
      integer i,ndist

C     INTRINSIC FUNCTIONS
      intrinsic sqrt

      inform = 0

      g(1) = 0.0d0
      g(2) = 0.0d0
      ndist= n/2-1
      do i = 1,ndist
          diff1 = x(1) - x(2*i+1)
          diff2 = x(2) - x(2*i+2)
          dist = sqrt(diff1**2+diff2**2)
          g(2*i+1) = - (diff1/dist) / ndist
          g(2*i+2) = - (diff2/dist) / ndist
          g(1) = g(1) - g(2*i+1)
          g(2) = g(2) - g(2*i+2)
      end do

      end


      subroutine proj(n,x,inform)

C     This subroutine computes the projection of an arbitrary 
C     point onto the feasible set. Since the feasible set can be 
C     described in many ways, its information is in a common 
C     block.
C
C     In this particular implementation of the projection subroutine
C     for the location problem, each point z_i in R^2 (stored at 
C     positions 2i-1 and 2i of x) must be projected onto polygon P_i.
C     See the Test Problems section for details.
C
C     On Entry:
C
C     n     integer,
C           size of the problem,
C
C     x     double precision x(n),
C           point that will be projected.
C
C     On Return
C
C     x     double precision x(n),
C           projected point,
C
C     inform integer,
C           termination parameter:
C           0= the projection was successfully done,
C           1= some error occurs in the projection.

C     PARAMETERS
      integer npmax
      parameter (npmax=50000)
      integer nvsmax
      parameter (nvsmax=30)

C     SCALAR ARGUMENTS
      integer n,inform

C     ARRAY ARGUMENTS
      double precision x(n)

C     COMMON SCALARS
      integer np,totnvs

C     COMMON ARRAYS 
      double precision edges(npmax*nvsmax*3),vert(npmax*nvsmax*2)
      integer nvs(npmax)

C     LOCAL SCALARS
      double precision xproj,yproj
      integer base,i

C     EXTERNAL SUBROUTINES
      external projp

C     COMMON BLOCKS
      common /polyg/nvs,vert,edges,np,totnvs

      inform = 0

      base = 0

      do i = 1,np

        ! PROJECT z_i ONTO P_i

          call projp(i,base,x(2*i-1),x(2*i),xproj,yproj,inform)

          x(2*i-1) = xproj
          x(2*i) = yproj

          base = base + nvs(i)

      end do

      end


      subroutine projp(p,base,x,y,xproj,yproj,inform)

C     This subroutine complements proj subroutine projecting
C     a point z_i in R^2 onto its corresponding polygon P_i.

C     On Entry:
C
C     p     integer,
C           index of the point z_p to be projected,
C
C     base  integer,
C           index which indicates the positions where the edges 
C           of the polygon  P_p are stored inside the "edges" 
C           vector which is part of the polygons description.
C
C     x     double precision,
C           first coordinate of z_i, i.e., z_i^1,
C
C     y     double precision,
C           second coordinate of z_i, i.e., z_i^2,
C
C     On Return
C
C     xproj double precision,
C           the projection of x,
C
C     yproj double precision,
C           the projection of y.
C
C     inform integer,
C           termination parameter:
C           0= the projection was successfully done,
C           1= some error occurs in the projection.

C     PARAMETERS
      integer npmax
      parameter (npmax=50000)
      integer nvsmax
      parameter (nvsmax=30)

C     SCALAR ARGUMENTS
      double precision x,xproj,y,yproj
      integer base,inform,p

C     COMMON SCALARS
      integer np,totnvs

C     COMMON ARRAYS
      double precision edges(npmax*nvsmax*3),vert(npmax*nvsmax*2)
      integer nvs(npmax)

C     LOCAL SCALARS
      double precision a,b,c,dist,mindis,nearvx,nearvy,px,py,v1x,v1y,
     +                 v2x,v2y,vx,vy
      integer i
      logical fact

C     LOCAL ARRAYS
      logical satisf(nvsmax)

C     INTRINSIC FUNCTIONS
      intrinsic max,min,sqrt

C     COMMON BLOCKS
      common /polyg/nvs,vert,edges,np,totnvs

C     TEST ALL EDGES OF THE POLYGON FOR SATISFIABILITY

      fact = .true.
      do i = base + 1,base + nvs(p)
          a = edges(3*i-2)
          b = edges(3*i-1)
          c = edges(3*i)
          if (a*x+b*y+c.le.0.0d0) then
              satisf(i-base) = .true.

          else
              fact = .false.
              satisf(i-base) = .false.
          end if
      end do

      if (fact) then
          xproj = x
          yproj = y
          return

      end if

C     SEARCH FOR THE CLOSEST VERTEX

      mindis = 1.0d+99
      do i = base + 1,base + nvs(p)
          vx = vert(2*i-1)
          vy = vert(2*i)
          dist = sqrt((vx-x)**2+ (vy-y)**2)
          if (dist.le.mindis) then
              mindis = dist
              nearvx = vx
              nearvy = vy
          end if

      end do

      xproj = nearvx
      yproj = nearvy

C     PROJECT ONTO THE VIOLATED CONSTRAINTS

      do i = base + 1,base + nvs(p)
          if (.not.satisf(i-base)) then
              a = edges(3*i-2)
              b = edges(3*i-1)
              c = edges(3*i)
              if (a.ne.0.0d0) then
                  py = (y-b* (c/a+x)/a)/ (b**2/a**2+1)
                  px = - (c+b*py)/a

              else
                  if (b.ne.0.0d0) then
                      px = (x-a* (c/b+y)/b)/ (a**2/b**2+1)
                      py = - (c+a*px)/b

                  else
                      write (*,fmt=*)
     +                  'ERROR IN PROBLEM DEFINITION (a=b=0 for a ',
     +                  'constraint of the type a x + b y + c = 0)'
                      inform= 1
                  end if

              end if

              v1x = vert(2*i-1)
              v1y = vert(2*i)
              if (i.ne.base+nvs(p)) then
                  v2x = vert(2*i+1)
                  v2y = vert(2*i+2)

              else
                  v2x = vert(2*base+1)
                  v2y = vert(2*base+2)
              end if

              if (min(v1x,v2x).le.px .and. px.le.max(v1x,v2x) .and.
     +            min(v1y,v2y).le.py .and. py.le.max(v1y,v2y)) then
                  dist = sqrt((px-x)**2+ (py-y)**2)
                  if (dist.le.mindis) then
                      mindis = dist
                      xproj = px
                      yproj = py
                  end if

              end if

          end if

      end do

      end


      subroutine genpro(nx,ny,xstep,ystep,prob,nvsvmi,nvsvma,rmin,rmax)

C     This subroutine generates a location problem (see the Figure).
C     First, a regular grid with nx horizontal points and ny vertical 
C     points in the positive orthant is considered. The points of the 
C     grid start at the origin with an horizontal distance of xstep and 
C     a vertical distance of ystept. This grid will be the space at which
C     the cities (represented by polygons) will be distributed. Before
C     building the cities, an area (rectangle) 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 nvsvmi, nvsvma, 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     nx    integer,
C           number of horizontal points in the grid,
C
C     ny    integer,
C           number of vertical points in the grid,
C
C     xstep double precision,
C           horizontal distance between points of the grid,
C
C     ystep double precision,
C           vertical distance between points of the grid,
C
C     prob  double precision, 
C           probability of defining a city at point of the grid
C           (0 <= prob <= 1),
C
C     nvsvmi integer,
C           parameter for the polygon generation (described
C           in genpol subroutine),
C
C     nvsvma integer,
C           parameter for the polygon generation (described
C           in genpol subroutine),
C
C     rmin  double precision,
C           parameter for the polygon generation (described
C           in genpol subroutine),
C
C     rmax  double precision,
C           parameter for the polygon generation (described
C           in genpol subroutine).
C
C     On output:
C
C     As described in the genpol subroutine, the output is
C     saved in the polyg common block.


C     PARAMETERS
      integer npmax
      parameter (npmax=50000)
      integer nvsmax
      parameter (nvsmax=30)

C     SCALAR ARGUMENTS
      double precision prob,rmax,rmin,xstep,ystep
      integer nvsvma,nvsvmi,nx,ny

C     COMMON SCALARS
      integer np,totnvs

C     COMMON ARRAYS
      double precision edges(npmax*nvsmax*3),vert(npmax*nvsmax*2)
      integer nvs(npmax)

C     LOCAL SCALARS
      double precision cx,cy,lx,ly,ux,uy
      integer i,j,seed

C     EXTERNAL FUNCTIONS
      real ran
      external ran

C     EXTERNAL SUBROUTINES
      external genpol

C     COMMON BLOCKS
      common /polyg/nvs,vert,edges,np,totnvs

C     SEED FOR THE RANDOM GENERATION

      seed = 760013

C     DEFINE BOX CONSTRAINTS FOR THE CENTRAL POINT

      lx = 0.40d0* (nx-1)*xstep
      ux = 0.60d0* (nx-1)*xstep
      ly = 0.40d0* (ny-1)*ystep
      uy = 0.60d0* (ny-1)*ystep

C     DEFINE CENTRAL POINT POLYGON

      nvs(1) = 4

      vert(1) = lx
      vert(2) = ly
      vert(3) = lx
      vert(4) = uy
      vert(5) = ux
      vert(6) = uy
      vert(7) = ux
      vert(8) = ly

      edges(1) = -1.0d0
      edges(2) = 0.0d0
      edges(3) = lx
      edges(4) = 0.0d0
      edges(5) = 1.0d0
      edges(6) = -uy
      edges(7) = 1.0d0
      edges(8) = 0.0d0
      edges(9) = -ux
      edges(10) = 0.0d0
      edges(11) = -1.0d0
      edges(12) = ly

C     DEFINE CITY-POLYGONS CENTERED AT THE GRID POINTS AND OUTSIDE
C     THE CENTRAL REGION

      np = 1
      totnvs = 4

      do i = 0,nx - 1
          do j = 0,ny - 1

              cx = i*xstep
              cy = j*ystep

              if ((cx.lt.lx-xstep.or.cx.gt.ux+xstep.or.
     +            cy.lt.ly-ystep.or.cy.gt.uy+ystep) .and.
     +            ran(seed).le.prob) then

                ! GENERATE A NEW POLYGON

                  np = np + 1
                  call genpol(cx,cy,nvsvmi,nvsvma,rmin,rmax,seed)
                  totnvs = totnvs + nvs(np)

              end if

          end do
      end do

      end


      subroutine genpol(cx,cy,nvsvmi,nvsvma,rmin,rmax,seed)

C     This subroutine generates a polygon in R^2 with its 
C     vertices in a sphere centered at point (cx,cy). The
C     number of vertices is randomly generated 
C     satisfying nvsvmi <= number of vertices  <= 
C     nvsvma. The ratio of the sphere is also randomly 
C     generated satisfying rmin <= ratio < rmax. The 
C     generated polygon is stored in the common block "polyg".
C 
C     On Entry:
C
C     cx    double precision,
C           first coordinate of the center of the sphere,
C
C     cy    double precision,
C           second coordinate of the center of the sphere,
C
C     nvsvmi integer,
C           minimum number of vertices,
C
C     nvsvma integer,
C           maximum number of vertices,
C
C     rmin double,
C           minimum ratio of the sphere,
C
C     rmax double,
C           maximum ratio of the sphere,
C
C     seed double,
C           seed for the random generation.
C
C     On Output:
C     
C     The generated polygon is stored in the polyg common block
C     described below.
C
C     Common block polyg:
C
C     common /polyg/nvs,vert,edges,np,totnvs
C 
C     This structure represents, at any time, np polygons.
C     Position i of array nvs indicates the number of vertices
C     of polygon i. Arrays vert and edges store the 
C     vertices and edges of the polygons.
C
C     For example, if nvs(1) = 3 it indicates that the first
C     polygon has 3 vertices (edges). Then, if the vertices
C     are (x1,y1), (x2,y2) and (x3,y3), we have that vert(1) 
C     = x1, vert(2) = y1, vert(3) = x2, vert(4) = y2, vert(5) 
C     = x3, and vert(6) = y3. And, if the edges (written as
C     ax + by + c = 0) are a1 x + b1 y + c1 = 0, a2 x + b2 y 
C     + c2 = 0, and a3 x + b3 y + c3 = 0 then edges(1) = a1,
C     edges(2) = b1, edges(3) = c1, edges(4) = a2, edges(5) = 
C     b2, edges(6) = c2, edges(7) = a3, edges(8) = b3 and
C     edges(9) = c3.
C
C     totnvs indicates the total number of vertices  
C     of the set of polygons. This information is used when
C     a new polygon is created to know the first free position
C     of arrays vert and edges at which the vertices and edges 
C     of the new polygon will be saved.
C
C     Two additional details: 
C
C     1) For each polygon, the vertices are ordered clockwise 
C     and edge i corresponds to the edge between vertices 
C     i and i+1 (0 if i=n).
C
C     2) For each edge of the form ax + bx + c = 0, constants
C     a, b and c are chosen in such a way that 
C     (|a| = 1 or |b| = 1) and (a cx + b cy + c <= 0).

C     PARAMETERS
      integer npmax
      parameter (npmax=50000)
      integer nvsmax
      parameter (nvsmax=30)
      double precision pi
      parameter (pi=3.1415926535898d0)

C     SCALAR ARGUMENTS
      double precision cx,cy,rmax,rmin
      integer nvsvma,nvsvmi,seed

C     COMMON SCALARS
      integer np,totnvs

C     COMMON  ARRAYS
      double precision edges(npmax*nvsmax*3),vert(npmax*nvsmax*2)
      integer nvs(npmax)

C     LOCAL SCALARS
      double precision r
      integer i

C     LOCAL ARRAYS
      double precision angl(nvsmax)

C     EXTERNAL FUNCTIONS
      real ran
      external ran

C     EXTERNAL SUBROUTINES
      external class,constr

C     INTRINSIC FUNCTIONS
      intrinsic dcos,dsin,int

C     COMMON BLOCKS
      common /polyg/nvs,vert,edges,np,totnvs

C     GENERATE THE NUMBER OF VERTICES 

      nvs(np) = nvsvmi + int((nvsvma-nvsvmi+1)*ran(seed))

C     GENERATE THE RATIO OF THE SPHERE

      r = rmin + (rmax-rmin)*ran(seed)

C     GENERATE ALL ANGLES SATISFYING 0 <= ANGLE_I < 2*PI

      do i = 1,nvs(np)
          angl(i) = 2*pi*ran(seed)
      end do

C     CLASSIFY THE ANGLES IN DECREASING ORDER

      call class(nvs(np),angl)

C     CONSTRUCT THE VERTICES

      do i = 1,nvs(np)
          vert(2* (totnvs+i)-1) = cx + r*dcos(angl(i))
          vert(2* (totnvs+i)) = cy + r*dsin(angl(i))
      end do

C     CONSTRUCT THE EDGES

      do i = totnvs + 1,totnvs + nvs(np) - 2
          call constr(vert(2*i-1),vert(2*i),vert(2*i+1),vert(2*i+2),
     +                vert(2*i+3),vert(2*i+4),edges(3*i-2),edges(3*i-1),
     +                edges(3*i))
      end do

      i = totnvs + nvs(np) - 1

      call constr(vert(2*i-1),vert(2*i),vert(2*i+1),vert(2*i+2),
     +            vert(2*totnvs+1),vert(2*totnvs+2),edges(3*i-2),
     +            edges(3*i-1),edges(3*i))

      i = totnvs + nvs(np)

      call constr(vert(2*i-1),vert(2*i),vert(2*totnvs+1),
     +            vert(2*totnvs+2),vert(2*totnvs+3),vert(2*totnvs+4),
     +            edges(3*i-2),edges(3*i-1),edges(3*i))

      end


      subroutine constr(x1,y1,x2,y2,x3,y3,a,b,c)

C     This subroutine computes the real constants a, b and c of 
C     the straight line ax + by + c = 0 in R^2 defined by the 
C     points (x1,y1) and (x2,y2); such that the point (x3,y3) 
C     satisfies the constraint ax + by + c <= 0.
C
C     On Entry:
C
C     x1    double precision,
C           first coordinate of point (x1,y1),
C
C     y1    double precision,
C           second coordinate of point (x1,y1),
C
C     x2    double precision,
C           first coordinate of point (x2,y2),
C
C     y2    double precision,
C           second coordinate of point (x2,y2),
C
C     x3    double precision,
C           first coordinate of point (x3,y3),
C
C     y3    double precision,
C           second coordinate of point (x3,y3).
C
C     On Return
C
C     a,b,c double precision
C           the desired constants.

C     SCALAR ARGUMENTS
      double precision a,b,c,x1,x2,x3,y1,y2,y3

      if (x1.eq.x2 .and. y1.eq.y2) then
          write (*,fmt=*)
     +      'ERROR IN FUNCTION CONSTRAINT: X1=X2 AND Y1=Y2'
      end if

      if (y1.ne.y2) then
          a = 1.0d0
          b = - (x2-x1)/ (y2-y1)
          c = - (x1+b*y1)

      else
          a = 0.0d0
          b = 1.0d0
          c = -y1
      end if

      if (a*x3+b*y3+c.gt.0.0d0) then
          a = -a
          b = -b
          c = -c
      end if

      end


      subroutine class(n,x)

C     This subroutine classifies the elements of a vector in
C     decreasing order, i.e., on output: x(1) >= x(2) >= 
C     ... >= x(n).
C
C     On Entry:
C
C     n     integer,
C           number of elements of the vector to be classified,
C
C     x     double precision x(n),
C           vector to be classified.
C
C     On Return
C
C     x     double precision x(n),
C           classified vector.

C     SCALAR ARGUMENTS
      integer n

C     ARRAY ARGUMENTS
      double precision x(n)

C     LOCAL SCALARS
      double precision aux,xmax
      integer i,j,pos

      do i = 1,n

          xmax = x(i)
          pos = i

          do j = i + 1,n
              if (x(j).gt.xmax) then
                  xmax = x(j)
                  pos = j
              end if
          end do

          if (pos.ne.i) then
              aux = x(i)
              x(i) = x(pos)
              x(pos) = aux
          end if

      end do

      end













SHAR_EOF
fi # end of overwriting check
if test -f 'res1'
then
	echo shar: will not over-write existing file "'res1'"
else
cat << "SHAR_EOF" > 'res1'
ITER=          0 F=  2.5000000000D+04 PGINFNORM=  1.0000000000D+02
ITER=          1 F=  2.4010000000D+04 PGINFNORM=  9.8000000000D+01
ITER=          2 F=  0.0000000000D+00 PGINFNORM=  0.0000000000D+00

 F=          0.0000000000D+00
 PGINFNORM=  0.0000000000D+00
 PGTWONORM=  0.0000000000D+00
 FLAG=       0

 ITER=           2
 FCNT=           3
 GCNT=           3
SHAR_EOF
fi # end of overwriting check
if test -f 'res2'
then
	echo shar: will not over-write existing file "'res2'"
else
cat << "SHAR_EOF" > 'res2'


---------------------------------------


 NUMBER OF POLYGONS:           86
 NUMBER OF VERTICES (AND EDGES):          343

 F=          1.8532259472D+02
 PGINFNORM=  1.3099906937D-08
 PGTWONORM=  1.3177603722D-08
 FLAG=       1

 ITER=          14
 FCNT=          15
 GCNT=          15

 TIME=         0.12 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:           86
 NUMBER OF VERTICES (AND EDGES):          700

 F=          1.9592273928D+02
 PGINFNORM=  4.5548063099D-09
 PGTWONORM=  4.6044528622D-09
 FLAG=       1

 ITER=          14
 FCNT=          15
 GCNT=          15

 TIME=         0.04 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:           85
 NUMBER OF VERTICES (AND EDGES):         1018

 F=          1.9144177470D+02
 PGINFNORM=  4.8637156169D-07
 PGTWONORM=  5.2286896014D-07
 FLAG=       1

 ITER=          14
 FCNT=          15
 GCNT=          15

 TIME=         0.06 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          174
 NUMBER OF VERTICES (AND EDGES):          701

 F=          1.9445764952D+02
 PGINFNORM=  3.8522756540D-07
 PGTWONORM=  4.2872333583D-07
 FLAG=       1

 ITER=          24
 FCNT=          27
 GCNT=          25

 TIME=         0.11 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          181
 NUMBER OF VERTICES (AND EDGES):         1450

 F=          2.0213902708D+02
 PGINFNORM=  4.1171244902D-07
 PGTWONORM=  6.4752916313D-07
 FLAG=       1

 ITER=          23
 FCNT=          26
 GCNT=          24

 TIME=         0.15 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          176
 NUMBER OF VERTICES (AND EDGES):         2104

 F=          1.9484973167D+02
 PGINFNORM=  2.9913155686D-07
 PGTWONORM=  5.2346409717D-07
 FLAG=       1

 ITER=          20
 FCNT=          21
 GCNT=          21

 TIME=         0.15 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          257
 NUMBER OF VERTICES (AND EDGES):         1033

 F=          1.9523279096D+02
 PGINFNORM=  7.3946395673D-07
 PGTWONORM=  9.5590499054D-07
 FLAG=       1

 ITER=          19
 FCNT=          20
 GCNT=          20

 TIME=         0.13 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          256
 NUMBER OF VERTICES (AND EDGES):         2103

 F=          1.9893980522D+02
 PGINFNORM=  3.0896751468D-07
 PGTWONORM=  4.2273280213D-07
 FLAG=       1

 ITER=          24
 FCNT=          25
 GCNT=          25

 TIME=         0.22 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          262
 NUMBER OF VERTICES (AND EDGES):         3216

 F=          1.9482860880D+02
 PGINFNORM=  4.6882860261D-07
 PGTWONORM=  4.7237807379D-07
 FLAG=       1

 ITER=          34
 FCNT=          36
 GCNT=          35

 TIME=         0.41 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          359
 NUMBER OF VERTICES (AND EDGES):         1437

 F=          1.9688544937D+02
 PGINFNORM=  3.4909385249D-07
 PGTWONORM=  3.6310805447D-07
 FLAG=       1

 ITER=          19
 FCNT=          20
 GCNT=          20

 TIME=         0.18 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          349
 NUMBER OF VERTICES (AND EDGES):         2854

 F=          1.9687491958D+02
 PGINFNORM=  4.7915796131D-07
 PGTWONORM=  5.3882284221D-07
 FLAG=       1

 ITER=          29
 FCNT=          30
 GCNT=          30

 TIME=         0.36 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          349
 NUMBER OF VERTICES (AND EDGES):         4328

 F=          1.9433238446D+02
 PGINFNORM=  4.8600622904D-07
 PGTWONORM=  6.7702589275D-07
 FLAG=       1

 ITER=          23
 FCNT=          24
 GCNT=          24

 TIME=         0.36 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          435
 NUMBER OF VERTICES (AND EDGES):         1743

 F=          1.9681184123D+02
 PGINFNORM=  6.1534200313D-07
 PGTWONORM=  9.1061741874D-07
 FLAG=       1

 ITER=          14
 FCNT=          15
 GCNT=          15

 TIME=         0.16 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          432
 NUMBER OF VERTICES (AND EDGES):         3522

 F=          1.9828448531D+02
 PGINFNORM=  5.6335181853D-07
 PGTWONORM=  6.0794455611D-07
 FLAG=       1

 ITER=          25
 FCNT=          26
 GCNT=          26

 TIME=         0.38 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          430
 NUMBER OF VERTICES (AND EDGES):         5257

 F=          1.9739060109D+02
 PGINFNORM=  6.8659574026D-07
 PGTWONORM=  8.4739347247D-07
 FLAG=       1

 ITER=          26
 FCNT=          29
 GCNT=          27

 TIME=         0.49 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          935
 NUMBER OF VERTICES (AND EDGES):         3729

 F=          1.2858865478D+03
 PGINFNORM=  1.3228140006D-07
 PGTWONORM=  1.3264041724D-07
 FLAG=       1

 ITER=          20
 FCNT=          21
 GCNT=          21

 TIME=         0.49 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          928
 NUMBER OF VERTICES (AND EDGES):         7399

 F=          1.3242945994D+03
 PGINFNORM=  5.9177546063D-07
 PGTWONORM=  6.6201299375D-07
 FLAG=       1

 ITER=          30
 FCNT=          33
 GCNT=          31

 TIME=         0.98 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:          940
 NUMBER OF VERTICES (AND EDGES):        11200

 F=          1.3041918432D+03
 PGINFNORM=  8.7161379270D-07
 PGTWONORM=  8.8617729654D-07
 FLAG=       1

 ITER=          28
 FCNT=          30
 GCNT=          29

 TIME=         1.15 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         1841
 NUMBER OF VERTICES (AND EDGES):         7362

 F=          1.3108100872D+03
 PGINFNORM=  3.9998940338D-07
 PGTWONORM=  4.5397686033D-07
 FLAG=       1

 ITER=          66
 FCNT=          82
 GCNT=          67

 TIME=         3.18 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         1825
 NUMBER OF VERTICES (AND EDGES):        14574

 F=          1.3147211801D+03
 PGINFNORM=  4.1348357627D-07
 PGTWONORM=  4.1631782967D-07
 FLAG=       1

 ITER=          52
 FCNT=          64
 GCNT=          53

 TIME=         3.30 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         1841
 NUMBER OF VERTICES (AND EDGES):        22049

 F=          1.2927581924D+03
 PGINFNORM=  6.7566725193D-07
 PGTWONORM=  8.9656623602D-07
 FLAG=       1

 ITER=          42
 FCNT=          47
 GCNT=          43

 TIME=         3.31 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         2766
 NUMBER OF VERTICES (AND EDGES):        11102

 F=          1.3076190662D+03
 PGINFNORM=  3.0531191442D-07
 PGTWONORM=  4.9931310285D-07
 FLAG=       1

 ITER=          33
 FCNT=          35
 GCNT=          34

 TIME=         2.38 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         2795
 NUMBER OF VERTICES (AND EDGES):        22511

 F=          1.2663399290D+03
 PGINFNORM=  4.3137015382D-07
 PGTWONORM=  4.8857486504D-07
 FLAG=       1

 ITER=          55
 FCNT=          68
 GCNT=          56

 TIME=         5.37 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         2831
 NUMBER OF VERTICES (AND EDGES):        34445

 F=          1.3122736233D+03
 PGINFNORM=  7.8603625298D-07
 PGTWONORM=  9.6854750508D-07
 FLAG=       1

 ITER=          86
 FCNT=         113
 GCNT=          87

 TIME=        10.66 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         3779
 NUMBER OF VERTICES (AND EDGES):        15104

 F=          1.3110867348D+03
 PGINFNORM=  2.5001872928D-07
 PGTWONORM=  4.1614136310D-07
 FLAG=       1

 ITER=          65
 FCNT=          92
 GCNT=          66

 TIME=         6.45 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         3815
 NUMBER OF VERTICES (AND EDGES):        30699

 F=          1.3009282790D+03
 PGINFNORM=  4.2950250645D-07
 PGTWONORM=  4.9035339224D-07
 FLAG=       1

 ITER=          50
 FCNT=          60
 GCNT=          51

 TIME=         6.59 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         3852
 NUMBER OF VERTICES (AND EDGES):        46704

 F=          1.3214063723D+03
 PGINFNORM=  5.7753186411D-07
 PGTWONORM=  6.7699253721D-07
 FLAG=       1

 ITER=          37
 FCNT=          42
 GCNT=          38

 TIME=         6.15 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         4735
 NUMBER OF VERTICES (AND EDGES):        18937

 F=          1.3049304651D+03
 PGINFNORM=  3.8924918044D-07
 PGTWONORM=  8.0643760240D-07
 FLAG=       1

 ITER=          85
 FCNT=         122
 GCNT=          86

 TIME=        10.60 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         4767
 NUMBER OF VERTICES (AND EDGES):        38247

 F=          1.3056273706D+03
 PGINFNORM=  2.5131157599D-07
 PGTWONORM=  2.7780731830D-07
 FLAG=       1

 ITER=          50
 FCNT=          60
 GCNT=          51

 TIME=         8.33 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         4836
 NUMBER OF VERTICES (AND EDGES):        58622

 F=          1.3087534158D+03
 PGINFNORM=  4.4804482968D-07
 PGTWONORM=  6.9564057350D-07
 FLAG=       1

 ITER=          44
 FCNT=          47
 GCNT=          45

 TIME=         9.33 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         9680
 NUMBER OF VERTICES (AND EDGES):        38790

 F=          1.9678500970D+03
 PGINFNORM=  4.9539175961D-07
 PGTWONORM=  6.5738040742D-07
 FLAG=       1

 ITER=          18
 FCNT=          20
 GCNT=          19

 TIME=         4.84 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         9644
 NUMBER OF VERTICES (AND EDGES):        77413

 F=          1.9658121929D+03
 PGINFNORM=  2.4842211133D-07
 PGTWONORM=  3.9109092003D-07
 FLAG=       1

 ITER=          24
 FCNT=          30
 GCNT=          25

 TIME=         8.63 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:         9639
 NUMBER OF VERTICES (AND EDGES):       116195

 F=          1.9594484247D+03
 PGINFNORM=  2.4242581276D-07
 PGTWONORM=  4.7378258335D-07
 FLAG=       1

 ITER=          20
 FCNT=          21
 GCNT=          21

 TIME=         8.98 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        19117
 NUMBER OF VERTICES (AND EDGES):        76550

 F=          1.9683454765D+03
 PGINFNORM=  4.2561168812D-08
 PGTWONORM=  1.0048593268D-07
 FLAG=       1

 ITER=          14
 FCNT=          15
 GCNT=          15

 TIME=         7.54 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        19093
 NUMBER OF VERTICES (AND EDGES):       153156

 F=          1.9640035411D+03
 PGINFNORM=  5.5163764046D-07
 PGTWONORM=  8.0557421950D-07
 FLAG=       1

 ITER=          15
 FCNT=          17
 GCNT=          16

 TIME=        10.71 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        19114
 NUMBER OF VERTICES (AND EDGES):       230190

 F=          1.9728371669D+03
 PGINFNORM=  3.3214200812D-07
 PGTWONORM=  7.1934387440D-07
 FLAG=       1

 ITER=          19
 FCNT=          21
 GCNT=          20

 TIME=        16.79 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        28770
 NUMBER OF VERTICES (AND EDGES):       115301

 F=          1.9730821434D+03
 PGINFNORM=  4.6466084314D-08
 PGTWONORM=  1.7614356564D-07
 FLAG=       1

 ITER=          14
 FCNT=          15
 GCNT=          15

 TIME=        11.49 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        28799
 NUMBER OF VERTICES (AND EDGES):       230878

 F=          1.9738282112D+03
 PGINFNORM=  1.4175066099D-07
 PGTWONORM=  5.6768062768D-07
 FLAG=       1

 ITER=          18
 FCNT=          20
 GCNT=          19

 TIME=        19.14 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        28767
 NUMBER OF VERTICES (AND EDGES):       346048

 F=          1.9735553698D+03
 PGINFNORM=  1.7218553694D-07
 PGTWONORM=  5.7205827017D-07
 FLAG=       1

 ITER=          13
 FCNT=          14
 GCNT=          14

 TIME=        17.36 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        38409
 NUMBER OF VERTICES (AND EDGES):       153784

 F=          1.9747875258D+03
 PGINFNORM=  2.5164081308D-07
 PGTWONORM=  9.3978722828D-07
 FLAG=       1

 ITER=          14
 FCNT=          15
 GCNT=          15

 TIME=        15.38 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        38568
 NUMBER OF VERTICES (AND EDGES):       309252

 F=          1.9797932200D+03
 PGINFNORM=  1.9212620828D-07
 PGTWONORM=  6.6858179959D-07
 FLAG=       1

 ITER=          25
 FCNT=          28
 GCNT=          26

 TIME=        35.51 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        38506
 NUMBER OF VERTICES (AND EDGES):       463718

 F=          1.9798134073D+03
 PGINFNORM=  2.7334908737D-07
 PGTWONORM=  7.8077288684D-07
 FLAG=       1

 ITER=          24
 FCNT=          26
 GCNT=          25

 TIME=        42.17 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        48004
 NUMBER OF VERTICES (AND EDGES):       192152

 F=          1.9731285861D+03
 PGINFNORM=  9.2753907666D-08
 PGTWONORM=  2.9721472708D-07
 FLAG=       1

 ITER=          18
 FCNT=          20
 GCNT=          19

 TIME=        24.33 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        48209
 NUMBER OF VERTICES (AND EDGES):       386121

 F=          1.9775617204D+03
 PGINFNORM=  5.9787862483D-07
 PGTWONORM=  7.5367590050D-07
 FLAG=       1

 ITER=          13
 FCNT=          14
 GCNT=          14

 TIME=        23.52 SECONDS


---------------------------------------


 NUMBER OF POLYGONS:        48126
 NUMBER OF VERTICES (AND EDGES):       578648

 F=          1.9802088407D+03
 PGINFNORM=  1.6024887373D-07
 PGTWONORM=  7.4644171077D-07
 FLAG=       1

 ITER=          17
 FCNT=          19
 GCNT=          18

 TIME=        36.68 SECONDS


---------------------------------------


SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
	mkdir 'Src'
fi
cd 'Src'
if test -f 'src.f'
then
	echo shar: will not over-write existing file "'src.f'"
else
cat << "SHAR_EOF" > 'src.f'
      subroutine spg(n,x,m,eps,eps2,maxit,maxfc,output,f,pginfn,pgtwon,
     +               iter,fcnt,gcnt,flag)

C     Subroutine SPG implements the Spectral Projected Gradient
C     Method (Version 2: "continuous projected gradient direction") 
C     to find the local minimizers of a given function with convex
C     constraints, described in 
C
C     E. G. Birgin, J. M. Martinez, and M. Raydan, "Nonmonotone 
C     spectral projected gradient methods on convex sets", SIAM 
C     Journal on Optimization 10, pp. 1196-1211, 2000.
C 
C     and
C 
C     E. G. Birgin, J. M. Martinez, and M. Raydan, "SPG: software 
C     for convex-constrained optimization", ACM Transactions on 
C     Mathematical Software, 2001 (to appear).
C
C     The user must supply the external subroutines evalf, evalg 
C     and proj to evaluate the objective function and its gradient 
C     and to project an arbitrary point onto the feasible region.
C
C     This version 17 JAN 2000 by E.G.Birgin, J.M.Martinez and M.Raydan.
C     Reformatted 03 OCT 2000 by Tim Hopkins.
C     Final revision 03 JUL 2001 by E.G.Birgin, J.M.Martinez and M.Raydan.
C
C     On Entry:
C
C     n     integer,
C           size of the problem,
C
C     x     double precision x(n),
C           initial guess,
C
C     m     integer,
C           number of previous function values to be considered 
C           in the nonmonotone line search,
C
C     eps   double precision,
C           stopping criterion: ||projected grad||_inf < eps,
C
C     eps2  double precision,
C           stopping criterion: ||projected grad||_2 < eps2,
C
C     maxit integer,
C           maximum number of iterations,
C
C     maxfc integer,
C           maximum number of function evaluations,
C
C     output logical,
C           true: print some information at each iteration,
C           false: no print.
C
C     On Return:
C
C     x     double precision x(n),
C           approximation to the local minimizer,
C
C     f     double precision,
C           function value at the approximation to the local
C           minimizer,
C
C     pginfn double precision,
C           ||projected grad||_inf at the final iteration,
C
C     pgtwon double precision,
C           ||projected grad||_2^2 at the final iteration,
C
C     iter  integer,
C           number of iterations,
C
C     fcnt  integer,
C           number of function evaluations,
C
C     gcnt  integer,
C           number of gradient evaluations,
C
C     flag  integer,
C           termination parameter:
C           0= convergence with projected gradient infinite-norm,
C           1= convergence with projected gradient 2-norm,
C           2= too many iterations,
C           3= too many function evaluations,
C           4= error in proj subroutine,
C           5= error in evalf subroutine,
C           6= error in evalg subroutine.

C     PARAMETERS
      double precision lmin
      parameter (lmin=1.0d-30)
      double precision lmax
      parameter (lmax=1.0d+30)
      integer nmax
      parameter (nmax=100000)
      integer mmax
      parameter (mmax=100)

C     SCALAR ARGUMENTS
      double precision pginfn,pgtwon,eps,eps2,f
      integer fcnt,flag,gcnt,iter,m,maxfc,maxit,n
      logical output

C     ARRAY ARGUMENTS
      double precision x(n)

C     LOCAL SCALARS
      double precision fbest,fnew,gtd,lambda,sts,sty
      integer i,lsflag,inform

C     LOCAL 
      double precision pg(nmax),g(nmax),gnew(nmax),s(nmax),y(nmax),
     +                 d(nmax),xbest(nmax),xnew(nmax),lastfv(0:mmax-1)

C     EXTERNAL SUBROUTINES
      external ls,evalf,evalg,proj

C     INTRINSIC FUNCTIONS
      intrinsic abs,max,min,mod

C     INITIALIZATION

      iter = 0
      fcnt = 0
      gcnt = 0

      do i = 0,m - 1
          lastfv(i) = -1.0d+99
      end do

C     PROJECT INITIAL GUESS

      call proj(n,x,inform)

      if (inform .ne. 0) then
        ! ERROR IN PROJ SUBROUTINE, STOP
          flag = 4
          go to 200
      end if

C     INITIALIZE BEST SOLUTION

      do i = 1,n
          xbest(i) = x(i)
      end do

C     EVALUATE FUNCTION AND GRADIENT

      call evalf(n,x,f,inform)
      fcnt = fcnt + 1

      if (inform .ne. 0) then
        ! ERROR IN EVALF SUBROUTINE, STOP
          flag = 5
          go to 200
      end if

      call evalg(n,x,g,inform)
      gcnt = gcnt + 1

      if (inform .ne. 0) then
        ! ERROR IN EVALG SUBROUTINE, STOP
          flag = 6
          go to 200
      end if

C     STORE FUNCTION VALUE FOR THE NONMONOTONE LINE SEARCH

      lastfv(0) = f

C     INITIALIZE BEST FUNCTION VALUE

      fbest = f

C     COMPUTE CONTINUOUS PROJECTED GRADIENT (AND ITS NORMS)

      do i = 1,n
          pg(i) = x(i) - g(i)
      end do

      call proj(n,pg,inform)

      if (inform .ne. 0) then
        ! ERROR IN PROJ SUBROUTINE, STOP
          flag = 4
          go to 200
      end if

      pgtwon = 0.0D0
      pginfn = 0.0D0
      do i = 1,n
          pg(i) = pg(i) - x(i)
          pgtwon = pgtwon + pg(i)**2
          pginfn = max(pginfn,abs(pg(i)))
      end do

C     PRINT ITERATION INFORMATION

      if (output) then
          write (*,fmt=9010) iter,f,pginfn
      end if

C     DEFINE INITIAL SPECTRAL STEPLENGTH

      if (pginfn .ne. 0.0d0) then
          lambda =  min(lmax,max(lmin,1.0d0/pginfn))
      end if

C     MAIN LOOP

C     TEST STOPPING CRITERIA

 100  continue

      if (pginfn .le. eps) then
        ! GRADIENT INFINITE-NORM STOPPING CRITERION SATISFIED, STOP
          flag = 0
          go to 200
      end if

      if (pgtwon .le. eps2**2) then
        ! GRADIENT 2-NORM STOPPING CRITERION SATISFIED, STOP
          flag = 1
          go to 200
      end if

      if (iter .gt. maxit) then
        ! MAXIMUM NUMBER OF ITERATIONS EXCEEDED, STOP
          flag = 2
          go to 200
      end if

      if (fcnt .gt. maxfc) then
        ! MAXIMUM NUMBER OF FUNCTION EVALUATIONS EXCEEDED, STOP
          flag = 3
          go to 200
      end if

C     DO AN ITERATION

      iter = iter + 1

C     COMPUTE THE SPECTRAL PROJECTED GRADIENT DIRECTION
C     AND <G,D>

      do i = 1,n
          d(i) = x(i) - lambda*g(i)
      end do

      call proj(n,d,inform)

      if (inform .ne. 0) then
        ! ERROR IN PROJ SUBROUTINE, STOP
          flag = 4
          go to 200
      end if

      gtd = 0.0d0
      do i = 1,n
          d(i) = d(i) - x(i)
          gtd = gtd + g(i)*d(i)
      end do

C     NONMONOTONE LINE SEARCH

      call ls(n,x,f,d,gtd,m,lastfv,maxfc,fcnt,fnew,xnew,lsflag)

      if (lsflag .eq. 3) then
        ! THE NUMBER OF FUNCTION EVALUATIONS WAS EXCEEDED 
        ! INSIDE  THE LINE SEARCH, STOP
          flag = 3
          go to 200
      end if

C     SET NEW FUNCTION VALUE AND SAVE IT FOR THE NONMONOTONE
C     LINE SEARCH

      f = fnew
      lastfv(mod(iter,m)) = f

C     COMPARE THE NEW FUNCTION VALUE AGAINST THE BEST FUNCTION 
C     VALUE  AND, IF SMALLER, UPDATE THE BEST FUNCTION 
C     VALUE AND THE CORRESPONDING BEST POINT

      if (f .lt. fbest) then
          fbest = f
          do i = 1,n
              xbest(i) = xnew(i)
          end do
      end if

C     EVALUATE THE GRADIENT AT THE NEW ITERATE

      call evalg(n,xnew,gnew,inform)
      gcnt = gcnt + 1

      if (inform .ne. 0) then
        ! ERROR IN EVALG SUBROUTINE, STOP
          flag = 6
          go to 200
      end if

C     COMPUTE S = XNEW - X, Y = GNEW - G, <S,S>, <S,Y>,
C     THE CONTINUOUS PROJECTED GRADIENT AND ITS NORMS

      sts = 0.0d0
      sty = 0.0d0
      do i = 1,n
          s(i) = xnew(i) - x(i)
          y(i) = gnew(i) - g(i)
          sts = sts + s(i)*s(i)
          sty = sty + s(i)*y(i)
          x(i) = xnew(i)
          g(i) = gnew(i)
          pg(i) = x(i) - g(i)
      end do

      call proj(n,pg,inform)

      if (inform .ne. 0) then
        ! ERROR IN PROJ SUBROUTINE, STOP
          flag = 4
          go to 200
      end if

      pgtwon = 0.0D0
      pginfn = 0.0D0
      do i = 1,n
          pg(i) = pg(i) - x(i)
          pgtwon = pgtwon + pg(i)**2
          pginfn = max(pginfn,abs(pg(i)))
      end do

C     PRINT ITERATION INFORMATION

      if (output) then
          write (*,fmt=9010) iter,f,pginfn
      end if

C     COMPUTE SPECTRAL STEPLENGTH

      if (sty .le. 0.0d0) then
          lambda = lmax

      else
          lambda = min(lmax,max(lmin,sts/sty))
      end if

C     FINISH THE ITERATION

      go to 100

C     STOP

 200  continue

C     SET X AND F WITH THE BEST SOLUTION AND ITS FUNCTION VALUE

      f = fbest
      do i = 1,n
          x(i) = xbest(i)
      end do

 9010 format ('ITER= ',I10,' F= ',1P,D17.10,' PGINFNORM= ',1P,D17.10)

      end


      subroutine ls(n,x,f,d,gtd,m,lastfv,maxfc,fcnt,fnew,xnew,flag)

C     Subroutine LS implements a nonmonotone line search with
C     safeguarded quadratic interpolation.
C
C     This version 17 JAN 2000 by E.G.Birgin, J.M.Martinez and M.Raydan.
C     Reformatted 03 OCT 2000 by Tim Hopkins.
C     Final revision 03 JUL 2001 by E.G.Birgin, J.M.Martinez and M.Raydan.
C
C     On Entry:
C
C     n     integer,
C           size of the problem,
C
C     x     double precision x(n),
C           initial guess,
C
C     f     double precision,
C           function value at the actual point,
C
C     d     double precision d(n),
C           search direction,
C
C     gtd   double precision,
C           internal product <g,d>, where g is the gradient at x,
C
C     m     integer,
C           number of previous function values to be considered 
C           in the nonmonotone line search,
C
C     lastfv double precision lastfv(m),
C           last m function values,
C
C     maxfc integer,
C           maximum number of function evaluations,
C
C     fcnt  integer,
C           actual number of function evaluations.
C
C     On Return:
C
C     fcnt  integer,
C           actual number of function evaluations,
C
C     fnew  double precision,
C           function value at the new point,
C
C     xnew  double precision xnew(n),
C           new point,
C
C     flag  integer,
C           0= convergence with nonmonotone Armijo-like criterion,
C           3= too many function evaluations,
C           5= error in evalf subroutine.

C     PARAMETERS
      double precision gamma
      parameter (gamma=1.0d-04)

C     SCALAR ARGUMENTS
      double precision f,fnew,gtd
      integer maxfc,fcnt,m,n,flag

C     ARRAY ARGUMENTS
      double precision d(n),lastfv(0:m-1),x(n),xnew(n)

C     LOCAL SCALARS
      double precision alpha,atemp,fmax
      integer i,inform

C     EXTERNAL SUBROUTINES
      external evalf

C     INTRINSIC FUNCTIONS
      intrinsic max

C     INITIALIZATION

C     COMPUTE THE MAXIMUM FUNCTIONAL VALUE OF THE LAST M ITERATIONS

      fmax = lastfv(0)
      do i = 1,m - 1
          fmax = max(fmax,lastfv(i))
      end do

C     COMPUTE FIRST TRIAL

      alpha = 1.0d0

      do i = 1,n
          xnew(i) = x(i) + d(i)
      end do

      call evalf(n,xnew,fnew,inform)
      fcnt = fcnt + 1

      if (inform .ne. 0) then
        ! ERROR IN EVALF SUBROUTINE, STOP
          flag = 5
          go to 200
      end if

C     MAIN LOOP

 100  continue

C     TEST STOPPING CRITERIA

      if (fnew .le. fmax + gamma*alpha*gtd) then
        ! NONMONOTONE ARMIJO-LIKE STOPPING CRITERION SATISFIED, STOP
          flag = 0
          go to 200
      end if

      if (fcnt .ge. maxfc) then
        ! MAXIMUM NUMBER OF FUNCTION EVALUATIONS EXCEEDED, STOP
          flag = 3
          go to 200
      end if 

C     DO AN ITERATION

C     SAFEGUARDED QUADRATIC INTERPOLATION

      if (alpha .le. 0.1d0) then
          alpha = alpha/2.0d0

      else
          atemp = (-gtd*alpha**2) / (2.0d0*(fnew-f-alpha*gtd))
          if (atemp .lt. 0.1d0 .or. atemp .gt. 0.9d0*alpha) then
              atemp = alpha/2.0d0
          end if
          alpha = atemp
      end if

C     COMPUTE TRIAL POINT 

      do i = 1,n
          xnew(i) = x(i) + alpha*d(i)
      end do

C     EVALUATE FUNCTION

      call evalf(n,xnew,fnew,inform)
      fcnt = fcnt + 1

      if (inform .ne. 0) then
        ! ERROR IN EVALF SUBROUTINE, STOP
          flag = 5
          go to 200
      end if

C     ITERATE

      go to 100

C     STOP

 200  continue

      end
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
#       End of shell archive
exit 0