	PROGRAM D9R1
C	Driver for routine SCRSHO
	EXTERNAL BESSJ0
	WRITE(*,*) 'Graph of the Bessel Function J0:'
	CALL SCRSHO(BESSJ0)
	END
	PROGRAM D8R9
C	Driver for routine SHELL
	DIMENSION A(100)
	OPEN(5,FILE='TARRAY.DAT',STATUS='OLD')
	READ(5,*) (A(I),I=1,100)
	CLOSE(5)
C	Print original array
	WRITE(*,*) 'Original array:'
	DO 11 I=1,10
		WRITE(*,'(1X,10F6.2)') (A(10*(I-1)+J),J=1,10)
11	CONTINUE
C	Sort array
	CALL SHELL(100,A)
C	Print sorted array
	WRITE(*,*) 'Sorted array:'
	DO 12 I=1,10
		WRITE(*,'(1X,10F6.2)') (A(10*(I-1)+J),J=1,10)
12	CONTINUE
	END
	PROGRAM D16R1
C	Driver for routine SHOOT
C	Solves for eigenvalues of Spheroidal Harmonics. Both
C	Prolate and Oblate case are handled simultaneously, leading
C	to six first-order equations. Unknown to SHOOT, these are 
C	actually two independent sets of three coupled equations,
C	one set with c^2 positive and the other with c^2 negative.
	PARAMETER(NVAR=6,N2=2,DELTA=1.0E-3,EPS=1.0E-6,DX=1.0E-4)
	DIMENSION V(2),DELV(2),F(2),DV(2)
	COMMON C2,M,N,FACTR
1	WRITE(*,*) 'Input M,N,C-Squared (999 to end)'
	READ(*,*) M,N,C2
	IF (C2.EQ.999.) STOP
	IF ((N.LT.M).OR.(M.LT.0).OR.(N.LT.0)) GOTO 1
	FACTR=1.0
	IF (M.NE.0) THEN
		Q1=N
		DO 11 I=1,M
			FACTR=-0.5*FACTR*(N+I)*(Q1/I)
			Q1=Q1-1.0
11		CONTINUE
	ENDIF
	V(1)=N*(N+1)-M*(M+1)+C2/2.0
	V(2)=N*(N+1)-M*(M+1)-C2/2.0
	DELV(1)=DELTA*V(1)
	DELV(2)=DELV(1)
	H1=0.1
	HMIN=0.0
	X1=-1.0+DX
	X2=0.0
	WRITE(*,'(1X,T12,A,T36,A)') 'Prolate','Oblate'
	WRITE(*,'(1X,T6,A,T17,A,T30,A,T41,A)') 
     *		'Mu(M,N)','Error Est.','Mu(M,N)','Error Est.'
2	CALL SHOOT(NVAR,V,DELV,N2,X1,X2,EPS,H1,HMIN,F,DV)
	WRITE(*,'(1X,4F12.6)') V(1),DV(1),V(2),DV(2)
	IF ((ABS(DV(1)).GT.ABS(EPS*V(1))).OR.
     *		((DV(2)).GT.ABS(EPS*V(2)))) GOTO 2
	END

	SUBROUTINE LOAD(X1,V,Y)
	COMMON C2,M,N,FACTR
	DIMENSION V(2),Y(6)
	Y(3)=V(1)
	Y(2)=-(Y(3)-C2)*FACTR/2.0/(M+1.0)
	Y(1)=FACTR+Y(2)*DX
	Y(6)=V(2)
	Y(5)=-(Y(6)+C2)*FACTR/2.0/(M+1.0)
	Y(4)=FACTR+Y(5)*DX
	RETURN
	END

	SUBROUTINE SCORE(X2,Y,F)
	COMMON C2,M,N
	DIMENSION Y(6),F(2)
	IF (MOD((N-M),2).EQ.0) THEN
		F(1)=Y(2)
		F(2)=Y(5)
	ELSE 
		F(1)=Y(1)
		F(2)=Y(4)
	ENDIF
	RETURN
	END

	SUBROUTINE DERIVS(X,Y,DYDX)
	COMMON C2,M,N
	DIMENSION Y(6),DYDX(6)
	DYDX(1)=Y(2)
	DYDX(3)=0.0
	DYDX(2)=(2.0*X*(M+1.0)*Y(2)-(Y(3)-C2*X*X)*Y(1))/(1.0-X*X)
	DYDX(4)=Y(5)
	DYDX(6)=0.0
	DYDX(5)=(2.0*X*(M+1.0)*Y(5)-(Y(6)+C2*X*X)*Y(4))/(1.0-X*X)
	RETURN
	END
	PROGRAM D16R2
C	Driver for routine SHOOTF
	PARAMETER(NVAR=3,N1=2,N2=1,DELTA=1.0E-3,EPS=1.0E-6,DXX=1.0E-4)
	DIMENSION V1(N2),DELV1(N2),V2(N1),DELV2(N1),
     *		DV1(N2),DV2(N1),F(NVAR)
	COMMON C2,M,N,FACTR,DX
