c subroutine GAUSSIAN, new common parameter sigma added 3.12.98 N.B.

      INCLUDE 'cont1.h'

      INTEGER i,iq,nq,nr,iw,st,iqx,iqy,iqz,steps
      REAL*8    a(-wmax:wmax),xt1(tmax),x2(20),ss(3,20)
      REAL*8    s0,s1,s2,alpha0,dw,CHI2,RAN
      CHARACTER*8 string

      CALL READIN(nq)
      WRITE(*,*)'Number of frequencies'
      READ(*,*)nw
      WRITE(*,*)'Delta frequency'
      READ(*,*)dw
      WRITE(*,*)'Number of annealing steps'
      READ(*,*)steps
      WRITE(*,*)'Starting Alpha'
      READ(*,*)alpha0
      WRITE(*,*)'Number of smoothed runs'
      READ(*,*)nr
      seed=53718551
      WRITE(*,*)'Random number seed'
      READ(*,*)seed

      DO 10 iw=-nw,nw
         w(iw)=iw*dw
 10   CONTINUE

      DO 50 iq=1,nq
         WRITE(*,*)'Starting Maxent ',iq         
         CALL GETDATA(iq,iqx,iqy,iqz)
         CALL CONTINUE(steps,alpha0,nr,a)
         OPEN(UNIT=10,FILE='aw.dat',STATUS='unknown')
         DO 20 i=1,100000
            READ(10,*,END=30)string
 20      CONTINUE
c 30      CONTINUE
 30      WRITE(10,*)iqx,iqy,iqz
         DO 40 iw=-nw,nw
            WRITE(10,80)w(iw),a(iw)
 40      CONTINUE
         CALL CALCXT(a,xt1)
         x2(iq)=CHI2(xt1)
         CALL SUMRULES(a,s0,s1,s2)
         ss(1,iq)=s0
         ss(2,iq)=s1
         ss(3,iq)=s2
 50   CONTINUE
 80   FORMAT(F9.5,'   ',F9.5)

      OPEN(UNIT=20,FILE='sums.dat',STATUS='new')
      DO 200 i=1,nq
         WRITE(20,310)i,ss(1,i),ss(2,i),ss(3,i)
 200  CONTINUE
      CLOSE(10)
 310  FORMAT(I4,'   Sums:  ',3F9.5)

      END

c------------------------------------------------------------------------------
      SUBROUTINE CONTINUE(steps,alpha0,nr,a)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER i,iq,nr,st,steps
      REAL*8    a(-wmax:wmax),mm(-wmax:wmax),xt1(tmax)
      REAL*8    aa(-wmax:wmax,5),x2(5),ss(3,5)
      REAL*8    n0,rfac,temp,s0,s1,s2,alpha,alpha0,ent,tr,ls,adiv
      REAL*8    CHI2,ENTROPY,RAN

      n0=1.
      CALL INITKER
      CALL INITF0
      CALL INITDLDA
