C Search for optimal packings circles in circle,
C optionally with different radii
C Author: Hugo Pfoertner http://www.pfoertner.org
C
C Version history:
C Jan 6 2006 Simplified single precision version
C            with omission of small circles
C Oct 22 2005 Output format for coordinate files changed
C Oct 21 2005 Optional rotation of coordinates read from file
C Oct 14 2005 Double precision
C Sep 30 2005 Additional output of inverse radius
C Jun 16 2005 Index error in single point move phase corrected,
C              inhibit moves with negative total effect.
C Jun 14 2005 Initial version, derived from sphinsph.for
      PARAMETER ( ZERO=0.0, ONE=1.0, TWO=2.0, VLARGE=1.0E35 )
      REAL MINCIR, MINCI1
      EXTERNAL MINCIR, MINCI1, FMINSI
CWatcom
      DOUBLEPRECISION GENRAND_REAL2
*$pragma aux genrand_real2 "genrand_real2_"
      PARAMETER ( MM=50 )
C NMAX = MM/3
      PARAMETER (NMAX=2*MM)
      REAL XX(NMAX), XS(NMAX), EPS(NMAX)
      REAL                      X(MM), Y(MM)
      COMMON /COORD/ MSTART, M, X,     Y
C Target sphere during single point shift phase of optimization
      COMMON /TARGET/ NTARGT
      CHARACTER FNAM*10
C
C Local optimization variables
      REAL  X3(2)
CAltix
CAltix      INTEGER SEED(2)
C STATEMENT FUNCTION
C***      D2(I,J) = SQRT((X(I)-X(J))**2+(Y(I)-Y(J))**2)
      FNAM = 'WCIC00.TXT'
C Name and version of program
      write (*,*) ' Program UCIRPCK, Version Date: January 6 2006'
C
C RANDOM LOOP
      M = 1
      LOOPS = 1
      BESTR = ZERO
      WRITE (*,*) 'N, LOOPS'
      READ (*,*) MREAD, LOOPS
      if ( abs(mread) .gt. 99 ) stop 'N exceeds maximum 99.'
      MD = MREAD
C Start big loop,
C including potentially multiple re-reading of config files
C
      DBEST = ZERO
      DO 900 LOOP = 1, LOOPS
      M = ABS ( MD )
C Small circles with sizes < MSTART are excluded from search
      MSTART = NINT ( 0.3 * FLOAT(M) )
C FREE VARIABLES
      N = M + M
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
      WRITE (FNAM(5:6),1002) M
1002  FORMAT ( I2.2)
      OPEN ( UNIT=10, FORM='FORMATTED', FILE=FNAM,
     &       STATUS='OLD', IOSTAT=IOS )
      IF ( IOS .NE. 0 ) GOTO 90
      REWIND ( UNIT=10 )
      READ (10,*) M, DBEST
      IF ( ABS(MD) .NE. M ) STOP 'NON-FITTING M'
      IF ( MD .LT. 0 ) THEN
C READ COORDINATES FROM FILE (NORMALIZED POSITION IS ASSUMED)
        L = 0
        DO 5 K = MSTART, M
        READ (10,*) I, X(K), Y(K)
        L = L + 1
        XX(L) = X(K)
        L = L + 1
        XX(L) = Y(K)
5       CONTINUE
        CLOSE ( UNIT=10 )
6       CONTINUE
        DBEST = -MINCIR(N,XX)
        WRITE (*,*) ' BEST SO FAR:', DBEST
C SKIP RANDOM SEARCH FOR START CONFIG
        GOTO 101
      ENDIF
      CLOSE ( UNIT=10 )
90    CONTINUE
C ZUFALLSZAHLENGENERATOR INITIALISIEREN
      IF ( LOOP .EQ. 1 ) then
        write (*,*) ' Seed Mersenne Twister.'
        CALL MTINIT
CAltix
CAltix        T0 = SECNDS ( 0.0 )
CAltix        SEED(1) = NINT ( 100.0 * T0 )
CAltix        SEED(2) = 1
CAltix        CALL RANDOM_SEED ( PUT=SEED(1:2) )
CAltix
      endif
C
      D = ZERO
C
C START LOOP OVER RANDOM CONFIGS
C REPETIONS
      write (*,*) ' '
      write (*,1001) loop, loops, M
1001  format ( ' Loop', I6, ' of', I6, ' for N=', I3 )
      write (*,*) ' Random Search:'
      write (*,*) ' Configuration   Best Radius'
      MTRY = 1000000
      ITRY = 0
100   CONTINUE
      ITRY = ITRY + 1
      L = 0
      LRAN = 0
      do 110 np = MSTART, M
111   continue
CAltix      IF ( LRAN .LT. 3 ) THEN
CAltix        LRAN = N+1
CAltix
CAltix        CALL RANDOM_NUMBER ( HARVEST=EPS(1:N) )
CAltix
CAltix      ENDIF
      XP = TWO * GENRAND_REAL2 () - ONE
      YP = TWO * GENRAND_REAL2 () - ONE