1	WRITE(*,*) 'INPUT M,N,C-SQUARED (999 TO STOP)'
	READ(*,*) M,N,C2
	IF (C2.EQ.999.) STOP
	IF ((N.LT.M).OR.(M.LT.0)) THEN
		WRITE(*,*) 'Improper arguments'
		GOTO 1
	ENDIF
	FACTR=1.0
	IF (M.NE.0) THEN
		Q1=N
		DO 11 I=1,M
			FACTR=-0.5*FACTR*(N+I)*(Q1/I)
			Q1=Q1-1.0
11		CONTINUE
	ENDIF
	DX=DXX
	V1(1)=N*(N+1)-M*(M+1)+C2/2.0
	IF (MOD(N-M,2).EQ.0) THEN
		V2(1)=SIGN(1.0,FACTR)
	ELSE
		V2(1)=-SIGN(1.0,FACTR)
	ENDIF	
	V2(2)=V1(1)+1.0
	DELV1(1)=DELTA*V1(1)
	DELV2(1)=DELTA*FACTR
	DELV2(2)=DELV1(1)
	H1=0.1
	HMIN=0.0
	X1=-1.0+DX
	X2=1.0-DX
	XF=0.0
	WRITE(*,'(/1X,T20,A,T40,A,T60,A)') 'Mu(-1)','Y(1-dx)','Mu(+1)'
2	CALL SHOOTF(NVAR,V1,V2,DELV1,DELV2,N1,N2,X1,X2,
     *		XF,EPS,H1,HMIN,F,DV1,DV2)
	WRITE(*,'(/4X,A,3F20.6)') 'V ',V1(1),V2(1),V2(2)
	WRITE(*,'(4X,A,3F20.6)') 'DV',DV1(1),DV2(1),DV2(2)
	IF (ABS(DV1(1)).GT.ABS(EPS*V1(1))) GOTO 2
	END

	SUBROUTINE LOAD1(X1,V1,Y)
	COMMON C2,M,N,FACTR,DX
	DIMENSION V1(1),Y(3)
	Y(3)=V1(1)
	Y(2)=-(Y(3)-C2)*FACTR/2.0/(M+1.0)
	Y(1)=FACTR+Y(2)*DX
	RETURN
	END

	SUBROUTINE LOAD2(X2,V2,Y)
	COMMON C2,M,N
	DIMENSION V2(2),Y(3)
	Y(3)=V2(2)
	Y(1)=V2(1)
	Y(2)=(Y(3)-C2)*Y(1)/2.0/(M+1.0)
	RETURN
	END

	SUBROUTINE SCORE(XF,Y,F)
	COMMON C2,M,N
	DIMENSION Y(3),F(3)
	DO 11 I=1,3
		F(I)=Y(I)
11	CONTINUE
	RETURN
	END

	SUBROUTINE DERIVS(X,Y,DYDX)
	COMMON C2,M,N
	DIMENSION Y(3),DYDX(3)
	DYDX(1)=Y(2)
	DYDX(3)=0.0
	DYDX(2)=(2.0*X*(M+1.0)*Y(2)-(Y(3)-C2*X*X)*Y(1))/(1.0-X*X)
	RETURN
	END
	PROGRAM D10R12
C	Driver for routine SIMPLX
C	Incorporates examples discussed in text
	PARAMETER(N=4,M=4,NP=5,MP=6,M1=2,M2=1,M3=1,NM1M2=N+M1+M2)
	DIMENSION A(MP,NP),IZROV(N),IPOSV(M),ANUM(NP)
	CHARACTER TXT(NM1M2)*2,ALPHA(NP)*2
	DATA TXT/'x1','x2','x3','x4','y1','y2','y3'/
	DATA A/0.0,740.0,0.0,0.5,9.0,0.0,1.0,-1.0,0.0,0.0,-1.0,0.0,
     *		1.0,0.0,-2.0,-1.0,-1.0,0.0,3.0,-2.0,0.0,1.0,-1.0,0.0,
     *		-0.5,0.0,7.0,-2.0,-1.0,0.0/
	CALL SIMPLX(A,M,N,MP,NP,M1,M2,M3,ICASE,IZROV,IPOSV)
	IF (ICASE.EQ.1) THEN
		WRITE(*,*) 'Unbounded objective function'
	ELSE IF (ICASE.EQ.-1) THEN
		WRITE(*,*) 'No solutions satisfy constraints given'
	ELSE
		JJ=1
		DO 11 I=1,N
			IF (IZROV(I).LE.(N+M1+M2)) THEN
				ALPHA(JJ)=TXT(IZROV(I))
				JJ=JJ+1
			ENDIF