c  NB generalized !!!!!!!!
      IF (sigma.lt.f0(0)) THEN
         CALL FLATMODEL(mm)
      ELSE
         CALL GAUSSIAN(mm)
      ENDIF
      CALL NORMALIZE(n0,mm)
      DO 10 i=-nw,nw
         a(i)=RAN(seed)
 10   CONTINUE
      CALL NORMALIZE(n0,a)
      temp=10.
      rfac=1.

      alpha=alpha0
 20   CALL MAXENT(mm,a,rfac,temp,alpha,steps)
      ent=ENTROPY(a,mm)
      CALL LAMBDA(a,alpha,tr)
      ls=-2.*ent*alpha
      WRITE(*,*)'-2*alpha*S = ',ls
      WRITE(*,*)'Trace      = ',tr
      temp=0.001
      rfac=0.05
      IF (ABS(ls/tr-1.).LT.0.05) THEN
          GOTO 30
       ELSE
          IF (tr/ls.LT.0.05) THEN
             alpha=alpha*0.05
          ELSE
            alpha=alpha*(tr/ls)*(1.+0.001*(RAN(seed)-0.5))
         ENDIF
         GOTO 20
      ENDIF

 30   DO 100 st=1,nr
         WRITE(*,*)'Starting smoothing ',st         
         CALL SMOOTH(a,3)
         CALL NORMALIZE(n0,a)
         temp=0.005
         rfac=0.005
         CALL MAXENT(mm,a,rfac,temp,alpha,steps)
 100  CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE SMOOTH(a,ns)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER i,j,i1,i2,ii,ns
      REAL*8    a(-wmax:wmax),a1(-wmax:wmax)
      
      DO 20 i=-nw,nw
         i1=i-ns
         IF (i1.LT.-nw) i1=-nw
         i2=i+ns
         IF (i2.GT.nw) i2=nw
         a1(i)=0.
         ii=0
         DO 10 j=i1,i2
            ii=ii+1
            a1(i)=a1(i)+a(j)
 10      CONTINUE
         a1(i)=a1(i)/FLOAT(ii)
 20   CONTINUE
      DO 30 i=-nw,nw
         a(i)=a1(i)
 30   CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------
       SUBROUTINE READIN(nq)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER xyz(3,20)
      REAL*8  axt1(tmax,20),sxt1(tmax,20)
      COMMON/BQDATA/axt1,sxt1,xyz

      INTEGER it,iit,iq,nq,nt1

      WRITE(*,*)'Number of momenta'
      READ(*,*)nq
      WRITE(*,*)'Number of times'
      READ(*,*)nt
ctttttttttt M.U. , 8/15/95
      WRITE(*,*)'Number of times in the MC sampling'
      READ(*,*)nt1

      WRITE(*,*)'Beta = , sigma = '
      READ(*,*)beta,sigma
      OPEN(UNIT=10,FILE='gav.dat',STATUS='old')
      DO 20 iq=1,nq
         READ(10,*)xyz(1,iq),xyz(2,iq),xyz(3,iq)
         DO 10 it=1,nt
            READ(10,*)iit,axt1(it,iq),sxt1(it,iq)
            axt1(it,iq)=-axt1(it,iq)
            t(it)=FLOAT(iit)*beta/FLOAT(nt1-1)
            IF (ABS(sxt1(it,iq)).LT.0.00001) sxt1(it,iq)=0.00001
            sxt1(it,iq)=1./sxt1(it,iq)**2
 10      CONTINUE
         WRITE(*,*)'Read momentum ',xyz(1,iq),xyz(2,iq),xyz(3,iq)
 20   CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE GETDATA(iq,iqx,iqy,iqz)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER xyz(3,20)
      REAL*8    axt1(tmax,20),sxt1(tmax,20)
      COMMON/BQDATA/axt1,sxt1,xyz

      INTEGER i,iq,iqx,iqy,iqz
      
      iqx=xyz(1,iq)
      iqy=xyz(2,iq)
      iqz=xyz(3,iq)
      DO 10 i=1,nt
         axt(i)=axt1(i,iq)
         sxt(i)=sxt1(i,iq)
 10   CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE MAXENT(mm,a,rfac,temp,alpha,steps)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER i,j,j1,j2,jj,it,steps,acc,try
      REAL*8    a(-wmax:wmax),da(-wmax:wmax),xt1(tmax),xt2(tmax)
      REAL*8    dj1,dj2,x1,x2,temp,alpha,wgt,eps,RAN
      REAL*8    CHI2,aa1,aa2,de,p,mm1,mm2,mm(-wmax:wmax),rfac,arat,am

      am=a(-nw)
      DO 10 i=-nw,nw         
         da(i)=rfac*(a(i)/10.+0.01)
         IF (a(i).GT.am) am=a(i)
 10   CONTINUE
      eps=1.E-12

      DO 200 i=1,steps
         CALL CALCXT(a,xt1)
         x1=CHI2(xt1)
         DO 100 jj=1,2*nw+1

 35         j1=MIN(INT(RAN(seed)*(2*nw+1)),2*nw)-nw
 40         dj1=da(j1)*(RAN(seed)-0.5)
            IF (a(j1)+dj1.LT.0.) GOTO 40
 50         j2=MIN(INT(RAN(seed)*(2*nw+1)),2*nw)-nw
            IF (j1.EQ.j2) GOTO 50
            dj2=-dj1+0.05*da(j2)*(RAN(seed)-0.5)            
            IF (a(j2)+dj2.LT.0.) GOTO 35
            IF (a(j1).GT.0.1*am) try=try+1

            DO 70 it=1,nt
               xt2(it)=xt1(it)+dj1*ker(j1,it)+dj2*ker(j2,it)
 70         CONTINUE
            x2=CHI2(xt2)
            aa1=a(j1)+dj1
            aa2=a(j2)+dj2
            mm1=mm(j1)
            mm2=mm(j2)
            de=0.
            IF (aa1.GT.eps) de=de-f0(j1)*aa1*LOG(aa1/mm1)
            IF (a(j1).GT.eps) de=de+f0(j1)*a(j1)*LOG(a(j1)/mm1)
            IF (aa2.GT.eps) de=de-f0(j2)*aa2*LOG(aa2/mm2)
            IF (a(j2).GT.eps) de=de+f0(j2)*a(j2)*LOG(a(j2)/mm2)            
            wgt=((x1-x2)/2.+alpha*de)/temp
            IF (wgt.GT.0) THEN
               p=1.
            ELSEIF (wgt.LT.10.) THEN               
               p=EXP(wgt)
            ELSE
               p=0.
            ENDIF
            IF (RAN(seed).LT.p) THEN
               IF (a(j1).GT.0.1*am) acc=acc+1
               a(j1)=aa1
               a(j2)=aa2
               DO 80 it=1,nt
                  xt1(it)=xt2(it)
 80            CONTINUE
               x1=x2
             ENDIF

 100      CONTINUE

          IF (MOD(i,100).EQ.0) THEN
             arat=FLOAT(acc)/FLOAT(try)
             IF (arat.GT.0.1) THEN
                IF (rfac.LT.0.01) rfac=rfac*1.5
             ELSE
                IF (rfac.GT.0.001) rfac=rfac/1.5
             ENDIF
             am=a(1)
             DO 110 j=-nw,nw
                da(j)=rfac*(0.1*a(j)+0.01*am)
                IF (a(j).GT.am) am=a(j)

 110         CONTINUE
             IF (MOD(i,steps/10).EQ.0) THEN
                WRITE(*,*)i,x1,acc,rfac
             ENDIF
             temp=temp/1.5
             acc=0
             try=0
          ENDIF
 200   CONTINUE
 300   FORMAT(I8,F12.6,I8,F12.6)
       RETURN
       END

