PROCEDURE adi(a,b,c,d,e,f,g: gljmax; VAR u: gljmax;
      jmax,k: integer; alpha,beta,eps: double);
(* Programs using routine ADI must define the type
TYPE
   gljmax = ARRAY [1..jmax,1..jmax] OF double;
in the main routine. *)
LABEL 99;
CONST
   jj=50;
   kk=6;
   nrr=32;      (* nrr=2 to the power (kk-1) *)
   maxits=100;
   zero=0.0;
   two=2.0;
   half=0.5;
TYPE
   gljjarray = ARRAY [1..jj] OF double;
VAR
   i,nr,nits,next,n,l,kits,k1,j,twopwr: integer;
   rfact,resid,disc,anormg,anorm,ab: double;
   aa,bb,cc,rr,uu: gljjarray;
   psi: ARRAY [1..jj,1..jj] OF double;
   alph,bet: ARRAY [1..kk] OF double;
   r: ARRAY [1..nrr] OF double;
   s: ARRAY [1..nrr,1..kk] OF double;
PROCEDURE tridag(a,b,c,r: gljjarray; VAR u: gljjarray; n: integer);
   (* This is a double precision version of TRIDAG for use with ADI,
      which defines the constant jj and the type gljjarray.   *)
   VAR
      j: integer;
      bet: double;
      gam: gljjarray;
   BEGIN
      IF (b[1] = 0.0) THEN BEGIN
         writeln('pause 1 in TRIDAG'); readln END;
      bet := b[1];
      u[1] := r[1]/bet;
      FOR j := 2 to n DO BEGIN
         gam[j] := c[j-1]/bet;
         bet := b[j]-a[j]*gam[j];
         IF (bet = 0.0) THEN BEGIN
            writeln('pause 2 in TRIDAG'); readln END;
         u[j] := (r[j]-a[j]*u[j-1])/bet
      END;
      FOR j := n-1 DOWNTO 1 DO u[j] := u[j]-gam[j+1]*u[j+1]
   END;
BEGIN
   IF (jmax > jj) THEN BEGIN
      writeln('Pause in routine ADI - increase jj'); readln
   END;
   IF (k > (kk-1)) THEN BEGIN
      writeln('Pause in routine ADI - increase kk'); readln
   END;
   k1 := k+1;
   nr := 1;
   FOR i := 1 to k DO nr := 2*nr;
   alph[1] := alpha;
   bet[1] := beta;
   FOR j := 1 to k DO BEGIN
      alph[j+1] := sqrt(alph[j]*bet[j]);
      bet[j+1] := half*(alph[j]+bet[j])
   END;
   s[1,1] := sqrt(alph[k1]*bet[k1]);
   FOR j := 1 to k DO BEGIN
      ab := alph[k1-j]*bet[k1-j];
      twopwr := 1;
      FOR i := 1 to (j-1) DO twopwr := 2*twopwr;
      FOR n := 1 to twopwr DO BEGIN
         disc := sqrt(sqr(s[n,j])-ab);
         s[2*n,j+1] := s[n,j]+disc;
         s[2*n-1,j+1] := ab/s[2*n,j+1]
      END
   END;
   FOR n := 1 to nr DO r[n] := s[n,k1];
   anormg := zero;
   FOR j := 2 to (jmax-1) DO BEGIN
      FOR l := 2 to (jmax-1) DO BEGIN
         anormg := anormg+abs(g[j,l]);
         psi[j,l] := -d[j,l]*u[j,l-1]
            +(r[1]-e[j,l])*u[j,l]-f[j,l]*u[j,l+1]
      END
   END;
   nits := maxits DIV nr;
   FOR kits := 1 to nits DO BEGIN
      FOR n := 1 to nr DO BEGIN
         IF (n = nr) THEN next := 1 ELSE next := n+1;
         rfact := r[n]+r[next];
         FOR l := 2 to (jmax-1) DO BEGIN
            FOR j := 2 to jmax-1 DO BEGIN
               aa[j-1] := a[j,l];
               bb[j-1] := b[j,l]+r[n];
               cc[j-1] := c[j,l];
               rr[j-1] := psi[j,l]-g[j,l]
            END;
            tridag(aa,bb,cc,rr,uu,jmax-2);
            FOR j := 2 to (jmax-1) DO BEGIN
               psi[j,l] := -psi[j,l]
                  +two*r[n]*uu[j-1]
            END
         END;
         FOR j := 2 to (jmax-1) DO BEGIN
            FOR l := 2 to (jmax-1) DO BEGIN
               aa[l-1] := d[j,l];
               bb[l-1] := e[j,l]+r[n];
               cc[l-1] := f[j,l];
               rr[l-1] := psi[j,l]
            END;
            tridag(aa,bb,cc,rr,uu,jmax-2);
            FOR l := 2 to (jmax-1) DO BEGIN
               u[j,l] := uu[l-1];
               psi[j,l] := -psi[j,l]+rfact*uu[l-1]
            END
         END
      END;
      anorm := zero;
      FOR j := 2 to (jmax-1) DO BEGIN
         FOR l := 2 to (jmax-1) DO BEGIN
            resid := a[j,l]*u[j-1,l]+(b[j,l]+e[j,l])*u[j,l]
               +c[j,l]*u[j+1,l]+d[j,l]*u[j,l-1]
               +f[j,l]*u[j,l+1]+g[j,l];
            anorm := anorm+abs(resid)
         END
      END;
      IF (anorm < (eps*anormg)) THEN GOTO 99
   END;
   writeln('Pause in routine ADI - too many iterations');