11		CONTINUE
		JMAX=JJ-1
		WRITE(*,'(/3X,5A10)') '  ',(ALPHA(JJ),JJ=1,JMAX)		
		DO 13 I=1,M+1
			IF (I.GT.1) THEN 
				ALPHA(1)=TXT(IPOSV(I-1))
			ELSE
				ALPHA(1)='  '
			ENDIF
			ANUM(1)=A(I,1)
			JJ=2
			DO 12 J=2,N+1
				IF (IZROV(J-1).LE.(N+M1+M2)) THEN
					ANUM(JJ)=A(I,J)
					JJ=JJ+1
				ENDIF
12			CONTINUE
			JMAX=JJ-1
			WRITE(*,'(1X,A3,(5F10.2))') ALPHA(1),
     *				(ANUM(JJ),JJ=1,JMAX)
13		CONTINUE
	ENDIF
	END
	PROGRAM D12R4
C	Driver for routine SINFT
	PARAMETER(EPS=1.0E-3,NP=16,WIDTH=30.0,PI=3.14159)
	DIMENSION DATA(NP),SIZE(NP)
1	WRITE(*,'(1X,A,I2,A)') 'Period of sinusoid in channels (2-',NP,')'
	READ(*,*) PER
	IF (PER.LE.0.) STOP
	DO 11 I=1,NP
		DATA(I)=SIN(2.0*PI*(I-1)/PER)
11	CONTINUE
	CALL SINFT(DATA,NP)
	BIG=-1.0E10
	SMALL=1.0E10
	DO 12 I=1,NP
		IF (DATA(I).LT.SMALL) SMALL=DATA(I)
		IF (DATA(I).GT.BIG) BIG=DATA(I)
12	CONTINUE
	SCAL=WIDTH/(BIG-SMALL)
	DO 13 I=1,NP
		NLIM=SCAL*(DATA(I)-SMALL)+EPS
		WRITE(*,'(1X,I4,1X,60A1)') I,('*',J=1,NLIM+1)
13	CONTINUE
	WRITE(*,*) 'press continue ...'
	READ(*,*)
	CALL SINFT(DATA,NP)
	BIG=-1.0E10
	SMALL=1.0E10
	DO 14 I=1,NP
		IF(DATA(I).LT.SMALL) SMALL=DATA(I)
		IF(DATA(I).GT.BIG) BIG=DATA(I)
14	CONTINUE
	SCAL=WIDTH/(BIG-SMALL)
	DO 15 I=1,NP
		NLIM=SCAL*(DATA(I)-SMALL)+EPS
		WRITE(*,'(1X,I4,1X,60A1)') I,('*',J=1,NLIM+1)
15	CONTINUE
	GOTO 1
	END
	PROGRAM D13R20
C	Driver for routine SMOOFT
	PARAMETER(MMAX=1024,N=100,HASH=0.05,SCALE=100.0,PTS=10.0)
	DIMENSION Y(MMAX)
	CHARACTER TEXT(64)*1
	IDUM=-7
	DO 11 I=1,N
		Y(I)=3.0*I/N*EXP(-3.0*I/N)
		Y(I)=Y(I)+HASH*GASDEV(IDUM)
11	CONTINUE
	DO 15 K=1,3
		NSTP=N/20
		WRITE(*,'(/3X,A,T16,A)') 'Data:','Graph:'
		DO 14 I=1,N,NSTP
			DO 12 J=1,64
				TEXT(J)=' '
12			CONTINUE
			IBAR=NINT(SCALE*Y(I))
			DO 13 J=1,64
				IF (J.LE.IBAR) TEXT(J)='*'
13			CONTINUE
			WRITE(*,'(1X,F10.6,4X,64A1)') Y(I),
     *				(TEXT(J),J=1,64)
14		CONTINUE
		WRITE(*,'(/1X,A)') ' press RETURN to smooth ...'
		READ(*,*)
		CALL SMOOFT(Y,N,PTS)
15	CONTINUE
	END
	PROGRAM D6R30
C	Driver for routine SNCNDN
	CHARACTER TEXT*26
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)',ERR=99) TEXT
	IF (TEXT.EQ.'Jacobian Elliptic Function') THEN
		READ(5,*) NVAL
		WRITE(*,*) TEXT
		WRITE(*,'(1X,T4,A,T13,A,T21,A,T38,A,T49,A,T60,A)')
     *			'Mc','U','Actual','SN','SN^2+CN^2',
     *			'(Mc)*(SN^2)+DN^2'
		DO 11 I=1,NVAL
			READ(5,*) EM,UU,VALUE
			EMMC=1.0-EM
			CALL SNCNDN(UU,EMMC,SN,CN,DN)
			RESULT1=SN*SN+CN*CN
			RESULT2=EM*SN*SN+DN*DN
			WRITE(*,'(1X,F5.2,F8.2,2E15.5,F12.5,F14.5)')
     *				EMMC,UU,VALUE,SN,RESULT1,RESULT2
11		CONTINUE
	ENDIF
	GOTO 10
99	CLOSE(5)
	END
	PROGRAM D17R1