c------------------------------------------------------------------------------
      SUBROUTINE NORMALIZE(n0,a)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER iw
      REAL*8    a(-wmax:wmax),n0,n1,f

      n1=0.
      DO 10 iw=-nw,nw
         n1=n1+a(iw)*f0(iw)
 10   CONTINUE
      f=n0/n1
      DO 20 iw=-nw,nw
         a(iw)=f*a(iw)
 20   CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE CALCXT(a,xt)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER it,iw
      REAL*8    a(-wmax:wmax),xt(tmax)

      DO 20 it=1,nt
         xt(it)=0.
         DO 10 iw=-nw,nw
            xt(it)=xt(it)+a(iw)*ker(iw,it)
 10      CONTINUE
 20   CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------
      REAL*8 FUNCTION CHI2(xt)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER i
      REAL*8    x2,xt(tmax)

      x2=0.
      DO 10 i=1,nt
          x2=x2+sxt(i)*(axt(i)-xt(i))**2
 10   CONTINUE
      CHI2=x2
      RETURN
      END

c------------------------------------------------------------------------------
      REAL*8 FUNCTION ENTROPY(a,mm)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER i
      REAL*8    a(-wmax:wmax),mm(-wmax:wmax),ent,eps

      ent=0.
      eps=1.E-12
      DO 10 i=-nw,nw
         IF (a(i).GT.eps.AND.mm(i).GT.eps) THEN
            ent=ent-a(i)*LOG(a(i)/mm(i))*f0(i)
         ENDIF
 10   CONTINUE
      ENTROPY=ent
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE LAMBDA(a,alpha,tr)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      REAL*8 dlda(2*wmax+1,2*wmax+1)
      COMMON/BLAM/dlda

      INTEGER i,j,iw,jw,indx(2*wmax+1)
      REAL*8    lam(2*wmax+1,2*wmax+1),alam(2*wmax+1,2*wmax+1)
      REAL*8    alami(2*wmax+1,2*wmax+1),a(-wmax:wmax)
      REAL*8    alpha,tr,dw

      dw=w(2)-w(1)
      DO 20 j=1,2*nw+1
         DO 10 i=1,2*nw+1
            iw=i-nw-1
            jw=j-nw-1
            lam(i,j)=SQRT(a(iw))*dlda(i,j)*SQRT(a(jw))*dw
            alam(i,j)=lam(i,j)
 10      CONTINUE
         alam(j,j)=lam(j,j)+alpha
 20   CONTINUE
      CALL INVERT(alam,alami,indx,2*nw+1,2*wmax+1)      
      tr=0.
      DO 40 i=1,2*nw+1
         DO 30 j=1,2*nw+1
            tr=tr+lam(i,j)*alami(j,i)
 30      CONTINUE
 40   CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE INITKER
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER it,iw
      REAL*8    bw,tw,dw

      dw=w(2)-w(1)
      DO 20 iw=-nw,nw
         DO 10 it=1,nt
            tw=EXP(-t(it)*w(iw))
            bw=1.+EXP(-beta*w(iw))
            ker(iw,it)=dw*tw/bw
 10      CONTINUE
 20   CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE INITDLDA
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      REAL*8 dlda(2*wmax+1,2*wmax+1)
      COMMON/BLAM/dlda

      INTEGER i,j,iw,jw,it
      REAL*8    dw
      
      dw=w(2)-w(1)
      DO 20 j=1,2*nw+1
      DO 20 i=1,2*nw+1
         iw=i-nw-1
         jw=j-nw-1
         dlda(i,j)=0.
         DO 10 it=1,nt
            dlda(i,j)=dlda(i,j)+ker(iw,it)*sxt(it)*ker(jw,it)/dw**2       
 10      CONTINUE
 20   CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE INITF0
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER iw
      REAL*8    dw

      dw=w(2)-w(1)
      DO 10 iw=-nw,nw
         f0(iw)=dw
         f1(iw)=dw*w(iw)
         f2(iw)=dw*w(iw)**2
 10   CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE FLATMODEL(mm)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER i
      REAL*8    mm(-wmax:wmax)

      DO 10 i=-nw,nw
         mm(i)=1.
 10   CONTINUE
      RETURN
      END

