* ------------------------------------------------------------------ * File maxi.f * * maxi mxdata funcon funmax * * This is a program to generate and solve optimization problems * of the form * * Given P, maximize d, the minimum diameter of P nonoverlapping * circles that will fit inside a unit square. * i.e. * maximize d2 * such that (x(p) - x(q))**2 + (y(p) - y(q))**2 >= d2, * p = 1..P, q = 1..p-1, * 0 <= x(p), y(p) <= 1, p = 1..P. * where d2 = d**2. * * maxi is the main program. * mxdata is the problem generator. * funcon defines the constraints and Jacobian. * funmax does the work for funcon. * snopt (subroutine version of SNOPT 4.00) gets called to * solve the problem. * * * 28 Aug 1992: First version, based on AMPL formulation sent to * MAS by David Gay. MINOS 5.4 (Jul 1992) had trouble. * 30 May 1994: Modified to call SNOPT by PEG. * 27 Nov 1996: Modified for SNOPT versions 4.3 and higher. * 26 Sep 1997: Modified for SNOPT versions 5.3 and higher. * ------------------------------------------------------------------ program maxi implicit double precision (a-h,o-z) parameter ( maxP = 50, $ maxm = (maxP - 1)*maxP / 2 + 1, $ maxn = 2*maxP + 1, $ maxne = 2*(maxP - 1)*maxP + maxm, $ nname = 1 ) character*8 Prob character*8 Names(nname) integer ha(maxne) , hs(maxm+maxn) integer ka(maxn+1) double precision a(maxne) , bl(maxm+maxn), bu(maxm+maxn), $ x(maxm+maxn), pi(maxm) , rc(maxm+maxn) parameter ( leniw = 100000 + 50*maxP) parameter ( lenrw = 100000 + 100*maxP) parameter ( lencw = 500 ) character*8 cw(lencw) integer iw(leniw) double precision rw(lenrw) integer P character*4 charP logical byname character*20 lfile external funcon, funobj * ------------------------------------------------------------------ * Specify some of the SNOPT files. * iSpecs is the Specs file (0 if none). * iPrint is the Print file (0 if none). * iSumm is the Summary file (0 if none). * nout is an output file used here. iSpecs = 4 iPrint = 15 iSumm = 6 nout = 6 byname = .true. if ( byname ) then * Unix and DOS systems. Open the Specs and print files. lfile = 'maxi.spc' open( iSpecs, file=lfile, status='OLD', err=800 ) lfile = 'maxi.out' open( iPrint, file=lfile, status='UNKNOWN', err=800 ) end if * ------------------------------------------------------------------ * Set options to their default values. * ------------------------------------------------------------------ call snInit( iPrint, iSumm, $ cw, lencw, iw, leniw, rw, lenrw ) * ------------------------------------------------------------------ * Read a Specs file. * ------------------------------------------------------------------ call snSpec( iSpecs, inform, $ cw, lencw, iw, leniw, rw, lenrw ) if (inform .ge. 2) then write(nout, *) 'iSpecs > 0 but no Specs file found' stop end if * ------------------------------------------------------------------ * Set P to be the number of circles. We assume the Specs file * does this via the line * Problem number P * The following call fetches nProb, which defines P. * ------------------------------------------------------------------ call sngeti( 'Problem number', nProb, $ inform, cw, lencw, iw, leniw, rw, lenrw ) P = nProb write(charP, '(i4)') P * Give names to the Problem, Objective, Rhs, Ranges and Bounds. Prob = 'maxi' // charP call mxdata( P, maxm, maxn, maxne, inform, $ m, n, ne, nnCon, nnObj, nnJac, $ a, ha, ka, bl, bu, x, pi, hs ) if (inform .ge. 1) then write(nout, *) 'Not enough storage to generate a problem ', $ 'with P =', P stop end if * ------------------------------------------------------------------ * Specify options that were not set in the Specs file. * i1 and i2 may refer to the Print and Summary file respectively. * Setting them to 0 suppresses printing. * ------------------------------------------------------------------ i1 = 0 i2 = 0 ltime = 2 call snseti( 'Timing level ', ltime, i1, i2, inform, $ cw, lencw, iw, leniw, rw, lenrw ) * ------------------------------------------------------------------ * Go for it, using a Cold start. * iObj = m specifies the linear objective row in a(*). * Objadd = 0.0 means there is no constant to be added to the * objective. * hs need not be set if a basis file is to be input. * Otherwise, each hs(1:n) should be 0, 1, 2, 3, 4, or 5. * The values are used by the Crash procedure m2crsh * to choose an initial basis B. * If hs(j) = 0 or 1, column j is eligible for B. * If hs(j) = 2, column j is initially superbasic (not in B). * If hs(j) = 3, column j is eligible for B and is given * preference over columns with hs(j) = 0 or 1. * If hs(j) = 4 or 5, column j is initially nonbasic. * ------------------------------------------------------------------ iObj = m ObjAdd = 0.0d0 call snopt ( 'Cold', m, n, ne, nName, $ nnCon, nnObj, nnJac, $ iObj, ObjAdd, Prob, $ funcon, funobj, $ a, ha, ka, bl, bu, Names, $ hs, x, pi, rc, $ inform, mincw, miniw, minrw, $ nS, nInf, sInf, Obj, $ cw, lencw, iw, leniw, rw, lenrw, $ cw, lencw, iw, leniw, rw, lenrw ) if (inform .eq. 42 .or. inform .eq. 43 .or. inform .eq. 44) then write(nout, *) ' ' write(nout, *) 'Estimate of required lencw:', mincw write(nout, *) 'Estimate of required leniw:', miniw write(nout, *) 'Estimate of required lenrw:', minrw go to 900 end if write(nout, *) ' ' write(nout, *) 'snopt finished.' write(nout, *) 'inform =', inform write(nout, *) 'nInf =', nInf write(nout, *) 'sInf =', sInf write(nout, *) 'Obj =', Obj if (inform .ge. 20) go to 900 stop * ------------------------------------------------------------------ * Error exit. * ------------------------------------------------------------------ 800 write(nout, 4000) 'Error while opening file', lfile stop 900 write(nout, *) ' ' write(nout, *) 'STOPPING because of error condition' stop 4000 format(/ a, 2x, a ) * end of main program maxi end