      SUBROUTINE ADI(A,B,C,D,E,F,G,U,JMAX,K,ALPHA,BETA,EPS)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(JJ=100,KK=6,NRR=32,MAXITS=100,ZERO=0.D0,TWO=2.D0,HALF=.5
     *D0)
      DIMENSION A(JMAX,JMAX),B(JMAX,JMAX),C(JMAX,JMAX),D(JMAX,JMAX),
     *    E(JMAX,JMAX),F(JMAX,JMAX),G(JMAX,JMAX),U(JMAX,JMAX),
     *    AA(JJ),BB(JJ),CC(JJ),RR(JJ),UU(JJ),PSI(JJ,JJ),
     *    ALPH(KK),BET(KK),R(NRR),S(NRR,KK)
      IF(JMAX.GT.JJ)PAUSE 'Increase JJ'
      IF(K.GT.KK-1)PAUSE 'Increase KK'
      K1=K+1
      NR=2**K
      ALPH(1)=ALPHA
      BET(1)=BETA
      DO 11 J=1,K
        ALPH(J+1)=SQRT(ALPH(J)*BET(J))
        BET(J+1)=HALF*(ALPH(J)+BET(J))
11    CONTINUE
      S(1,1)=SQRT(ALPH(K1)*BET(K1))
      DO 13 J=1,K
        AB=ALPH(K1-J)*BET(K1-J)
        DO 12 N=1,2**(J-1)
          DISC=SQRT(S(N,J)**2-AB)
          S(2*N,J+1)=S(N,J)+DISC
          S(2*N-1,J+1)=AB/S(2*N,J+1)
12      CONTINUE
13    CONTINUE
      DO 14 N=1,NR
        R(N)=S(N,K1)
14    CONTINUE
      ANORMG=ZERO
      DO 16 J=2,JMAX-1
        DO 15 L=2,JMAX-1
          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)
15      CONTINUE
16    CONTINUE
      NITS=MAXITS/NR
      DO 27 KITS=1,NITS
        DO 24 N=1,NR
          IF(N.EQ.NR)THEN
            NEXT=1
          ELSE
            NEXT=N+1
          ENDIF
          RFACT=R(N)+R(NEXT)
          DO 19 L=2,JMAX-1
            DO 17 J=2,JMAX-1
              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)
17          CONTINUE
            CALL TRIDAG(AA,BB,CC,RR,UU,JMAX-2)
            DO 18 J=2,JMAX-1
              PSI(J,L)=-PSI(J,L)+TWO*R(N)*UU(J-1)
18          CONTINUE
19        CONTINUE
          DO 23 J=2,JMAX-1
            DO 21 L=2,JMAX-1
              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)
21          CONTINUE
            CALL TRIDAG(AA,BB,CC,RR,UU,JMAX-2)
            DO 22 L=2,JMAX-1
              U(J,L)=UU(L-1)
              PSI(J,L)=-PSI(J,L)+RFACT*UU(L-1)
22          CONTINUE
23        CONTINUE
24      CONTINUE
        ANORM=ZERO
        DO 26 J=2,JMAX-1
          DO 25 L=2,JMAX-1
            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)
25        CONTINUE
26      CONTINUE
        IF(ANORM.LT.EPS*ANORMG)RETURN
27    CONTINUE
      PAUSE 'MAXITS exceeded'
      END
      SUBROUTINE AMOEBA(P,Y,MP,NP,NDIM,FTOL,FUNK,ITER)
      PARAMETER (NMAX=20,ALPHA=1.0,BETA=0.5,GAMMA=2.0,ITMAX=500)
      DIMENSION P(MP,NP),Y(MP),PR(NMAX),PRR(NMAX),PBAR(NMAX)
      MPTS=NDIM+1
      ITER=0
1     ILO=1
      IF(Y(1).GT.Y(2))THEN
        IHI=1
        INHI=2
      ELSE
        IHI=2
        INHI=1
      ENDIF
      DO 11 I=1,MPTS
        IF(Y(I).LT.Y(ILO)) ILO=I
        IF(Y(I).GT.Y(IHI))THEN
          INHI=IHI
          IHI=I
        ELSE IF(Y(I).GT.Y(INHI))THEN
          IF(I.NE.IHI) INHI=I
        ENDIF
11    CONTINUE
      RTOL=2.*ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))
      IF(RTOL.LT.FTOL)RETURN
      IF(ITER.EQ.ITMAX) PAUSE 'Amoeba exceeding maximum iterations.'
      ITER=ITER+1
      DO 12 J=1,NDIM
        PBAR(J)=0.
12    CONTINUE
      DO 14 I=1,MPTS
        IF(I.NE.IHI)THEN
          DO 13 J=1,NDIM
            PBAR(J)=PBAR(J)+P(I,J)
13        CONTINUE
        ENDIF
14    CONTINUE
      DO 15 J=1,NDIM
        PBAR(J)=PBAR(J)/NDIM
        PR(J)=(1.+ALPHA)*PBAR(J)-ALPHA*P(IHI,J)
