	PROGRAM D9R10
C	Driver for routine LAGUER
	PARAMETER(M=4,MP1=M+1,NTRY=21,EPS=1.0E-6)
	COMPLEX A(MP1),Y(NTRY),X
	LOGICAL POLISH
	DATA A/(0.0,2.0),(0.0,0.0),(-1.0,-2.0),(0.0,0.0),(1.0,0.0)/
	WRITE(*,'(/1X,A)') 'Roots of polynomial x^4-(1+2i)*x^2+2i'
	WRITE(*,'(/1X,T16,A,T29,A/)') 'Real','Complex'
	N=0
	POLISH=.FALSE.
	DO 12 I=1,NTRY
		X=CMPLX((I-11.0)/10.0,(I-11.0)/10.0)
		CALL LAGUER(A,M,X,EPS,POLISH)
		IF (N.EQ.0) THEN
			N=1
			Y(1)=X
			WRITE(*,'(1X,I5,2F15.6)') N,X
		ELSE
			IFLAG=0
			DO 11 J=1,N
				IF (CABS(X-Y(J)).LE.
     *					EPS*CABS(X)) IFLAG=1
11			CONTINUE
			IF (IFLAG.EQ.0) THEN
				N=N+1
				Y(N)=X
				WRITE(*,'(1X,I5,2F15.6)') N,X
			ENDIF
		ENDIF
12	CONTINUE
	END
	PROGRAM D14R2
C	Driver for routine LFIT
	PARAMETER(NPT=100,SPREAD=0.1,NTERM=3)
	DIMENSION X(NPT),Y(NPT),SIG(NPT),A(NTERM),COVAR(NTERM,NTERM),
     *			LISTA(NTERM)
	EXTERNAL FUNCS
	IDUM=-911
	DO 12 I=1,NPT
		X(I)=0.1*I
		Y(I)=FLOAT(NTERM)
		DO 11 J=NTERM-1,1,-1
			Y(I)=J+Y(I)*X(I)
11		CONTINUE
		Y(I)=Y(I)+SPREAD*GASDEV(IDUM)
		SIG(I)=SPREAD
12	CONTINUE
	MFIT=NTERM
	DO 13 I=1,MFIT
		LISTA(I)=I
13	CONTINUE
	CALL LFIT(X,Y,SIG,NPT,A,NTERM,LISTA,MFIT,COVAR,NTERM,CHISQ,FUNCS)
	WRITE(*,'(/1X,T4,A,T22,A)') 'Parameter','Uncertainty'
	DO 14 I=1,NTERM
		WRITE(*,'(1X,T5,A,I1,A,F8.6,F11.6)') 'A(',I,') = ',
     *			A(I),SQRT(COVAR(I,I))
14	CONTINUE
	WRITE(*,'(/3X,A,E12.6)') 'Chi-squared = ',CHISQ
	WRITE(*,'(/3X,A)') 'Full covariance matrix'
	DO 15 I=1,NTERM
		WRITE(*,'(1X,4E12.2)') (COVAR(I,J),J=1,NTERM)
15	CONTINUE
	WRITE(*,'(/1X,A)') 'press RETURN to continue...'
	READ(*,*)
C	Now test the LISTA feature
	DO 16 I=1,NTERM
		LISTA(I)=NTERM+1-I
16	CONTINUE
	CALL LFIT(X,Y,SIG,NPT,A,NTERM,LISTA,MFIT,COVAR,NTERM,CHISQ,FUNCS)
	WRITE(*,'(/1X,T4,A,T22,A)') 'Parameter','Uncertainty'
	DO 17 I=1,NTERM
		WRITE(*,'(1X,T5,A,I1,A,F8.6,F11.6)') 'A(',I,') = ',
     *			A(I),SQRT(COVAR(I,I))
17	CONTINUE
	WRITE(*,'(/3X,A,E12.6)') 'Chi-squared = ',CHISQ
	WRITE(*,'(/3X,A)') 'Full covariance matrix'
	DO 18 I=1,NTERM
		WRITE(*,'(1X,4E12.2)') (COVAR(I,J),J=1,NTERM)
18	CONTINUE
	WRITE(*,'(/1X,A)') 'press RETURN to continue...'
	READ(*,*)
C	Now check results of restricting fit parameters
	II=1
	DO 19 I=1,NTERM
		IF (MOD(I,2).EQ.1) THEN
			LISTA(II)=I
			II=II+1
		ENDIF