99:   END;
PROCEDURE amoeba(VAR p: glmpnp; VAR y: glmp; ndim: integer;
       ftol: real; VAR iter: integer);
(* Programs using routine AMOEBA must supply an external function
func(pr:glnp):real whose minimum is to be found. They must
also define types
TYPE
   glmpnp = ARRAY [1..mp,1..np] OF real;
   glmp = ARRAY [1..mp] OF real;
   glnp = ARRAY [1..np] OF real;
where mp and np are physical dimensions *)
LABEL 99;
CONST
   alpha=1.0;
   beta=0.5;
   gamma=2.0;
   itmax=500;
VAR
   mpts,j,inhi,ilo,ihi,i: integer;
   yprr,ypr,rtol: real;
   pr,prr,pbar: glnp;
BEGIN
   mpts := ndim+1;
   iter := 0;
   WHILE true DO BEGIN
      ilo := 1;
      IF (y[1] > y[2]) THEN BEGIN
         ihi := 1;
         inhi := 2
      END ELSE BEGIN
         ihi := 2;
         inhi := 1
      END;
      FOR i := 1 to mpts DO BEGIN
         IF (y[i] < y[ilo]) THEN  ilo := i;
         IF (y[i] > y[ihi]) THEN BEGIN
            inhi := ihi;
            ihi := i
         END ELSE IF  (y[i] > y[inhi]) THEN
            IF (i <> ihi) THEN  inhi := i
      END;
      rtol := 2.0*abs(y[ihi]-y[ilo])/(abs(y[ihi])+abs(y[ilo]));
      IF (rtol < ftol) THEN GOTO 99;
      IF (iter = itmax) THEN BEGIN
         writeln('pause in AMOEBA - too many iterations'); readln
      END;
      iter := iter+1;
      FOR j := 1 to ndim DO pbar[j] := 0.0;
      FOR i := 1 to mpts DO
         IF (i <> ihi) THEN FOR j := 1 to ndim DO pbar[j] := pbar[j]+p[i,j];
      FOR j := 1 to ndim DO BEGIN
         pbar[j] := pbar[j]/ndim;
         pr[j] := (1.0+alpha)*pbar[j]-alpha*p[ihi,j]
      END;
      ypr := func(pr);
      IF (ypr <= y[ilo]) THEN BEGIN
         FOR j := 1 to ndim DO prr[j] := gamma*pr[j]+(1.0-gamma)*pbar[j];
         yprr := func(prr);
         IF (yprr < y[ilo]) THEN BEGIN
            FOR j := 1 to ndim DO p[ihi,j] := prr[j];
            y[ihi] := yprr
         END ELSE BEGIN
            FOR j := 1 to ndim DO p[ihi,j] := pr[j];
            y[ihi] := ypr
         END
      END ELSE IF  (ypr >= y[inhi]) THEN BEGIN
         IF (ypr < y[ihi]) THEN BEGIN
            FOR j := 1 to ndim DO p[ihi,j] := pr[j];
            y[ihi] := ypr
         END;
         FOR j := 1 to ndim DO prr[j] := beta*p[ihi,j]+(1.0-beta)*pbar[j];
         yprr := func(prr);
         IF (yprr < y[ihi]) THEN BEGIN
            FOR j := 1 to ndim DO p[ihi,j] := prr[j];
            y[ihi] := yprr
         END ELSE BEGIN
            FOR i := 1 to mpts DO
               IF (i <> ilo) THEN BEGIN
                  FOR j := 1 to ndim DO BEGIN
                     pr[j] := 0.5*(p[i,j]+p[ilo,j]);
                     p[i,j] := pr[j]
                  END;
                  y[i] := func(pr)
               END
         END
      END ELSE BEGIN
         FOR j := 1 to ndim DO p[ihi,j] := pr[j];
         y[ihi] := ypr
      END
   END;
