	PROGRAM D2R7
C	Driver for routine TOEPLZ
	PARAMETER(N=5,N2=2*N)
	DIMENSION X(N),Y(N),R(N2)
	DO 11 I=1,N
		Y(I)=0.1*I
11	CONTINUE
	DO 12 I=1,2*N-1
		R(I)=0.1*I
12	CONTINUE
	CALL TOEPLZ(R,X,Y,N)
	WRITE(*,*) 'Solution vector:'
	DO 13 I=1,N
		WRITE(*,'(5X,A2,I1,A4,E13.6)') 'X(',I,') = ',X(I)
13	CONTINUE
	WRITE(*,'(/1X,A)') 'Test of solution:'
	WRITE(*,'(1X,T6,A,T19,A)') 'mtrx*soln','original'
	DO 15 I=1,N
		SUM=0.0
		DO 14 J=1,N
			SUM=SUM+R(N+I-J)*X(J)
14		CONTINUE
		WRITE(*,'(1X,2F12.4)') SUM,Y(I)
15	CONTINUE
	END
	PROGRAM D13R6B
C	Driver for routine TPTEST
C	Compare two correlated distributions vs. two 
C	uncorrelated distributions
	PARAMETER(NPTS=500,EPS=0.01,NSHFT=10,ANOISE=0.3)
	DIMENSION DATA1(NPTS),DATA2(NPTS),DATA3(NPTS)
	IDUM=-5
	WRITE(*,'(1X,T18,A,T46,A)') 'Correlated:','Uncorrelated:'
	WRITE(*,'(1X,T4,A,T18,A,T25,A,T46,A,T53,A)') 
     *		'Shift','T','Probability','T','Probability'
	OFFSET=(NSHFT/2)*EPS
	DO 11 J=1,NPTS
		GAUSS=GASDEV(IDUM)
		DATA1(J)=GAUSS
		DATA2(J)=GAUSS+ANOISE*GASDEV(IDUM)
		DATA3(J)=GASDEV(IDUM)+ANOISE*GASDEV(IDUM)
11	CONTINUE
	CALL AVEVAR(DATA1,NPTS,AVE1,VAR1)
	CALL AVEVAR(DATA2,NPTS,AVE2,VAR2)
	CALL AVEVAR(DATA3,NPTS,AVE3,VAR3)
	DO 12 J=1,NPTS
		DATA1(J)=DATA1(J)-AVE1+OFFSET
		DATA2(J)=DATA2(J)-AVE2
		DATA3(J)=DATA3(J)-AVE3		
12	CONTINUE
	DO 14 I=1,NSHFT+1
		SHIFT=I*EPS
		DO 13 J=1,NPTS
			DATA2(J)=DATA2(J)+EPS
			DATA3(J)=DATA3(J)+EPS
13		CONTINUE
		CALL TPTEST(DATA1,DATA2,NPTS,T1,PROB1)
		CALL TPTEST(DATA1,DATA3,NPTS,T2,PROB2)
		WRITE(*,'(1X,F6.2,2X,2F12.4,4X,2F12.4)') 
     *			SHIFT,T1,PROB1,T2,PROB2
14	CONTINUE
	END
	PROGRAM D11R4
C	Driver for routine TQLI
	PARAMETER(NP=10,TINY=1.0E-6)
	DIMENSION A(NP,NP),C(NP,NP),D(NP),E(NP),F(NP)
	DATA C/5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,-4.0,
     *		4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,
     *		3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,
     *		2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,
     *		1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,
     *		0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,
     *		-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,
     *		-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,
     *		-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,
     *		-4.0,-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0/
	DO 12 I=1,NP
		DO 11 J=1,NP
			A(I,J)=C(I,J)
11		CONTINUE
12	CONTINUE
	CALL TRED2(A,NP,NP,D,E)
	CALL TQLI(D,E,NP,NP,A)
	WRITE(*,'(/1X,A)') 'Eigenvectors for a real symmetric matrix'
	DO 16 I=1,NP
		DO 14 J=1,NP
			F(J)=0.0
			DO 13 K=1,NP
				F(J)=F(J)+C(J,K)*A(K,I)
13			CONTINUE
14		CONTINUE
		WRITE(*,'(/1X,A,I3,A,F10.6)') 'Eigenvalue',I,' =',D(I)
		WRITE(*,'(/1X,T7,A,T17,A,T31,A)') 'Vector','Mtrx*Vect.','Ratio'
		DO 15 J=1,NP
			IF (ABS(A(J,I)).LT.TINY) THEN
				WRITE(*,'(1X,2F12.6,A12)') A(J,I),F(J),'div. by 0'
			ELSE
				WRITE(*,'(1X,2F12.6,E14.6)') A(J,I),F(J),
     *					F(J)/A(J,I)
			ENDIF