19	CONTINUE
	MFIT=II-1
	CALL LFIT(X,Y,SIG,NPT,A,NTERM,LISTA,MFIT,COVAR,NTERM,CHISQ,FUNCS)
	WRITE(*,'(/1X,T4,A,T22,A)') 'Parameter','Uncertainty'
	DO 21 I=1,NTERM
		WRITE(*,'(1X,T5,A,I1,A,F8.6,F11.6)') 'A(',I,') = ',A(I),
     *			SQRT(COVAR(I,I))
21	CONTINUE
	WRITE(*,'(/3X,A,E12.6)') 'Chi-squared = ',CHISQ
	WRITE(*,'(/3X,A)') 'Full covariance matrix'
	DO 22 I=1,NTERM
		WRITE(*,'(1X,4E12.2)') (COVAR(I,J),J=1,NTERM)
22	CONTINUE
	END

	SUBROUTINE FUNCS(X,AFUNC,MA)
	DIMENSION AFUNC(MA)
	AFUNC(1)=1.0
	DO 11 I=2,MA
		AFUNC(I)=X*AFUNC(I-1)
11	CONTINUE
	END
	PROGRAM D10R7
C	Driver for routine LINMIN
	PARAMETER(NDIM=3,PIO2=1.5707963)
	DIMENSION P(NDIM),XI(NDIM)
	WRITE(*,'(/1X,A)') 'Minimum of a 3-D quadratic centered'
	WRITE(*,'(1X,A)') 'at (1.0,1.0,1.0). Minimum is found'
	WRITE(*,'(1X,A)') 'along a series of radials.'
	WRITE(*,'(/1X,T10,A,T22,A,T34,A,T42,A/)') 'x','y','z','minimum'
	DO 11 I=0,10
		X=PIO2*I/10.0
		SR2=SQRT(2.0)
		XI(1)=SR2*COS(X)
		XI(2)=SR2*SIN(X)
		XI(3)=1.0
		P(1)=0.0
		P(2)=0.0
		P(3)=0.0
		CALL LINMIN(P,XI,NDIM,FRET)
		WRITE(*,'(1X,4F12.6)') (P(J),J=1,3),FRET
11	CONTINUE
	END

	FUNCTION FUNC(X)
	DIMENSION X(3)
	FUNC=0.0
	DO 11 I=1,3
		FUNC=FUNC+(X(I)-1.0)**2
11	CONTINUE
	END
	PROGRAM D3R5
C	Driver for routine LOCATE
	PARAMETER(N=100)
	DIMENSION XX(N)
C	Create array to be searched
	DO 11 I=1,N
		XX(I)=EXP(I/20.0)-74.0
11	CONTINUE
	WRITE(*,*) 'Result of:   j=0 indicates x too small'
	WRITE(*,*) '           j=100 indicates x too large'
	WRITE(*,'(T5,A7,T17,A1,T24,A5,T34,A7)') 'locate ','j'
     *  ,'xx(j)','xx(j+1)'
C	Do test
	DO 12 I=1,19
		X=-100.0+200.0*I/20.0
		CALL LOCATE(XX,N,X,J)
		IF (J.EQ.0) THEN
			WRITE(*,'(1X,F10.4,I6,A12,F12.6)') X,J,'lower lim',XX(J+1)
		ELSE IF (J.EQ.N) THEN
			WRITE(*,'(1X,F10.4,I6,F12.6,A12)') X,J,XX(J),'upper lim'
		ELSE
			WRITE(*,'(1X,F10.4,I6,2F12.6)') X,J,XX(J),XX(J+1)
		ENDIF
12	CONTINUE
	END
	PROGRAM D2R3
C	Driver for routine LUBKSB
	PARAMETER (NP=20)
	DIMENSION A(NP,NP),B(NP,NP),INDX(NP)
	DIMENSION C(NP,NP),X(NP)
	CHARACTER TXT*3
	OPEN(5,FILE='MATRX1.DAT',STATUS='OLD')
	READ(5,*)
10	READ(5,*)
	READ(5,*) N,M
	READ(5,*)
	READ(5,*) ((A(K,L), L=1,N), K=1,N)
	READ(5,*)
	READ(5,*) ((B(K,L), K=1,N), L=1,M)
