	PROGRAM D17R2
C	Driver for routine ADI
	IMPLICIT REAL*8(A-H,O-Z)
	PARAMETER(JMAX=11,PI=3.1415926)
	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)
	DO 12 I=1,JMAX
		DO 11 J=1,JMAX
			A(I,J)=-1.0
			B(I,J)=2.0
			C(I,J)=-1.0
			D(I,J)=-1.0
			E(I,J)=2.0
			F(I,J)=-1.0
			G(I,J)=0.0
			U(I,J)=0.0
11		CONTINUE
12	CONTINUE
	MID=JMAX/2+1
	G(MID,MID)=2.0
	ALPHA=2.0*(1.0-COS(PI/JMAX))
	BETA=2.0*(1.0-COS((JMAX-1)*PI/JMAX))
	ALIM=LOG(4.0*JMAX/PI)
	K=0
1	K=K+1
	IF (2**K .LT. ALIM) GOTO 1
	EPS=1.0E-4
	CALL ADI(A,B,C,D,E,F,G,U,JMAX,K,ALPHA,BETA,EPS)
	WRITE(*, '(1X,A)') 'ADI Solution:'
	DO 13 I=1,JMAX
		WRITE(*,'(1X,11F7.2)') (U(I,J),J=1,JMAX)
13	CONTINUE
	WRITE(*,'(/1X,A)') 'Test that solution satisfies Difference Eqns:'
	DO 15 I=2,JMAX-1
		DO 14 J=2,JMAX-1
			G(I,J)=-4.0*U(I,J)+U(I+1,J)
     *				+U(I-1,J)+U(I,J-1)+U(I,J+1)
14		CONTINUE
		WRITE(*,'(8X,9F7.2)') (G(I,J),J=2,JMAX-1)
15	CONTINUE
	END

	SUBROUTINE TRIDAG(A,B,C,R,U,N)
C	This is a double precision version for use with ADI
	IMPLICIT REAL*8(A-H,O-Z)
	PARAMETER (NMAX=100)
	DIMENSION GAM(NMAX),A(N),B(N),C(N),R(N),U(N)
	IF(B(1).EQ.0.)PAUSE
	BET=B(1)
	U(1)=R(1)/BET
	DO 11 J=2,N
		GAM(J)=C(J-1)/BET
		BET=B(J)-A(J)*GAM(J)
		IF(BET.EQ.0.)PAUSE
		U(J)=(R(J)-A(J)*U(J-1))/BET
11	CONTINUE
	DO 12 J=N-1,1,-1
		U(J)=U(J)-GAM(J+1)*U(J+1)
12	CONTINUE
	RETURN
	END
	PROGRAM D10R5
C	Driver for routine AMOEBA
	EXTERNAL FAMOEB
	PARAMETER(NP=3,MP=4,FTOL=1.0E-6)
	DIMENSION P(MP,NP),X(NP),Y(MP)
	DATA P/0.0,1.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,1.0/
	NDIM=NP
	DO 12 I=1,MP
		DO 11 J=1,NP
			X(J)=P(I,J)
11		CONTINUE
		Y(I)=FAMOEB(X)
12	CONTINUE
	CALL AMOEBA(P,Y,MP,NP,NDIM,FTOL,FAMOEB,ITER)
	WRITE(*,'(/1X,A,I3)') 'Iterations: ',ITER
	WRITE(*,'(/1X,A)') 'Vertices of final 3-D simplex and'
	WRITE(*,'(1X,A)') 'function values at the vertices:'
	WRITE(*,'(/3X,A,T11,A,T23,A,T35,A,T45,A/)') 'I',
     *		'X(I)','Y(I)','Z(I)','FUNCTION'
	DO 13 I=1,MP
		WRITE(*,'(1X,I3,4F12.6)') I,(P(I,J),J=1,NP),Y(I)
13	CONTINUE
	WRITE(*,'(/1X,A)') 'True minimum is at (0.5,0.6,0.7)'
	END

	FUNCTION FAMOEB(X)
	DIMENSION X(3)
	FAMOEB=0.6-BESSJ0((X(1)-0.5)**2+(X(2)-0.6)**2+(X(3)-0.7)**2)
	END
	PROGRAM D13R5
C	Driver for routine AVEVAR
	PARAMETER(NPTS=1000, EPS=0.1)
	DIMENSION DATA(NPTS)
C	Generate Gaussian distributed data
	IDUM=-5
	WRITE(*,'(1X,T4,A,T14,A,T26,A)') 'Shift','Average','Variance'
	DO 12 I=1,11
		SHIFT=(I-1)*EPS
		DO 11 J=1,NPTS
			DATA(J)=SHIFT+I*GASDEV(IDUM)
11		CONTINUE
		CALL AVEVAR(DATA,NPTS,AVE,VAR)
		WRITE(*,'(1X,F6.2,2F12.2)') SHIFT,AVE,VAR
12	CONTINUE
	END