99:   END;
PROCEDURE anneal(x,y : cityarray; VAR iorder: iarray; ncity: integer);
(* Programs using routine ANNEAL must define types
   cityarray : ARRAY [1..ncity] OF real;
   iarray : ARRAY [1..ncity] OF integer;
in the main routine. *)
LABEL 10,20,99;
CONST
   tfactr = 0.9;
TYPE
   nsix = ARRAY [1..6] OF integer;
VAR
   ans : boolean;
   path,de,t : real;
   nover,nlimit,i1,i2,idum,iseed: integer;
   i,j,k,nsucc,nn,idec : integer;
   n : nsix;
 
FUNCTION alen(x1,x2,y1,y2: real): real;
BEGIN
   alen := sqrt(sqr(x2-x1)+sqr(y2-y1))
END;

PROCEDURE revcst(x,y: cityarray; iorder: iarray; ncity: integer;
      VAR n: nsix; VAR de: real);
VAR
   xx,yy : ARRAY [1..6] OF real;
   j,ii : integer;
BEGIN
   n[3] := 1 + ((n[1]+ncity-2) MOD ncity);
   n[4] := 1 + (n[2] MOD ncity);
   FOR j := 1 to 4 DO BEGIN
      ii := iorder[n[j]];
      xx[j] := x[ii];
      yy[j] := y[ii]
   END;
   de := -alen(xx[1],xx[3],yy[1],yy[3])-alen(xx[2],xx[4],yy[2],yy[4])
      +alen(xx[1],xx[4],yy[1],yy[4])+alen(xx[2],xx[3],yy[2],yy[3])
END;

PROCEDURE reverse(VAR iorder: iarray; ncity: integer; n: nsix);
VAR
   nn,j,k,l,itmp : integer;
BEGIN
   nn := (1+((n[2]-n[1]+ncity) MOD ncity)) DIV 2;
   FOR j := 1 to nn DO BEGIN
      k := 1 + ((n[1]+j-2) MOD ncity);
      l := 1 + ((n[2]-j+ncity) MOD ncity);
      itmp := iorder[k];
      iorder[k] := iorder[l];
      iorder[l] := itmp
   END
END;

PROCEDURE trncst(x,y: cityarray; iorder: iarray; ncity: integer;
      VAR n: nsix; VAR de: real);
VAR
   xx,yy : ARRAY [1..6] OF real;
   j,ii : integer;
BEGIN
   n[4] := 1 + (n[3] MOD ncity);
   n[5] := 1 + ((n[1]+ncity-2) MOD ncity);
   n[6] := 1 + (n[2] MOD ncity);
   FOR j := 1 to 6 DO BEGIN
      ii := iorder[n[j]];
      xx[j] := x[ii];
      yy[j] := y[ii]
   END;
   de := -alen(xx[2],xx[6],yy[2],yy[6])-alen(xx[1],xx[5],yy[1],yy[5])
      -alen(xx[3],xx[4],yy[3],yy[4])+alen(xx[1],xx[3],yy[1],yy[3])
      +alen(xx[2],xx[4],yy[2],yy[4])+alen(xx[5],xx[6],yy[5],yy[6])
END;
 
PROCEDURE trnspt(VAR iorder: iarray; ncity: integer; n: nsix);
CONST
   maxcity=1000;