CAltix
CAltix      LRAN = LRAN - 1
CAltix      XP = EPS(LRAN) + EPS(LRAN) - 1.0
CAltix      LRAN = LRAN - 1
CAltix      YP = EPS(LRAN) + EPS(LRAN) - 1.0
CAltix
C Reject points outside container sphere
      IF ( XP*XP + YP*YP .GT. ONE ) goto 111
      L = L + 1
      XX(L) = XP
      L = L + 1
      XX(L) = YP
110   CONTINUE
C
      DTRY = MINCIR ( N, XX )
      IF ( DTRY .LT. D ) THEN
C SAVE BEST CONFIG FOUND SO FAR
        DO 130 L = 1, N
        XS(L) = XX(L)
130     CONTINUE
        D = DTRY
        WRITE (*,*) ITRY, -D
      ENDIF
C
C TRY NEXT CONFIG
      IF ( ITRY .LE. MTRY ) GOTO 100
C
C WRITE COORDINATES OF BEST CONFIGS
C RESTORE
        DO 140 L = 1, N
        XX(L) = XS(L)
140     CONTINUE
C
C TARGET LABEL IF START CONFIG IS READ FROM FILE
101   CONTINUE
C
      D = MINCIR ( N, XX )
C***  write (*,*) ' BESTR=', BESTR
      WRITE (*,*) ' Non-Linear Optimization:'
      write (*,*) '     Restart     Best Radius'
C
C BESETZUNG DER GENAUIGKEITSSCHRANKEN
      DO 24 I = 1, N
      EPS(I) = ZERO
24    CONTINUE
C
C MEHRFACHE OPTIMIERUNGSVERSUCHE
C
      FALT = VLARGE
      NREP = 0
      NOPFIN = 0
      MOVFIN = 0
C
C Target label after single point moves and exchanges
299   continue
      NREP = NREP + 1
      if ( NREP .GT. M .or. 
     &     (nopfin .eq. 1 .and. movfin .eq. 1000*M+1) ) goto 900
C
      DO 300 NOPT = 1, 100000
      IER = 0
C
      CALL FMINSI ( N, XX, MINCIR, EPS, START, IER )
C
C ABFRAGE AUF VERBESSERUNG
      IF ( START .LT. FALT ) THEN
C
        FALT = START
        DO 203 I = 1, N
        XS(I) = XX(I)
203     CONTINUE
      IF ( NOPT .EQ. 1 .OR. MOD (NOPT,500) .EQ. 0 )
     &  WRITE (*,*) NOPT, ABS(START)
      ELSE
C
C KEINE WEITERE VERBESSERUNG ERREICHT
        NOPFIN = NOPT
        WRITE (*,*) ' No improvement after', NOPT, ' restarts'
        GOTO 202
      ENDIF
C
C ENDE DER SCHLEIFE UEBER MEHRFACHE OPTIMIERUNGSVERSUCHE
300   CONTINUE
202   CONTINUE
      START = -MINCIR (N, XS)
      BESTR = MAX ( BESTR, START )
      WRITE (*,*) ' Best radius found:', START
      IF ( START .GE. DBEST ) THEN
        OPEN ( UNIT=10, FORM='FORMATTED', FILE=FNAM,
     &         STATUS='UNKNOWN', IOSTAT=IOS )
        WRITE (*,*) ' IMPROVED:', DBEST, '->', START
        REWIND (10)
        WRITE (10,*) M, START, ONE/START
        DO 310 I = MSTART, M
        WRITE (10,1010) I, X(I), Y(I)
1010    FORMAT ( I3, 2F21.17 )
310     CONTINUE
C        DO 311 K = 1, M
C        XBEST(K) = X(K)
C        YBEST(K) = Y(K)
C311     CONTINUE
        CLOSE (UNIT=10)
        DBEST = START
      ENDIF
C
C
Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
C Move single points to maximally distant locations from neighbors
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
C      write (*,*) ' Move single points:'
      do 401 i = 1, n
      xx(i) = xs(i)
401   continue
C Best result found by shaking
      FALT = MINCIR ( n, xx )
C***  write (*,*) ' FALT=', FALT
C
C Counter for print frequency
      NOPTP = 0
C      lastsc = 0
C Counter for unsuccessful passes
      nosuc = 0
C
      DO 400 NOPT = 1, 100000*M
      NOPTP = NOPTP + 1
C
      deltot = ZERO
      do 410 NTARGT = MSTART, M
      K1 = 2*(NTARGT-MSTART)+1
      K2 = K1 + 1
      x3(1) = xx(K1)
      x3(2) = xx(K2)
      FLOC = MINCI1(2,x3)
      xxs = x3(1)
      yys = x3(2)
      IER = 0
      START = FLOC
      CALL FMINSI ( 2, X3, MINCI1, EPS, START, IER )
      if ( IER .eq. 0 .and. start .lt. FLOC ) then