c new N.B. 3.12.98 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

c------------------------------------------------------------------------------
      SUBROUTINE GAUSSIAN(mm)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER i
      REAL*8    mm(-wmax:wmax), f

      f = 0.5/(sigma*sigma)
      DO 10 i=-nw,nw
         mm(i)=exp(-f*w(i)*w(i))
 10   CONTINUE
      RETURN
      END



c------------------------------------------------------------------------------
      SUBROUTINE SUMRULES(a,s0,s1,s2)
c------------------------------------------------------------------------------

      INCLUDE 'cont1.h'

      INTEGER i
      REAL*8    a(-wmax:wmax),s0,s1,s2

      s0=0.
      s1=0.
      s2=0.
      DO 10 i=-nw,nw
         s0=s0+a(i)*f0(i)
         s1=s1+a(i)*f1(i)
         s2=s2+a(i)*f2(i)
 10   CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------
      SUBROUTINE INVERT(a,y,indx,n,np)
c------------------------------------------------------------------------------
      REAL*8 a(np,np),y(np,np)
      DIMENSION indx(np)
      DO 12 i=1,n
         DO 11 j=1,n
            y(i,j)=0.
 11      CONTINUE
         y(i,i)=1.
 12   CONTINUE
      CALL LUDCMP(a,n,np,indx,d)
      DO 13 j=1,n
         CALL LUBKSB(a,n,np,indx,y(1,j))
 13   CONTINUE
      RETURN
      END