VAR
   jorder : ARRAY [1..maxcity] OF integer;
   m1,m2,m3,nn,j,jj : integer;
BEGIN
   m1 := 1 + ((n[2]-n[1]+ncity) MOD ncity);
   m2 := 1 + ((n[5]-n[4]+ncity) MOD ncity);
   m3 := 1 + ((n[3]-n[6]+ncity) MOD ncity);
   nn := 1;
   FOR j := 1 to m1 DO BEGIN
      jj := 1 + ((j+n[1]-2) MOD ncity);
      jorder[nn] := iorder[jj];
      nn := nn+1
   END;
   IF (m2>0) THEN BEGIN
      FOR j := 1 to m2 DO BEGIN
         jj := 1+((j+n[4]-2) MOD ncity);
         jorder[nn] := iorder[jj];
         nn := nn+1
      END
   END;
   IF (m3>0) THEN BEGIN
      FOR j := 1 to m3 DO BEGIN
         jj := 1 + ((j+n[6]-2) MOD ncity);
         jorder[nn] := iorder[jj];
         nn := nn+1
      END
   END;
   FOR j := 1 to ncity DO BEGIN
      iorder[j] := jorder[j] 
   END
END;
 
PROCEDURE metrop(de,t: real; VAR ans: boolean);
(* Programs using routine METROP must declare the variable 
VAR
   gljdum : integer;
and initialize its value to
   gljdum := 1;
in the main routine. *)
BEGIN
   ans := (de<0.0) OR (ran3(gljdum)<exp(-de/t))
END;
 
BEGIN
   nover := 100*ncity;
   nlimit := 10*ncity;
   path := 0.0;
   t := 0.5;
   FOR i := 1 to (ncity-1) DO BEGIN
      i1 := iorder[i];
      i2 := iorder[i+1];
      path := path+alen(x[i1],x[i2],y[i1],y[i2])
   END;
   i1 := iorder[ncity];
   i2 := iorder[1];
   path := path+alen(x[i1],x[i2],y[i1],y[i2]);
   idum := -1;
   iseed := 111;
   FOR j := 1 to 100 DO BEGIN
      nsucc := 0;
      FOR k := 1 to nover DO BEGIN
10:      n[1] := 1+trunc(ncity*ran3(idum));
         n[2] := 1+trunc((ncity-1)*ran3(idum));
         IF (n[2]>=n[1]) THEN n[2] := n[2]+1;
         nn := 1+((n[1]-n[2]+ncity-1) MOD ncity);
         IF (nn<3) THEN goto 10;
         idec := irbit1(iseed);
         IF (idec=0) THEN BEGIN
            n[3] := n[2]+trunc(abs(nn-2)*ran3(idum))+1;
            n[3] := 1+((n[3]-1) MOD ncity);
            trncst(x,y,iorder,ncity,n,de);
            metrop(de,t,ans);
            IF ans THEN BEGIN
               nsucc := nsucc+1;
               path := path+de;
               trnspt(iorder,ncity,n)
            END
         END ELSE BEGIN
            revcst(x,y,iorder,ncity,n,de);
            metrop(de,t,ans);
            IF ans THEN BEGIN
               nsucc := nsucc+1;
               path := path+de;
               reverse(iorder,ncity,n)
            END
         END;
         IF (nsucc>=nlimit) THEN goto 20
      END;
20:   writeln;
      writeln('T =',t:10:6,'    Path Length =',path:12:6);
      writeln('Successful Moves: ',nsucc:6);
      t := t*tfactr;
      IF (nsucc=0) THEN goto 99
   END;
99:
END;
PROCEDURE avevar(data: glnparray; n: integer; VAR ave,svar: real);
(* Programs using routine AVEVAR must define type
TYPE
   glnparray = ARRAY [1..np] OF real;
in the main routine, with np >= n.   *)
VAR
   j: integer;
   s: real;
BEGIN
   ave := 0.0;
   svar := 0.0;
   FOR j := 1 to n DO BEGIN
      ave := ave+data[j];
   END;
   ave := ave/n;
   FOR j := 1 to n DO BEGIN
      s := data[j]-ave;
      svar := svar+s*s
   END;
   svar := svar/(n-1)
END;