C	Save matrix A for later testing
	DO 12 L=1,N
		DO 11 K=1,N
			C(K,L)=A(K,L)
11		CONTINUE
12	CONTINUE
C	Do LU decomposition
	CALL LUDCMP(C,N,NP,INDX,P)
C	Solve equations for each right-hand vector
	DO 16 K=1,M
		DO 13 L=1,N
			X(L)=B(L,K)
13		CONTINUE
		CALL LUBKSB(C,N,NP,INDX,X)
C	Test results with original matrix
		WRITE(*,*) 'Right-hand side vector:'
		WRITE(*,'(1X,6F12.6)') (B(L,K), L=1,N)
		WRITE(*,*) 'Result of matrix applied to sol''n vector'
		DO 15 L=1,N
			B(L,K)=0.0
			DO 14 J=1,N
				B(L,K)=B(L,K)+A(L,J)*X(J)
14			CONTINUE
15		CONTINUE
		WRITE(*,'(1X,6F12.6)') (B(L,K), L=1,N)
		WRITE(*,*) '***********************************'
16	CONTINUE
	WRITE(*,*) 'Press RETURN for next problem:'
	READ(*,*)
	READ(5,'(A3)') TXT
	IF (TXT.NE.'END') GOTO 10
	CLOSE(5)
	END
	PROGRAM D2R2
C	Driver for routine LUDCMP
	PARAMETER(NP=20)
	DIMENSION A(NP,NP),XL(NP,NP),XU(NP,NP)
	DIMENSION INDX(NP),JNDX(NP),X(NP,NP)
	CHARACTER TXT*3
	OPEN (5,FILE='MATRX1.DAT',STATUS='OLD')
	READ(5,*)
10	READ(5,*)
	READ(5,*) N,M
	READ(5,*)
	READ(5,*) ((A(K,L), L=1,N), K=1,N)
	READ(5,*)
	READ(5,*) ((X(K,L), K=1,N), L=1,M)
C	Print out A-matrix for comparison with product of lower
C	and upper decomposition matrices.
	WRITE(*,*) 'Original matrix:'
	DO 11 K=1,N
		WRITE(*,'(1X,6F12.6)') (A(K,L), L=1,N)
11	CONTINUE
C	Perform the decomposition
	CALL LUDCMP(A,N,NP,INDX,D)
C	Compose separately the lower and upper matrices
	DO 13 K=1,N
		DO 12 L=1,N
			IF (L.GT.K) THEN
				XU(K,L)=A(K,L)
				XL(K,L)=0.0
			ELSE IF (L.LT.K) THEN
				XU(K,L)=0.0
				XL(K,L)=A(K,L)
			ELSE
				XU(K,L)=A(K,L)
				XL(K,L)=1.0
			ENDIF
12		CONTINUE
13	CONTINUE
C	Compute product of lower and upper matrices for
C	comparison with original matrix.
	DO 16 K=1,N
		JNDX(K)=K
		DO 15 L=1,N
			X(K,L)=0.0
			DO 14 J=1,N
				X(K,L)=X(K,L)+XL(K,J)*XU(J,L)
14			CONTINUE
15		CONTINUE
16	CONTINUE
	WRITE(*,*) 'Product of lower and upper matrices (unscrambled):'
	DO 17 K=1,N
		DUM=JNDX(INDX(K))
		JNDX(INDX(K))=JNDX(K)
		JNDX(K)=DUM
17	CONTINUE
	DO 19 K=1,N
		DO 18 J=1,N
			IF (JNDX(J).EQ.K) THEN
				WRITE(*,'(1X,6F12.6)') (X(J,L), L=1,N)
			ENDIF
18		CONTINUE
19	CONTINUE
	WRITE(*,*) 'Lower matrix of the decomposition:'
	DO 21 K=1,N
		WRITE(*,'(1X,6F12.6)') (XL(K,L), L=1,N)
21	CONTINUE
	WRITE(*,*) 'Upper matrix of the decomposition:'
	DO 22 K=1,N
		WRITE(*,'(1X,6F12.6)') (XU(K,L), L=1,N)
22	CONTINUE
	WRITE(*,*) '***********************************'
	WRITE(*,*) 'Press RETURN for next problem:'
	READ(*,*)
	READ(5,'(A3)') TXT
	IF (TXT.NE.'END') GOTO 10
	CLOSE(5)
	END