15    CONTINUE
      YPR=FUNK(PR)
      IF(YPR.LE.Y(ILO))THEN
        DO 16 J=1,NDIM
          PRR(J)=GAMMA*PR(J)+(1.-GAMMA)*PBAR(J)
16      CONTINUE
        YPRR=FUNK(PRR)
        IF(YPRR.LT.Y(ILO))THEN
          DO 17 J=1,NDIM
            P(IHI,J)=PRR(J)
17        CONTINUE
          Y(IHI)=YPRR
        ELSE
          DO 18 J=1,NDIM
            P(IHI,J)=PR(J)
18        CONTINUE
          Y(IHI)=YPR
        ENDIF
      ELSE IF(YPR.GE.Y(INHI))THEN
        IF(YPR.LT.Y(IHI))THEN
          DO 19 J=1,NDIM
            P(IHI,J)=PR(J)
19        CONTINUE
          Y(IHI)=YPR
        ENDIF
        DO 21 J=1,NDIM
          PRR(J)=BETA*P(IHI,J)+(1.-BETA)*PBAR(J)
21      CONTINUE
        YPRR=FUNK(PRR)
        IF(YPRR.LT.Y(IHI))THEN
          DO 22 J=1,NDIM
            P(IHI,J)=PRR(J)
22        CONTINUE
          Y(IHI)=YPRR
        ELSE
          DO 24 I=1,MPTS
            IF(I.NE.ILO)THEN
              DO 23 J=1,NDIM
                PR(J)=0.5*(P(I,J)+P(ILO,J))
                P(I,J)=PR(J)
23            CONTINUE
              Y(I)=FUNK(PR)
            ENDIF
24        CONTINUE
        ENDIF
      ELSE
        DO 25 J=1,NDIM
          P(IHI,J)=PR(J)
25      CONTINUE
        Y(IHI)=YPR
      ENDIF
      GO TO 1
      END
      SUBROUTINE ANNEAL(X,Y,IORDER,NCITY)
C Rev. 12/13/85
      DIMENSION X(NCITY),Y(NCITY),IORDER(NCITY),N(6)
      LOGICAL ANS
      ALEN(X1,X2,Y1,Y2)=SQRT((X2-X1)**2+(Y2-Y1)**2)
      NOVER=100*NCITY
      NLIMIT=10*NCITY
      TFACTR=0.9
      PATH=0.0
      T=0.5
      DO 11 I=1,NCITY-1
          I1=IORDER(I)
          I2=IORDER(I+1)
          PATH=PATH+ALEN(X(I1),X(I2),Y(I1),Y(I2))
11    CONTINUE
      I1=IORDER(NCITY)
      I2=IORDER(1)
      PATH=PATH+ALEN(X(I1),X(I2),Y(I1),Y(I2))
      IDUM=-1
      ISEED=111
      DO 15 J=1,100
          NSUCC=0
          DO 13 K=1,NOVER
12            N(1)=1+INT(NCITY*RAN3(IDUM))
              N(2)=1+INT((NCITY-1)*RAN3(IDUM))
              IF (N(2).GE.N(1)) N(2)=N(2)+1
              IDEC=IRBIT1(ISEED)
              NN=1+MOD((N(1)-N(2)+NCITY-1),NCITY)
              IF (NN.LT.3) GOTO 12
              IF (IDEC.EQ.0) THEN
                  N(3)=N(2)+INT(ABS(NN-2)*RAN3(IDUM))+1
                  N(3)=1+MOD(N(3)-1,NCITY)
                  CALL TRNCST(X,Y,IORDER,NCITY,N,DE)
                  CALL METROP(DE,T,ANS)
                  IF (ANS) THEN
                      NSUCC=NSUCC+1
                      PATH=PATH+DE
                      CALL TRNSPT(IORDER,NCITY,N)
                  ENDIF
              ELSE
                  CALL REVCST(X,Y,IORDER,NCITY,N,DE)
                  CALL METROP(DE,T,ANS)
                  IF (ANS) THEN
                      NSUCC=NSUCC+1
                      PATH=PATH+DE
                      CALL REVERS(IORDER,NCITY,N)
                  ENDIF
              ENDIF
              IF (NSUCC.GE.NLIMIT) GOTO 14
13        CONTINUE
14        WRITE(*,*)
          WRITE(*,'(1X,A,F10.6,A,F10.6)') 'T =',T,
     *            '   Path Length =',PATH
          WRITE(*,'(1X,A,I6)') 'Successful Moves: ',NSUCC
          T=T*TFACTR
          IF (NSUCC.EQ.0) RETURN
15    CONTINUE
      RETURN
      END
      SUBROUTINE AVEVAR(DATA,N,AVE,VAR)
      DIMENSION DATA(N)
      AVE=0.0
      VAR=0.0
      DO 11 J=1,N
        AVE=AVE+DATA(J)
11    CONTINUE
      AVE=AVE/N
      DO 12 J=1,N
        S=DATA(J)-AVE
        VAR=VAR+S*S
12    CONTINUE
      VAR=VAR/(N-1)
      RETURN
      END