c------------------------------------------------------------------------------
      SUBROUTINE LUDCMP(A,N,NP,INDX,D)
c------------------------------------------------------------------------------
      PARAMETER (NMAX=400,TINY=1.0E-20)
      REAL*8 A(NP,NP),VV(NMAX)
      DIMENSION INDX(N)
      D=1.
      DO 12 I=1,N
        AAMAX=0.
        DO 11 J=1,N
          IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
11      CONTINUE
        IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'
        VV(I)=1./AAMAX
12    CONTINUE
      DO 19 J=1,N
        IF (J.GT.1) THEN
          DO 14 I=1,J-1
            SUM=A(I,J)
            IF (I.GT.1)THEN
              DO 13 K=1,I-1
                SUM=SUM-A(I,K)*A(K,J)
13            CONTINUE
              A(I,J)=SUM
            ENDIF
14        CONTINUE
        ENDIF
        AAMAX=0.
        DO 16 I=J,N
          SUM=A(I,J)
          IF (J.GT.1)THEN
            DO 15 K=1,J-1
              SUM=SUM-A(I,K)*A(K,J)
15          CONTINUE
            A(I,J)=SUM
          ENDIF
          DUM=VV(I)*ABS(SUM)
          IF (DUM.GE.AAMAX) THEN
            IMAX=I
            AAMAX=DUM
          ENDIF
16      CONTINUE
        IF (J.NE.IMAX)THEN
          DO 17 K=1,N
            DUM=A(IMAX,K)
            A(IMAX,K)=A(J,K)
            A(J,K)=DUM
17        CONTINUE
          D=-D
          VV(IMAX)=VV(J)
        ENDIF
        INDX(J)=IMAX
        IF(J.NE.N)THEN
          IF(A(J,J).EQ.0.)A(J,J)=TINY
          DUM=1./A(J,J)
          DO 18 I=J+1,N
            A(I,J)=A(I,J)*DUM
18        CONTINUE
        ENDIF
19    CONTINUE
      IF(A(N,N).EQ.0.)A(N,N)=TINY
      RETURN
      END
c------------------------------------------------------------------------------
      SUBROUTINE LUBKSB(A,N,NP,INDX,B)
c------------------------------------------------------------------------------
      REAL*8 A(NP,NP),B(N)
      DIMENSION INDX(N)
      II=0
      DO 12 I=1,N
        LL=INDX(I)
        SUM=B(LL)
        B(LL)=B(I)
        IF (II.NE.0)THEN
          DO 11 J=II,I-1
            SUM=SUM-A(I,J)*B(J)
11        CONTINUE
        ELSE IF (SUM.NE.0.) THEN
          II=I
        ENDIF
        B(I)=SUM
12    CONTINUE
      DO 14 I=N,1,-1
        SUM=B(I)
        IF(I.LT.N)THEN
          DO 13 J=I+1,N
            SUM=SUM-A(I,J)*B(J)
13        CONTINUE
        ENDIF
        B(I)=SUM/A(I,I)
14    CONTINUE
      RETURN
      END

c------------------------------------------------------------------------------


c       USE THESE COMMENTED OUT LINES IF REAL*8 DESIRED.
      REAL*8 FUNCTION RAN(IDUM)
      IMPLICIT REAL*8(A-H,O-Z)
c      FUNCTION RAN(IDUM)
      save
      PARAMETER (M=714025,IA=1366,IC=150889,RM=1.4005112E-6)
      DIMENSION IR(97)
      DATA IFF /0/
      IF(IDUM.LT.0.OR.IFF.EQ.0)THEN
        IFF=1
        IDUM=MOD(IC-IDUM,M)
        DO 11 J=1,97
          IDUM=MOD(IA*IDUM+IC,M)
          IR(J)=IDUM
11      CONTINUE
        IDUM=MOD(IA*IDUM+IC,M)
        IY=IDUM
      ENDIF
      J=1+(97*IY)/M
      IF(J.GT.97.OR.J.LT.1)PAUSE
      IY=IR(J)
      RAN=IY*RM
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
      RETURN
      END