15		CONTINUE
		WRITE(*,'(/1X,A)') 'press ENTER to continue...'
		READ(*,*)
16	CONTINUE
	END
	PROGRAM D4R1
C	Driver for routine TRAPZD
	PARAMETER(NMAX=14, PIO2=1.5707963)
	EXTERNAL FUNC,FINT
	A=0.0
	B=PIO2
	WRITE(*,'(1X,A)') 'Integral of FUNC with 2^(n-1) points'
	WRITE(*,'(1X,A,F10.6)') 'Actual value of integral is',FINT(B)-FINT(A)
	WRITE(*,'(1X,T7,A,T16,A)') 'n','Approx. Integral'
	DO 11 I=1,NMAX
		CALL TRAPZD(FUNC,A,B,S,I)
		WRITE(*,'(1X,I6,F20.6)') I,S
11	CONTINUE
	END

	FUNCTION FUNC(X)
	FUNC=(X**2)*(X**2-2.0)*SIN(X)
	END

	FUNCTION FINT(X)
C	Integral of FUNC
	FINT=4.0*X*((X**2)-7.0)*SIN(X)
     *		-((X**4)-14.0*(X**2)+28.0)*COS(X)
	END
	PROGRAM D11R3
C	Driver for routine TRED2
	PARAMETER(NP=10)
	DIMENSION A(NP,NP),C(NP,NP),D(NP),E(NP),F(NP,NP)
	DATA C/5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,-4.0,
     *		4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,-3.0,
     *		3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,-2.0,
     *		2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,-1.0,
     *		1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,0.0,
     *		0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,1.0,
     *		-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,2.0,
     *		-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,3.0,
     *		-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0,4.0,
     *		-4.0,-3.0,-2.0,-1.0,0.0,1.0,2.0,3.0,4.0,5.0/
	DO 12 I=1,NP
		DO 11 J=1,NP
			A(I,J)=C(I,J)
11		CONTINUE
12	CONTINUE
	CALL TRED2(A,NP,NP,D,E)
	WRITE(*,'(/1X,A)') 'Diagonal elements'
	WRITE(*,'(1X,5F12.6)') (D(I),I=1,NP)
	WRITE(*,'(/1X,A)') 'Off-diagonal elements'
	WRITE(*,'(1X,5F12.6)') (E(I),I=2,NP)
C	Check transformation matrix
	DO 16 J=1,NP
		DO 15 K=1,NP
		F(J,K)=0.0
			DO 14 L=1,NP
				DO 13 M=1,NP
					F(J,K)=F(J,K)
     *					+A(L,J)*C(L,M)*A(M,K)
13				CONTINUE
14			CONTINUE
15		CONTINUE
16	CONTINUE
C	How does it look?
	WRITE(*,'(/1X,A)') 'Tridiagonal matrix'
	DO 17 I=1,NP
		WRITE(*,'(1X,10F7.2)') (F(I,J),J=1,NP)
17	CONTINUE
	END
	PROGRAM D2R4
C	Driver for routine TRIDAG
	PARAMETER (NP=20)
	DIMENSION DIAG(NP),SUPERD(NP),SUBD(NP),RHS(NP),U(NP)
	CHARACTER TXT*3
	OPEN(UNIT=5,FILE='MATRX2.DAT',STATUS='OLD')
10	READ(5,'(A3)') TXT
	IF (TXT.EQ.'END') GOTO 99
	READ(5,*)
	READ(5,*) N
	READ(5,*)
	READ(5,*) (DIAG(K), K=1,N)
	READ(5,*)
	READ(5,*) (SUPERD(K), K=1,N-1)
	READ(5,*)
	READ(5,*) (SUBD(K), K=2,N)
	READ(5,*)
	READ(5,*) (RHS(K), K=1,N)
C	Carry out solution
	CALL TRIDAG(SUBD,DIAG,SUPERD,RHS,U,N)
	WRITE(*,*) 'The solution vector is:'
	WRITE(*,'(1X,6F12.6)') (U(K), K=1,N)
C	Test solution
	WRITE(*,*) '(matrix)*(sol''n vector) should be:'
	WRITE(*,'(1X,6F12.6)') (RHS(K), K=1,N)
	WRITE(*,*) 'Actual result is:'
	DO 11 K=1,N
		IF (K.EQ.1) THEN
			RHS(K)=DIAG(1)*U(1) + SUPERD(1)*U(2)
		ELSE IF (K.EQ.N) THEN
			RHS(K)=SUBD(N)*U(N-1) + DIAG(N)*U(N)
		ELSE
			RHS(K)=SUBD(K)*U(K-1) + DIAG(K)*U(K)
     *				+ SUPERD(K)*U(K+1)
		ENDIF