C	Driver for routine SOR
	PARAMETER(JMAX=11,PI=3.1415926)
	IMPLICIT REAL*8(A-H,O-Z)
	DIMENSION A(JMAX,JMAX),B(JMAX,JMAX),C(JMAX,JMAX),
     *		D(JMAX,JMAX),E(JMAX,JMAX),F(JMAX,JMAX),U(JMAX,JMAX)
	DO 12 I=1,JMAX
	      	DO 11 J=1,JMAX
			A(I,J)=1.0
			B(I,J)=1.0
			C(I,J)=1.0
			D(I,J)=1.0
			E(I,J)=-4.0
			F(I,J)=0.0
			U(I,J)=0.0
11		CONTINUE
12	CONTINUE
	MIDL=JMAX/2+1
	F(MIDL,MIDL)=2.0
	RJAC=COS(PI/JMAX)
	CALL SOR(A,B,C,D,E,F,U,JMAX,RJAC)
	WRITE(*,'(1X,A)') 'SOR Solution:'
	DO 13 I=1,JMAX
		WRITE(*,'(1X,11F6.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
			F(I,J)=U(I+1,J)+U(I-1,J)+U(I,J+1)
     *				+U(I,J-1)-4.0*U(I,J)
14		CONTINUE
		WRITE(*,'(7X,11F6.2)') (F(I,J),J=2,JMAX-1)
15	CONTINUE
	END
	PROGRAM D8R3
C	Driver for routine SORT
	DIMENSION A(100)
	OPEN(5,FILE='TARRAY.DAT',STATUS='OLD')
	READ(5,*) (A(I),I=1,100)
	CLOSE(5)
C	Print original array
	WRITE(*,*) 'Original array:'
	DO 11 I=1,10
		WRITE(*,'(1X,10F6.2)') (A(10*(I-1)+J),J=1,10)
11	CONTINUE
C	Sort array
	CALL SORT(100,A)
C	Print sorted array
	WRITE(*,*) 'Sorted array:'
	DO 12 I=1,10
		WRITE(*,'(1X,10F6.2)') (A(10*(I-1)+J),J=1,10)
12	CONTINUE
	END
	PROGRAM D8R4
C	Driver for routine SORT2
	DIMENSION A(100),B(100)
	OPEN(5,FILE='TARRAY.DAT',STATUS='OLD')
	READ(5,*) (A(I),I=1,100)
	CLOSE(5)
C	Generate B-array
	DO 11 I=1,100
		B(I)=I
11	CONTINUE
C	Sort A and mix B
	CALL SORT2(100,A,B)
	WRITE(*,*) 'After sorting A and mixing B, array A is:'
	DO 12 I=1,10
		WRITE(*,'(1X,10F6.2)') (A(10*(I-1)+J), J=1,10)
12	CONTINUE
	WRITE(*,*) '...and array B is:'
	DO 13 I=1,10
		WRITE(*,'(1X,10F6.2)') (B(10*(I-1)+J), J=1,10)
13	CONTINUE
	WRITE(*,*) 'press RETURN to continue...'
	READ(*,*)
C	Sort B and mix A
	CALL SORT2(100,B,A)
	WRITE(*,*) 'After sorting B and mixing A, array A is:'
	DO 14 I=1,10
		WRITE(*,'(1X,10F6.2)') (A(10*(I-1)+J), J=1,10)
14	CONTINUE
	WRITE(*,*) '...and array B is:'
	DO 15 I=1,10
		WRITE(*,'(1X,10F6.2)') (B(10*(I-1)+J), J=1,10)
15	CONTINUE
	END
	PROGRAM D8R6
C	Driver for routine SORT3
	PARAMETER(NLEN=64)
	DIMENSION A(NLEN),IB(NLEN),IC(NLEN),WKSP(NLEN),INDX(NLEN)
	CHARACTER MSG1*33,MSG2*31
	CHARACTER MSG*64,AMSG(64)*1,BMSG(64)*1,CMSG(64)*1
	EQUIVALENCE(MSG,AMSG(1)),(MSG1,AMSG(1)),(MSG2,AMSG(34))
	DATA MSG1/'I''d rather have a bottle in front'/
	DATA MSG2/' of me than a frontal lobotomy.'/
	WRITE(*,*) 'Original message:'
	WRITE(*,'(1X,64A1,/)') (AMSG(J),J=1,64)
C	Read array of random numbers
	OPEN(5,FILE='TARRAY.DAT',STATUS='OLD')
	READ(5,*) (A(I),I=1,NLEN)
	CLOSE(5)
C	Create array IB and array IC
	DO 11 I=1,NLEN
		IB(I)=I
		IC(I)=NLEN+1-I
11	CONTINUE
C	Sort array A while mixing IB and IC
	CALL SORT3(NLEN,A,IB,IC,WKSP,INDX)
C	Scramble message according to array IB
	DO 12 I=1,NLEN
		J=IB(I)
		BMSG(I)=AMSG(J)
12	CONTINUE
	WRITE(*,*) 'Scrambled message:'
	WRITE(*,'(1X,64A1,/)') (BMSG(J),J=1,64)
C	Unscramble according to array C
	DO 13 I=1,NLEN
		J=IC(I)
		CMSG(J)=BMSG(I)
13	CONTINUE
	WRITE(*,*) 'Mirrored message:'
	WRITE(*,'(1X,64A1,/)') (CMSG(J),J=1,64)
	END
	PROGRAM D2R10
C	Driver for SPARSE
	PARAMETER(N=20)
	COMMON M
	DIMENSION B(N),X(N),BCMP(N)
	EXTERNAL ASUB,ATSUB
	M=N
	DO 11 I=1,N
		X(I)=0.0
		B(I)=1.0
11	CONTINUE
	B(1)=3.0
	B(N)=-1.0
	CALL SPARSE(B,N,ASUB,ATSUB,X,RSQ)
	WRITE(*,'(/1X,A,E15.6)') 'Sum-squared residual:',RSQ
	WRITE(*,'(/1X,A)') 'Solution vector:'
	WRITE(*,'(1X,5F12.6)') (X(I),I=1,N)
	CALL ASUB(X,BCMP)
	WRITE(*,'(/1X,A)') 'press RETURN to continue...'
	READ(*,*)
	WRITE(*,'(1X,A/T8,A,T22,A)') 'Test of solution vector:','a*x','b'
	DO 12 I=1,N
		WRITE(*,'(1X,2F12.6)') BCMP(I),B(I)
12	CONTINUE
	END


	SUBROUTINE ASUB(XIN,XOUT)
	COMMON N
	DIMENSION XIN(N),XOUT(N)
	XOUT(1)=XIN(1)+2.0*XIN(2)
	XOUT(N)=-2.0*XIN(N-1)+XIN(N)
	DO 11 I=2,N-1
		XOUT(I)=-2.0*XIN(I-1)+XIN(I)+2.0*XIN(I+1)
11	CONTINUE
	RETURN
	END


	SUBROUTINE ATSUB(XIN,XOUT)
	COMMON N
	DIMENSION XIN(N),XOUT(N)
	XOUT(1)=XIN(1)-2.0*XIN(2)
	XOUT(N)=2.0*XIN(N-1)+XIN(N)
	DO 11 I=2,N-1
		XOUT(I)=2.0*XIN(I-1)+XIN(I)-2.0*XIN(I+1)
11	CONTINUE
	RETURN
	END
	PROGRAM D12R8
C	Driver for routine SPCTRM
	PARAMETER(M=16,M4=4*M)
	DIMENSION P(M),Q(M),W1(M4),W2(M)
	LOGICAL OVRLAP
	OPEN(9,FILE='SPCTRL.DAT',STATUS='OLD')
	K=8
	OVRLAP=.TRUE.
	CALL SPCTRM(P,M,K,OVRLAP,W1,W2)
	REWIND(9)
	K=16
	OVRLAP=.FALSE.
	CALL SPCTRM(Q,M,K,OVRLAP,W1,W2)
	CLOSE(9)
	WRITE(*,*) 'Spectrum of data in file SPCTRL.DAT'
	WRITE(*,'(1X,T14,A,T29,A)') 'Overlapped','Non-Overlapped'
	DO 11 J=1,M
		WRITE(*,'(1X,I4,2F17.6)') J,P(J),Q(J)
11	CONTINUE
	END
	PROGRAM D13R16
C	Driver for routine SPEAR
	PARAMETER(NDAT=20,NMON=12)
	DIMENSION DATA1(NDAT),DATA2(NDAT),RAYS(NDAT,NMON)
	DIMENSION WKSP1(NDAT),WKSP2(NDAT),AVE(NDAT),ZLAT(NDAT)
	CHARACTER CITY(NDAT)*15,MON(NMON)*4,TEXT*64
	OPEN(5,FILE='TABLE2.DAT',STATUS='OLD')
	READ(5,*)
	READ(5,'(A)') TEXT
	READ(5,'(15X,12A4/)') (MON(I),I=1,12)
	DO 11 I=1,NDAT
		READ(5,'(A15,12F4.0,F6.0,F6.1)')
     *			CITY(I),(RAYS(I,J),J=1,12),AVE(I),ZLAT(I)
11	CONTINUE
	CLOSE(5)
	WRITE(*,*) TEXT
	WRITE(*,'(1X,15X,12A4)') (MON(I),I=1,12)
	DO 12 I=1,NDAT
		WRITE(*,'(1X,A,12I4,I6,F6.1)') CITY(I),
     *			(NINT(RAYS(I,J)),J=1,12)
12	CONTINUE
C	Check temperature correlations between different months
	WRITE(*,'(/1X,A)') 
     *		'Are sunny summer places also sunny winter places?'
	WRITE(*,'(1X,2A)') 'Check correlation of sampled U.S. solar ',
     *  'radiation (july with other months)'
	WRITE(*,'(/1X,A,T16,A,T23,A,T37,A,T49,A,T63,A/)') 
     *		'Month','D','St. Dev.','PROBD',
     *		'Spearman R','PROBRS'
	DO 13 I=1,NDAT
		DATA1(I)=RAYS(I,1)
13	CONTINUE
	DO 15 J=1,12
		DO 14 I=1,NDAT
			DATA2(I)=RAYS(I,J)
14		CONTINUE
		CALL SPEAR(DATA1,DATA2,NDAT,WKSP1,WKSP2,D,ZD,PROBD,RS,PROBRS)
		WRITE(*,'(1X,A,F13.2,2F12.6,3X,2F12.6)') 
     *			MON(J),D,ZD,PROBD,RS,PROBRS
15	CONTINUE
	END
	PROGRAM D3R12
C	Driver for routine SPLIE2
	PARAMETER(M=10,N=10)
	DIMENSION X1(M),X2(N),Y(M,N),Y2(M,N)
	DO 11 I=1,M
		X1(I)=0.2*I
11	CONTINUE
	DO 12 I=1,N
		X2(I)=0.2*I
12	CONTINUE
	DO 14 I=1,M
		DO 13 J=1,N
			X1X2=X1(I)*X2(J)
			Y(I,J)=X1X2**2
13		CONTINUE
14	CONTINUE
	CALL SPLIE2(X1,X2,Y,M,N,Y2)
	WRITE(*,'(/1X,A)') 'Second derivatives from SPLIE2'
	WRITE(*,'(1X,A/)') 'Natural spline assumed'
	DO 15 I=1,5
		WRITE(*,'(1X,5F12.6)') (Y2(I,J),J=1,5)
15	CONTINUE
	WRITE(*,'(/1X,A/)') 'Actual second derivatives'
	DO 17 I=1,5
		DO 16 J=1,5
			Y2(I,J)=2.0*(X1(I)**2)
16		CONTINUE
		WRITE(*,'(1X,5F12.6)') (Y2(I,J),J=1,5)
17	CONTINUE
	END
	PROGRAM D3R13
C	Driver for routine SPLIN2
	PARAMETER(M=10,N=10)
	DIMENSION X1(M),X2(N),Y(M,N),Y2(M,N)
	DO 11 I=1,M
		X1(I)=0.2*I
11	CONTINUE
	DO 12 I=1,N
		X2(I)=0.2*I
12	CONTINUE
	DO 14 I=1,M
		DO 13 J=1,N
			X1X2=X1(I)*X2(J)
			Y(I,J)=X1X2*EXP(-X1X2)
13		CONTINUE
14	CONTINUE
	CALL SPLIE2(X1,X2,Y,M,N,Y2)
	WRITE(*,'(/1X,T9,A,T21,A,T31,A,T43,A)') 
     *		'x1','x2','splin2','actual'
	DO 15 I=1,10
		XX1=0.1*I
		XX2=XX1**2
		CALL SPLIN2(X1,X2,Y,Y2,M,N,XX1,XX2,F)
		X1X2=XX1*XX2
		FF=X1X2*EXP(-X1X2)
		WRITE(*,'(1X,4F12.6)') XX1,XX2,F,FF
15	CONTINUE
	END
	PROGRAM D3R3
C	Driver for routine SPLINE
	PARAMETER(N=20,PI=3.141593)
	DIMENSION X(N),Y(N),Y2(N)
	WRITE(*,*) 'Second-derivatives for sin(x) from 0 to PI'
C	Generate array for interpolation
	DO 11 I=1,20
		X(I)=I*PI/N
		Y(I)=SIN(X(I))
11	CONTINUE
C	Calculate 2nd derivative with SPLINE
	YP1=COS(X(1))
	YPN=COS(X(N))
	CALL SPLINE(X,Y,N,YP1,YPN,Y2)
C	Test result
	WRITE(*,'(T19,A,T35,A)') 'spline','actual'
	WRITE(*,'(T6,A,T17,A,T33,A)') 'angle','2nd deriv','2nd deriv'
	DO 12 I=1,N
		WRITE(*,'(1X,F8.2,2F16.6)') X(I),Y2(I),-SIN(X(I))
12	CONTINUE
	END
	PROGRAM D3R4
C	Driver for routine SPLINT, which calls SPLINE
	PARAMETER(NP=10,PI=3.141593)
	DIMENSION XA(NP),YA(NP),Y2(NP)
	DO 14 NFUNC=1,2
		IF (NFUNC.EQ.1) THEN
			WRITE(*,*) 'Sine function from 0 to pi'
			DO 11 I=1,NP
				XA(I)=I*PI/NP
				YA(I)=SIN(XA(I))
11			CONTINUE
			YP1=COS(XA(1))
			YPN=COS(XA(NP))
		ELSE IF (NFUNC.EQ.2) THEN
			WRITE(*,*) 'Exponential function from 0 to 1'
			DO 12 I=1,NP
				XA(I)=1.0*I/NP
				YA(I)=EXP(XA(I))
12			CONTINUE
			YP1=EXP(XA(1))
			YPN=EXP(XA(NP))
		ELSE
			STOP
		ENDIF
C	Call SPLINE to get second derivatives
		CALL SPLINE(XA,YA,NP,YP1,YP2,Y2)
C	Call SPLINT for interpolations
		WRITE(*,'(1X,T10,A1,T20,A4,T28,A13)') 'x','f(x)','interpolation'
		DO 13 I=1,10
			IF (NFUNC.EQ.1) THEN
				X=(-0.05+I/10.0)*PI
				F=SIN(X)
			ELSE IF (NFUNC.EQ.2) THEN
				X=-0.05+I/10.0
				F=EXP(X)
			ENDIF
			CALL SPLINT(XA,YA,Y2,NP,X,Y)
			WRITE(*,'(1X,3F12.6)') X,F,Y
13		CONTINUE
		WRITE(*,*) '***********************************'
		WRITE(*,*) 'Press RETURN'
		READ(*,*)
14	CONTINUE
	END
	PROGRAM D2R8
C	Driver for routine SVBKSB, which calls routine SVDCMP
	PARAMETER(NP=20)
	DIMENSION A(NP,NP),B(NP,NP),U(NP,NP),W(NP)
	DIMENSION V(NP,NP),C(NP),X(NP)
	CHARACTER DUMMY*3
	OPEN(5,FILE='MATRX1.DAT',STATUS='OLD')
10	READ(5,'(A)') DUMMY
	IF (DUMMY.EQ.'END') GOTO 99
	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	Copy A into U
	DO 12 K=1,N
		DO 11 L=1,N
			U(K,L)=A(K,L)
11		CONTINUE
12	CONTINUE
C	Decompose matrix A
	CALL SVDCMP(U,N,N,NP,NP,W,V)
C	Find maximum singular value
	WMAX=0.0
	DO 13 K=1,N
		IF (W(K).GT.WMAX) WMAX=W(K)
13	CONTINUE
C	Define "small"
	WMIN=WMAX*(1.0E-6)
C	Zero the "small" singular values
	DO 14 K=1,N
		IF (W(K).LT.WMIN) W(K)=0.0
14	CONTINUE
C	Backsubstitute for each right-hand side vector
	DO 18 L=1,M
		WRITE(*,'(1X,A,I2)') 'Vector number ',L
		DO 15 K=1,N
			C(K)=B(K,L)
15		CONTINUE
		CALL SVBKSB(U,W,V,N,N,NP,NP,C,X)
		WRITE(*,*) '    Solution vector is:'
		WRITE(*,'(1X,6F12.6)') (X(K), K=1,N)
		WRITE(*,*) '    Original right-hand side vector:'
		WRITE(*,'(1X,6F12.6)') (C(K), K=1,N)
		WRITE(*,*) '    Result of (matrix)*(sol''n vector):'
		DO 17 K=1,N
			C(K)=0.0
			DO 16 J=1,N
				C(K)=C(K)+A(K,J)*X(J)
16			CONTINUE
17		CONTINUE
		WRITE(*,'(1X,6F12.6)') (C(K), K=1,N)
18	CONTINUE
	WRITE(*,*) '***********************************'
	WRITE(*,*) 'Press RETURN for next problem'
	READ(*,*)
	GOTO 10
99	CLOSE(5)
	END
	PROGRAM D2R9
C	Driver for routine SVDCMP
	PARAMETER(NP=20)
	DIMENSION A(NP,NP),U(NP,NP),W(NP),V(NP,NP)
	CHARACTER DUMMY*3
	OPEN(5,FILE='MATRX3.DAT',STATUS='OLD')
10	READ(5,'(A)') DUMMY
	IF (DUMMY.EQ.'END') GOTO 99
	READ(5,*)
	READ(5,*) M,N
	READ(5,*)
C	Copy original matrix into U
	DO 12 K=1,M
		READ(5,*) (A(K,L), L=1,N)
		DO 11 L=1,N
			U(K,L)=A(K,L)
11		CONTINUE
12	CONTINUE
	IF (N.GT.M) THEN
		DO 14 K=M+1,N
			DO 13 L=1,N
				A(K,L)=0.0
				U(K,L)=0.0
13			CONTINUE
14		CONTINUE
		M=N
	ENDIF
C	Perform decomposition
	CALL SVDCMP(U,M,N,NP,NP,W,V)
C	Print results
	WRITE(*,*) 'Decomposition Matrices:'
	WRITE(*,*) 'Matrix U'
	DO 15 K=1,M
		WRITE(*,'(1X,6F12.6))') (U(K,L),L=1,N)
15	CONTINUE
	WRITE(*,*) 'Diagonal of Matrix W'
	WRITE(*,'(1X,6F12.6))') (W(K),K=1,N)
	WRITE(*,*) 'Matrix V-Transpose'
	DO 16 K=1,N
		WRITE(*,'(1X,6F12.6))') (V(L,K),L=1,N)
16	CONTINUE
	WRITE(*,*) 'Check product against original matrix:'
	WRITE(*,*) 'Original Matrix:'
	DO 17 K=1,M
		WRITE(*,'(1X,6F12.6))') (A(K,L),L=1,N)
17	CONTINUE
	WRITE(*,*) 'Product U*W*(V-Transpose):'
	DO 21 K=1,M
		DO 19 L=1,N
			A(K,L)=0.0
			DO 18 J=1,N
				A(K,L)=A(K,L)+U(K,J)*W(J)*V(L,J)
18			CONTINUE
19		CONTINUE
		WRITE(*,'(1X,6F12.6))') (A(K,L),L=1,N)
21	CONTINUE
	WRITE(*,*) '***********************************'
	WRITE(*,*) 'Press RETURN for next problem'
	READ(*,*)
	GOTO 10
99	CLOSE(5)
	END
	PROGRAM D14R4
C	Driver for routine SVDFIT
	EXTERNAL FPOLY,FLEG
	PARAMETER(NPT=100,SPREAD=0.02,NPOL=5)
	DIMENSION X(NPT),Y(NPT),SIG(NPT),A(NPOL),CVM(NPOL,NPOL)
	DIMENSION U(NPT,NPOL),V(NPOL,NPOL),W(NPOL)
C	Polynomial fit
	IDUM=-911
	MP=NPT
	NP=NPOL
	DO 11 I=1,NPT
		X(I)=0.02*I
		Y(I)=1.0+X(I)*(2.0+X(I)*(3.0+X(I)*(4.0+X(I)*5.0)))
		Y(I)=Y(I)*(1.0+SPREAD*GASDEV(IDUM))
		SIG(I)=Y(I)*SPREAD
11	CONTINUE
	CALL SVDFIT(X,Y,SIG,NPT,A,NPOL,U,V,W,MP,NP,CHISQ,FPOLY)
	CALL SVDVAR(V,NPOL,NP,W,CVM,NPOL)
	WRITE(*,*) 'Polynomial fit:'
	DO 12 I=1,NPOL
		WRITE(*,'(1X,F12.6,A,F10.6)') A(I),'  +-',SQRT(CVM(I,I))
12	CONTINUE
	WRITE(*,'(1X,A,F12.6/)') 'Chi-squared',CHISQ
	CALL SVDFIT(X,Y,SIG,NPT,A,NPOL,U,V,W,MP,NP,CHISQ,FLEG)
	CALL SVDVAR(V,NPOL,NP,W,CVM,NPOL)
	WRITE(*,*) 'Legendre polynomial fit'
	DO 13 I=1,NPOL
		WRITE(*,'(1X,F12.6,A,F10.6)') A(I),'  +-',SQRT(CVM(I,I))
13	CONTINUE
	WRITE(*,'(1X,A,F12.6/)') 'Chi-squared',CHISQ
	END
	PROGRAM D14R5
C	Driver for routine SVDVAR
	PARAMETER(MP=6,MA=3,NCVM=MA)
	DIMENSION V(MP,MP),W(MP),CVM(NCVM,NCVM),TRU(MA,MA)
	DATA W/0.0,1.0,2.0,3.0,4.0,5.0/
	DATA V/1.0,2.0,3.0,4.0,5.0,6.0,1.0,2.0,3.0,4.0,5.0,6.0,
     *		1.0,2.0,3.0,4.0,5.0,6.0,1.0,2.0,3.0,4.0,5.0,6.0,
     *		1.0,2.0,3.0,4.0,5.0,6.0,1.0,2.0,3.0,4.0,5.0,6.0/
	DATA TRU/1.25,2.5,3.75,2.5,5.0,7.5,3.75,7.5,11.25/
	WRITE(*,'(/1X,A)') 'Matrix V'
	DO 11 I=1,MP
		WRITE(*,'(1X,6F12.6)') (V(I,J),J=1,MP)
11	CONTINUE
	WRITE(*,'(/1X,A)') 'Vector W'
	WRITE(*,'(1X,6F12.6)') (W(I),I=1,MP)
	CALL SVDVAR(V,MA,MP,W,CVM,NCVM)
	WRITE(*,'(/1X,A)') 'Covariance matrix from SVDVAR'
	DO 12 I=1,MA
		WRITE(*,'(1X,3F12.6)') (CVM(I,J),J=1,MA)
12	CONTINUE
	WRITE(*,'(/1X,A)') 'Expected covariance matrix'
	DO 13 I=1,MA
		WRITE(*,'(1X,3F12.6)') (TRU(I,J),J=1,MA)
13	CONTINUE
	END