C overwrite global coordinates
        xx(K1) = x3(1)
        xx(K2) = x3(2)
        dummy = MINCI1(2,x3)
        if ( dummy .gt. FALT ) then
C Undo move
          xx(K1) = xxs
          xx(K2) = yys
        else
C accumulate move distance
          if ( dummy .lt. falt )
     &        DELTOT = DELTOT + abs(x3(1)-xxs) + abs(x3(2)-yys)

        endif
C Diag
C***    write (*,1011) ntargt, start, floc, minsph(n,xx)
C***1011    format ( I3, 3F14.9 )
      endif
C End of optimization loop over single points
410   continue
C terminate if all points stay at their previous locations
      if ( DELTOT .EQ. ZERO ) goto 402
C Check if global distance has been improved
      start = MINCIR (n,xx)
C
      IF ( START .LT. FALT ) THEN
C
        FALT = START
        DO 403 I = 1, N
        XS(I) = XX(I)
403     CONTINUE
C print result of first step
        IF ( NOPT .EQ. 1 ) write (*,*) nopt, deltot, abs(start)
C print next step after a multiple of 10000
        IF ( NOPTP .GE. 10000 ) THEN
          write (*,*) nopt, deltot, abs(start)
          NOPTP = NOPTP - 10000
404       continue
          if ( NOPTP .GE. 10000 ) THEN
            NOPTP = NOPTP - 10000
            GOTO 404
          endif
        ENDIF
C        if ( nopt - lastsc .gt. nosuc ) then
C          nosuc = nopt - lastsc
C          write (*,*) nopt, nosuc
C        endif
C        lastsc = nopt
C Reset termination counter
        nosuc = 0
      ELSE
C
C KEINE WEITERE VERBESSERUNG ERREICHT
        nosuc = nosuc + 1
C Check for termination
        if ( nosuc .gt. 1000*M ) goto 402
C        WRITE (*,*) ' No improvement after', NOPT, ' restarts'
C        GOTO 402
      ENDIF
C
C ENDE DER SCHLEIFE UEBER MEHRFACHE OPTIMIERUNGSVERSUCHE
400   CONTINUE
402   CONTINUE
      write (*,*) ' Single point moves: ', nopt
      MOVFIN = NOPT
      START = -MINCIR (N, XS)
      WRITE (*,*) ' Best radius found:', START
      FIMPRV = START
C Check if circles have free space
      do 500 NTARGT = M-1, MSTART, -1
      K1 = 2*(NTARGT-MSTART)+1
      K2 = K1 + 1
      x3(1) = xs(K1)
      x3(2) = xs(K2)
      FLOC0 = MINCI1(2,x3)
      IER = 0
      CALL FMINSI ( 2, X3, MINCI1, EPS, FLOC, IER )
      if ( IER .eq. 0 .and. FLOC .le. FLOC0 ) then
        KTGT = min(INT(-NTARGT*FLOC/START),M)
C        write (*,*) ntargt, -FLOC/START, -NTARGT*FLOC/START, KTGT
        if ( KTGT .ge. NTARGT+1 ) then
C Exchange of circle locations
        j1 = 2*(KTGT-MSTART)+1
        xxs = xs(j1)
        yys = xs(j1+1)
        xs(k1) = xs(j1)
        xs(k2) = xs(j1+1)
        xs(j1) = x3(1)
        xs(j1+1) = x3(2)
        if ( -mincir(n,xs) .ge. START ) then
          write (*,*) ' Exchange ', NTARGT, KTGT
          goto 501
        else
          xs(j1) = xxs
          xs(j1+1) = yys
          xs(k1) = x3(1)
          xs(k2) = x3(2)
        endif
        endif
      else
        write (*,*) ntargt, FLOC0/START, -NTARGT*FLOC0/START
      endif
500   CONTINUE
501   CONTINUE
C
      BESTR = MAX ( BESTR, FIMPRV )
      IF ( START .GT. DBEST ) THEN
        OPEN ( UNIT=10, FORM='FORMATTED', FILE=FNAM,
     &         STATUS='UNKNOWN', IOSTAT=IOS )
        WRITE (*,*) ' IMPROVED:', DBEST, '->', FIMPRV
        REWIND (10)
        WRITE (10,*) M, FIMPRV, ONE/FIMPRV
        DO 450 I = MSTART, M
        WRITE (10,1010) I, X(I), Y(I)
450     CONTINUE
        CLOSE (UNIT=10)
        write (*,*) ' Results saved to file ', FNAM
        DBEST = START
      ENDIF
C Repeat optimization of all circle locations
      do 460 I = 1, N
      xx(i) = xs(i)
460   continue
      GOTO 299
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
C END OF BIG LOOP
900   CONTINUE
      END