11	CONTINUE
	WRITE(*,'(1X,6F12.6)') (RHS(K), K=1,N)
	WRITE(*,*) '***********************************'
	WRITE(*,*) 'Press RETURN for next problem:'
	READ(*,*)
	GOTO 10
99	CLOSE(5)
	END
	PROGRAM D13R4
C	Driver for routine TTEST
	PARAMETER(NPTS=1024, MPTS=512, EPS=0.02,NSHFT=10)
	DIMENSION DATA1(NPTS),DATA2(MPTS)
C	Generate Gaussian distributed data
	IDUM=-5
	DO 11 I=1,NPTS
		DATA1(I)=GASDEV(IDUM)
11	CONTINUE
	IDUM=-11
	DO 12 I=1,MPTS
		DATA2(I)=(NSHFT/2.0)*EPS+GASDEV(IDUM)
12	CONTINUE
	WRITE(*,'(/1X,T4,A,T18,A,T25,A)') 'Shift','T','Probability'
	DO 14 I=1,NSHFT+1
		CALL TTEST(DATA1,NPTS,DATA2,MPTS,T,PROB)
		SHIFT=(I-1)*EPS
		WRITE(*,'(1X,F6.2,2F12.2)') SHIFT,T,PROB
		DO 13 J=1,NPTS
			DATA1(J)=DATA1(J)+EPS
13		CONTINUE
14	CONTINUE
	END
	PROGRAM D13R6A
C	Driver for routine TUTEST
	PARAMETER(NPTS=5000,MPTS=1000,EPS=0.02,VAR1=1.0,
     *		VAR2=4.0,NSHFT=10)
	DIMENSION DATA1(NPTS),DATA2(MPTS)
C	Generate two Gaussian distributions of different variance
	IDUM=-51773
	FCTR1=SQRT(VAR1)
	DO 11 I=1,NPTS
		DATA1(I)=FCTR1*GASDEV(IDUM)
11	CONTINUE
	FCTR2=SQRT(VAR2)
	DO 12 I=1,MPTS
		DATA2(I)=(NSHFT/2.0)*EPS+FCTR2*GASDEV(IDUM)
12	CONTINUE
	WRITE(*,'(1X,A,F6.2)') 'Distribution #1 : variance = ',VAR1
	WRITE(*,'(1X,A,F6.2/)') 'Distribution #2 : variance = ',VAR2
	WRITE(*,'(1X,T4,A,T18,A,T25,A)') 'Shift','T','Probability'
	DO 14 I=1,NSHFT+1
		CALL TUTEST(DATA1,NPTS,DATA2,MPTS,T,PROB)
		SHIFT=(I-1)*EPS
		WRITE(*,'(1X,F6.2,2F12.2)') SHIFT,T,PROB
		DO 13 J=1,NPTS
			DATA1(J)=DATA1(J)+EPS
13		CONTINUE
14	CONTINUE
	END
	PROGRAM D12R2
C	Driver for routine TWOFFT
	PARAMETER(N=32,N2=2*N,PER=8.0,PI=3.14159)
	DIMENSION DATA1(N),DATA2(N),FFT1(N2),FFT2(N2)
	DO 11 I=1,N
		X=2.0*PI*I/PER
		DATA1(I)=NINT(COS(X))
		DATA2(I)=NINT(SIN(X))
11	CONTINUE
	CALL TWOFFT(DATA1,DATA2,FFT1,FFT2,N)
	WRITE(*,*) 'Fourier transform of first function:'
	CALL PRNTFT(FFT1,N2)
	WRITE(*,*) 'Fourier transform of second function:'
	CALL PRNTFT(FFT2,N2)
C	Invert transform
	ISIGN=-1
	CALL FOUR1(FFT1,N,ISIGN)
	WRITE(*,*) 'Inverted transform = first function:'
	CALL PRNTFT(FFT1,N2)
	CALL FOUR1(FFT2,N,ISIGN)
	WRITE(*,*) 'Inverted transform = second function:'
	CALL PRNTFT(FFT2,N2)
	END

	SUBROUTINE PRNTFT(DATA,N2)
	DIMENSION DATA(N2)
	WRITE(*,'(1X,T7,A,T13,A,T24,A,T35,A,T47,A)') 
     *		'n','Real(n)','Imag.(n)','Real(N-n)','Imag.(N-n)'
	WRITE(*,'(1X,I6,4F12.6)') 0,DATA(1),DATA(2),DATA(1),DATA(2)
	DO 11 I=3,(N2/2)+1,2
		M=(I-1)/2
		NN2=N2+2-I
		WRITE(*,'(1X,I6,4F12.6)') M,DATA(I),DATA(I+1),
     *			DATA(NN2),DATA(NN2+1)
11	CONTINUE
	WRITE(*,'(/1X,A)') ' press RETURN to continue ...'
	READ(*,*)
	RETURN
	END
