      parameter(nmax=80,npar=2,n2max=npar*(npar+1)/2)
      implicit double precision (a-h,o-z)
      double precision s(nmax),y(nmax),p(npar),funcm1,var(n2max),
     +step(npar),stopcr,simp,a,b,k,s1,s2,zero,funcm2,funcm3,k0,b0,
     + pval12,pval32,lrt1v2,lrt3v2
     +,zchi,pleft,pright,pdf,seratab,ratab
      integer type(nmax),n
      common s,y,type,n
      external loglk1,loglk2,loglk3
      data zero/0.0d0/
c
c       Define the input and output files
c       Both files are text files
c       Input file comprises three columns
c       col 1:  substrate concentration
c       col 2:  rate of methane production
c       col 3:  type of methane: must be 12 for 12CH4 and 13 for 13CH4
c
c       while the numbers in columns 1 and 2 can contain decimal points
c       the numbers in column 3 must be 12 or 13 - 12.0 or 13.0 etc
c       will cause an error
c       columns should be separated by one or more spaces
c
c
c       Now set some constants needed by the minimisation routine
c
      nop=2
      n2=nop*(nop+1)/2
      step(1)=1.0d-4
      step(2)=step(1)
      stopcr=1.0d-6
      simp=1.0d-6
      iquad=1
      iprint=1
      nloop=10
      max=50000
c
c       Now open data files and output files
c
      open(7,file='ch4mm.dat')
      open(8,file='output.txt')
      open(2,file='resids.dat')
      open(4,file='logfile.txt')
c
c
      write(6,401)
      write(4,401)
c
c      Now read in data and then close input channel
c
      i=0
  100 i=i+1
      read(7,*,end=200)s(i),y(i),type(i)
      goto 100
  200 n=i-1
      close(7)
      write(4,402)n
      do 111 i=1,n
      write(4,403)s(i),y(i),type(i)
  111 continue
      write(4,404)
      write(6,701)
      write(6,702)
      read(5,*)b0
      write(6,703)
      read(5,*)k0
      write(4,704)b0,k0
c
c      Now successively call MINIM to mimimise the negative profile
c      log-likelihood for models I,III and II
c
c      (done in that order so that the last fitting is of model II)
c
      p(1)=b0
      p(2)=k0
      call MINIM(P,STEP,NOP,N2,FUNCM1,MAX,IPRINT,STOPCR,NLOOP,
     1  IQUAD,SIMP,VAR,loglk1,IFAULT)
      if(ifault.ne.0)write(6,811)ifault
c      write(6,801)p(1),p(2),-funcm1,ifault
      p(1)=b0
      p(2)=k0
      call MINIM(P,STEP,NOP,N2,FUNCM3,MAX,IPRINT,STOPCR,NLOOP,
     1  IQUAD,SIMP,VAR,loglk3,IFAULT)
      if(ifault.ne.0)write(6,813)ifault
c      write(6,803)p(1),p(2),-funcm3,ifault
      p(1)=b0
      p(2)=k0
      call MINIM(P,STEP,NOP,N2,FUNCM2,MAX,IPRINT,STOPCR,NLOOP,
     1  IQUAD,SIMP,VAR,loglk2,IFAULT)
      if(ifault.ne.0)write(6,812)ifault
      b=p(1)
      k=p(2)
c      write(6,802)b,k,-funcm2,ifault
c
c      Now calculate hypothesis tests for testing model I vs model II
c      and model II versus model III
c
      lrt1v2=2.0d0*(funcm1-funcm2)
      lrt3v2=2.0d0*(funcm2-funcm3)
      zchi=sqrt(lrt1v2)
      call NPROB(zchi,pleft,pright,pdf)
      pval12=2.0d0*pright
      zchi=sqrt(lrt3v2)
      call NPROB(zchi,pleft,pright,pdf)
      pval32=2.0d0*pright
      write(6,501)lrt1v2,pval12,lrt3v2,pval32
      write(4,501)lrt1v2,pval12,lrt3v2,pval32
      pause
c
c
c      Having got B and K for Model II, now get other parameters
c      starting with A
c
      a=zero
      nc12=0
      do 300 i=1,n
      if(type(i).eq.12)then
       nc12=nc12+1
       a=a+y(i)+b*s(i)/(k+s(i))
      end if
  300 continue
      a=a/dble(nc12)
c
c     Now get the residual variances for the 12C and 13C data
c
c
      s12=zero
      s13=zero
      do 301 i=1,n
      if(type(i).eq.12)then
       s12=s12+(y(i)-a+b*s(i)/(k+s(i)))**2
      else if(type(i).eq.13)then
       s13=s13+(y(i)-b*s(i)/(k+s(i)))**2
      end if
  301 continue
      s12=s12/dble(nc12)
      s13=s13/dble(n-nc12)
c
c     Now write out the residuals to a file for detailed model
c     checking
c
      do 250 i=1,n
      if(type(i).eq.12)then
	write(2,48)y(i)-a+b*s(i)/(k+s(i)),s(i),y(i),type(i)
      else if(type(i).eq.13)then
	write(2,48)y(i)-b*s(i)/(k+s(i)),s(i),y(i),type(i)
      end if
  250 continue
      close(2)
      call getvar(var,a,b,k,s12,s13)
      write(6,47)a,sqrt(var(1)),b,sqrt(var(3)),k,sqrt(var(6)),s12,s13
      write(4,47)a,sqrt(var(1)),b,sqrt(var(3)),k,sqrt(var(6)),s12,s13
      ratab=b/a*100.0d0
      serat=var(3)/(a*a)+(b*b*var(1))/(a**4)-2.0d0*b*var(2)/(a**3)
      serat=100.0d0*sqrt(serat)
      write(6,46)ratab,serat
      write(4,46)ratab,serat
      close(4)
      pause
      stop
   48 format(1x,3f8.4,i4)
   47 format(/,/,1x,'**** Results from fitting Model II ****',/,/,
     +1x,'Parameter    Estimate    Standard error',/,/,
     +1x,'   A     ',f12.4,f15.4,/,
     +1x,'   B     ',f12.4,f15.4,/,
     +1x,'   K     ',f12.4,f15.4,/,/,
     +1x,' Residual variance for 12C data is ',f18.8,/
     +1x,' Residual variance for 13C data is ',f18.8,/,/,/)
   46 format(/,/,1x,'Percentage of contribution of acetoclastic',/,
     +1x,'methanogenesis to total methane production is 100B/A',/,
     +1x,'This has value ',f6.1,'% and standard error ',f6.1,'%')
   45 format(1x,'Twice log-likelihood',f10.5)
   44 format(1x,'Twice log likelihood ratio',f10.5)
  501 format(/,/,1x,'**** Likelihood ratio tests ****',/,
     +1x,'Test for equal residual variances:',/,
     +1x,'likelihood ratio statistic = ',f15.4,' P-value = ',f8.4,/,/,
     +1x,'Test for non-constant total methane production:',/,
     +1x,'likelihood ratio statistic = ',f15.4,' P-value = ',f8.4,/)
  701 format(1x,'The fitting procedure requires the user to supply ',
     +/,1x,'starting values for the numerical method.',/,/,
     +1x,'Please see the relevant document on the website for',/,
     +1x,'guidance on choosing these.',/,/,
     +1x,'For the test data try starting values of 0.15 (B) and 0.05 (K)
     +')
  702 format(1x,'Please enter starting value for B  ',/)
  703 format(/,1x,'Please enter starting value for K  ',/)
  704 format(/,1x,'Starting value for B = ',f15.6,/,
     +1x,'Starting value for K = ',f15.6,/)
  801 format(/,/,1x,'Parameter values for model 1 are',/,
     +'in the order B, K, Profile log-lik. and ifault',/
     +,2f10.5,f20.8,i4)
  802 format(/,/,1x,'Parameter values for model 2 are',/,
     +'in the order B, K, Profile log-lik. and ifault',/
     +,2f10.5,f20.8,i4)
  803 format(/,/,1x,'Parameter values for model 3 are',/,
     +'in the order B, K, Profile log-lik. and ifault',/
     +,2f10.5,f20.8,i4)
  811 format(/,/,1x,'When fitting model 1, a non-zero ifault (=',i3,
     +') was obtained: results may be unreliable')
  812 format(/,/,1x,'When fitting model 2, a non-zero ifault (=',i3,
     +') was obtained: results may be unreliable')
  813 format(/,/,1x,'When fitting model 3, a non-zero ifault (=',i3,
     +') was obtained: results may be unreliable')
  401 format(1x,'**********   MMLINK   **********',/,/,
     +1x,'A program to fit linked equations of Michaelis-Menten type',/,
     +1x,'to data on methane production',/,/,/)
  402 format(1x,'Analysis is based on ',i3,' data points',/,/,
     +1x,'These are listed below in order:',
     +1x,'substrate concentration, CH4 production rate and carbon type',
     +/,/)
  403 format(1x,2f20.6,i4)
  404 format(1x,/,/)
      end

      subroutine getvar(var,a,b,k,s12,s13)
c
c      This routine calculates the expected information matrix
c      for A, B and K from model II
c      The estimates of the parameters are in a,b,k and
c      s12, s13 are the residual variances (*not* SDs)
c
      parameter(nmax=80,npar2=3,n2max=npar2*(npar2+1)/2)
      implicit double precision (a-h,o-z)
      double precision s(nmax),y(nmax),var(n2max),zero,g,
     +wk2(npar2),rmax,a,b,k,vart(n2max,2),one,s12,s13
      integer n,type(nmax)
c
c      vart(i,j) contains the ith element of the information matrix
c      based on the 12C data (j=1) or 13C data (j=2), held in lower
c      triangular form (but for constants)
c      var(i) contains the information matrix and is then inverted
c      to give the dispersion matrix in lower triangular form.
c
      common s,y,type,n
      data zero,one/0.0d0,1.0d0/
      nop=npar2
      do 1 j=1,n2max
      vart(j,1)=zero
      vart(j,2)=zero
    1 var(j)=zero
      do 100 i=1,n
      g=s(i)/(k+s(i))
      if(type(i).eq.12)then
       vart(1,1)=vart(1,1)+1.0d0
       vart(2,1)=vart(2,1)+g
       vart(3,1)=vart(3,1)+g*g
       vart(4,1)=vart(4,1)+g*(one-g)
       vart(5,1)=vart(5,1)+g*g*(one-g)
       vart(6,1)=vart(6,1)+g*g*(one-g)*(one-g)
      else if(type(i).eq.13)then
       vart(3,2)=vart(3,2)+g*g
       vart(5,2)=vart(5,2)+g*g*(one-g)
       vart(6,2)=vart(6,2)+g*g*(one-g)*(one-g)
      end if
  100 continue
      var(1)=vart(1,1)/s12
      var(2)=-vart(2,1)/s12
      var(3)=vart(3,1)/s12+vart(3,2)/s13
      var(4)=vart(4,1)*b/(k*s12)
      var(5)=-(b/k)*(vart(5,1)/s12+vart(5,2)/s13)
      var(6)=(vart(6,1)/s12+vart(6,2)/s13)*(b/k)**2
      call SYMINV(var,nop,var,wk2,NULLTY,IFAULT,RMAX)
c      write(6,49)NULLTY
      if(nullty.ne.0)then
	write(6,49)nullty
	close(4)
	pause
	stop
      end if
   49 format(1x,'Nullity is ',i2, ' which is not zero so information',/,
     +1x,'matrix is singular and standard errors are not available.',/,
     +1x,'This should only happen for ill-conditioned or, more likely',/
     +,1x,'small datasets.  Too few points for one of the types of ',/,
     +1x,'methane production will make data set small in this sense.',/,
     +1x,'The program now stops.')
      return
      end

      subroutine loglk1(p,func)
c
c     This routine calculates the negative of the profile log-likelihood
c     of B and K for model I (equal variances)
c
      parameter(nmax=80,npar=2)
      double precision p(npar),func,s(nmax),y(nmax),wk,ybar,fbar,
     + zero,b,k,a,func12,func13
      integer type(nmax),n
      common  s,y,type,n
      data zero/0.0d0/
      b=p(1)
      k=p(2)
      func=zero
      fbar=zero
      ybar=zero
      nc12=0
      do 50 i=1,n
      if(type(i).eq.12) then
	nc12=nc12+1
	ybar=ybar+y(i)
	fbar=fbar+b*s(i)/(s(i)+k)
      end if
   50 continue
      a=(ybar+fbar)/dble(nc12)
      func12=zero
      func13=zero
      do 100 i=1,n
      wk=b*s(i)/(k+s(i))
      if(type(i).eq.12) then
       func12=func12+(y(i)-a+wk)**2
      else if (type(i).eq.13) then
       func13=func13+(y(i)-wk)**2
      end if
  100 continue
      func=0.5d0*dble(n)*log(func12+func13)
c      write(6,49)a,b,k,sqrt((func12+func13)/dble(n))
      return
   49 format(1x,'a= ',f8.4,'b= ',f8.4,'k= ',f8.4,'sigma ',f8.4)
      end

      subroutine loglk2(p,func)
c
c      This routine calculates the negative of the profile
c      log-likelihood for B and K
c      for model II: different variances and constant total methane
c      production
c
      parameter(nmax=80,npar=2)
      double precision p(npar),func,s(nmax),y(nmax),wk,ybar,fbar,
     + zero,b,k,a,func12,func13
      integer type(nmax),n
      common  s,y,type,n
      data zero/0.0d0/
      b=p(1)
      k=p(2)
      func=zero
      fbar=zero
      ybar=zero
      nc12=0
      do 50 i=1,n
      if(type(i).eq.12) then
	nc12=nc12+1
	ybar=ybar+y(i)
	fbar=fbar+b*s(i)/(s(i)+k)
      end if
   50 continue
      a=(ybar+fbar)/dble(nc12)
      func12=zero
      func13=zero
      do 100 i=1,n
      wk=b*s(i)/(k+s(i))
      if(type(i).eq.12) then
       func12=func12+(y(i)-a+wk)**2
      else if (type(i).eq.13) then
       func13=func13+(y(i)-wk)**2
      end if
  100 continue
      func=0.5d0*(dble(nc12)*log(func12)+dble(n-nc12)*log(func13))
      return
      end


      subroutine loglk3(p,func)
      parameter(nmax=80,npar=2)
c
c        This routine computes the log-likelihood for model III
c        the model with non-constant total methane production
c
c
      double precision p(npar),func,s(nmax),y(nmax),wk,ybar,fbar,
     + zero,b,k,a,func12,func13,sbar,s2bar,fsbar,ysbar,f,rn,c
      integer type(nmax),n
      common  s,y,type,n
      data zero/0.0d0/
      b=p(1)
      k=p(2)
      func=zero
      fbar=zero
      ybar=zero
      sbar=zero
      s2bar=zero
      fsbar=zero
      ysbar=zero
      nc12=0
      do 50 i=1,n
      if(type(i).eq.12) then
	f=b*s(i)/(k+s(i))
	nc12=nc12+1
	ybar=ybar+y(i)
	fbar=fbar+f
	sbar=sbar+s(i)
	s2bar=s2bar+s(i)**2
	fsbar=fsbar+f*s(i)
	ysbar=ysbar+y(i)*s(i)
      end if
   50 continue
      rn=dble(nc12)
      ybar=ybar/rn
      fbar=fbar/rn
      sbar=sbar/rn
      s2bar=s2bar/rn
      fsbar=fsbar/rn
      ysbar=ysbar/rn
      c=(ysbar+fsbar-ybar*sbar-fbar*sbar)/(s2bar-sbar**2)
      a=ybar+fbar-c*sbar
      func12=zero
      func13=zero
      do 100 i=1,n
      wk=b*s(i)/(k+s(i))
      if(type(i).eq.12) then
       func12=func12+(y(i)-a-c*s(i)+wk)**2
      else if (type(i).eq.13) then
       func13=func13+(y(i)-wk)**2
      end if
  100 continue
      func=0.5d0*(dble(nc12)*log(func12)+dble(n-nc12)*log(func13))
c      write(6,49)a,c,b,k,sqrt(func12/dble(nc12)),sqrt(func13/
c     +dble(n-nc12))
      return
   49 format(1x,'a= ',f8.5,' c= ',f8.5,' b= ',f8.5,' k= ',f8.5,
     + 's12= ',f8.5,'s13= ',f8.5)
      end
c
	SUBROUTINE MINIM(P,STEP,NOP,N2,FUNC,MAX,IPRINT,STOPCR,NLOOP,
     1  IQUAD,SIMP,VAR,FUNCTN,IFAULT)
C
C       A PROGRAM FOR FUNCTION MINIMIZATION USING THE SIMPLEX METHOD.
C       The minimum found will often be a local, not a global, minimum.
C
C       FOR DETAILS, SEE NELDER & MEAD, THE COMPUTER JOURNAL, JANUARY 1965
C
C       PROGRAMMED BY D.E.SHAW,
C       CSIRO, DIVISION OF MATHEMATICS & STATISTICS
C       P.O. BOX 218, LINDFIELD, N.S.W. 2070
C
C       WITH AMENDMENTS BY R.W.M.WEDDERBURN
C       ROTHAMSTED EXPERIMENTAL STATION
C       HARPENDEN, HERTFORDSHIRE, ENGLAND
C
C       Further amended by Alan Miller,
C       CSIRO, Division of Mathematics & Statistics
C       Private Bag 10, CLAYTON, VIC. 3168
C
C       ARGUMENTS:-
C       P()     = INPUT, STARTING VALUES OF PARAMETERS
C                 OUTPUT, FINAL VALUES OF PARAMETERS
C       STEP()  = INPUT, INITIAL STEP SIZES
C       NOP     = INPUT, NO. OF PARAMETERS, INCL. ANY TO BE HELD FIXED
C       FUNC    = OUTPUT, THE FUNCTION VALUE CORRESPONDING TO THE FINAL
C                 PARAMETER VALUES
C       MAX     = INPUT, THE MAXIMUM NO. OF FUNCTION EVALUATIONS ALLOWED
C       IPRINT  = INPUT, PRINT CONTROL PARAMETER
C                       < 0 NO PRINTING
C                       = 0 PRINTING OF PARAMETER VALUES AND THE FUNCTION
C                           VALUE AFTER INITIAL EVIDENCE OF CONVERGENCE.
C                       > 0 AS FOR IPRINT = 0 PLUS PROGRESS REPORTS AFTER
C                           EVERY IPRINT EVALUATIONS, PLUS PRINTING FOR THE
C                           INITIAL SIMPLEX.
C       STOPCR  = INPUT, STOPPING CRITERION
C       NLOOP   = INPUT, THE STOPPING RULE IS APPLIED AFTER EVERY NLOOP
C                 FUNCTION EVALUATIONS.
C       IQUAD   = INPUT, = 1 IF THE FITTING OF A QUADRATIC SURFACE IS REQUIRED
C                        = 0 IF NOT
C       SIMP    = INPUT, CRITERION FOR EXPANDING THE SIMPLEX TO OVERCOME
C                 ROUNDING ERRORS BEFORE FITTING THE QUADRATIC SURFACE.
C       VAR()   = OUTPUT, CONTAINS THE DIAGONAL ELEMENTS OF THE INVERSE OF
C                 THE INFORMATION MATRIX.
C       FUNCTN  = INPUT, NAME OF THE USER'S SUBROUTINE - ARGUMENTS (P,FUNC)
C                 WHICH RETURNS THE FUNCTION VALUE FOR A GIVEN SET OF
C                 PARAMETER VALUES IN ARRAY P.
C****   FUNCTN MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM.
C       IFAULT  = OUTPUT, = 0 FOR SUCCESSFUL TERMINATION
C                         = 1 IF MAXIMUM NO. OF FUNCTION EVALUATIONS EXCEEDED
C                         = 2 IF INFORMATION MATRIX IS NOT +VE SEMI-DEFINITE
C                         = 3 IF NOP < 1
C                         = 4 IF NLOOP < 1
C
C       Advice on usage:
C       If the function minimized can be expected to be smooth in the vicinity
C       of the minimum, users are strongly urged to use the quadratic-surface
C       fitting option.   This is the only satisfactory way of testing that the
C       minimum has been found.   The value of SIMP should be set to at least
C       1000 times the rounding error in calculating the fitted function.
C       e.g. in double precision on a micro- or mini-computer with about 16
C       decimal digit representation of floating-point numbers, the rounding
C       errors in calculating the objective function may be of the order of
C       1.E-12 say in a particular case.   A suitable value for SIMP would then
C       be 1.E-08.   However, if numerical integration is required in the
C       calculation of the objective function, it may only be accurate to say
C       1.E-05 and an appropriate value for SIMP would be about 0.1.
C       If the fitted quadratic surface is not +ve definite (and the function
C       should be smooth in the vicinity of the minimum), it probably means
C       that the search terminated prematurely and you have not found the
C       minimum.
C
C       N.B. P, STEP AND VAR (IF IQUAD = 1) MUST HAVE DIMENSION AT LEAST NOP
C            IN THE CALLING PROGRAM.
C       THE DIMENSIONS BELOW ARE FOR A MAXIMUM OF 20 PARAMETERS.
C       The dimension of BMAT should be at least NOP*(NOP+1)/2.
C
C****   N.B. This version is in DOUBLE PRECISION throughout
C
C       LATEST REVISION - 18 October 1985
C
C*****************************************************************************
C
	implicit double precision (a-h, o-z)
	external FUNCTN
	DIMENSION P(NOP),STEP(NOP),VAR(N2)
	DIMENSION G(21,20),H(21),PBAR(20),PSTAR(20),PSTST(20),AVAL(20),
     1  BMAT(210),PMIN(20),VC(210),TEMP(20)
	DATA ZERO/0.D0/, ONE/1.D0/, TWO/2.D0/, THREE/3.D0/, HALF/0.5D0/
C
C       A = REFLECTION COEFFICIENT, B = CONTRACTION COEFFICIENT, AND
C       C = EXPANSION COEFFICIENT.
C
	DATA A,B,C/1.D0, 0.5D0, 2.D0/
C
C       SET LOUT = LOGICAL UNIT NO. FOR OUTPUT
C
	DATA LOUT/8/
C
C       IF PROGRESS REPORTS HAVE BEEN REQUESTED, PRINT HEADING
C
	IF(IPRINT.GT.0) WRITE(LOUT,1000) IPRINT
 1000   FORMAT(' PROGRESS REPORT EVERY',I4,' FUNCTION EVALUATIONS'/,
     1  ' EVAL.  FUNC.',15X,'PARAMETER VALUES')
C
C       CHECK INPUT ARGUMENTS
C
	IFAULT=0
	IF(NOP.LE.0) IFAULT=3
	IF(NLOOP.LE.0) IFAULT=4
	IF(IFAULT.NE.0) RETURN
C
C       SET NAP = NO. OF PARAMETERS TO BE VARIED, I.E. WITH STEP.NE.0
C
	NAP=0
	LOOP=0
	IFLAG=0
	DO 10 I=1,NOP
	  IF(STEP(I).NE.ZERO) NAP=NAP+1
   10   CONTINUE
C
C       IF NAP = 0 EVALUATE FUNCTION AT THE STARTING POINT AND RETURN
C
	IF(NAP.GT.0) GO TO 30
	CALL FUNCTN(P,FUNC)
	RETURN
C
C       SET UP THE INITIAL SIMPLEX
C
   30   DO 40 I=1,NOP
   40   G(1,I)=P(I)
	IROW=2
	DO 60 I=1,NOP
	  IF(STEP(I).EQ.ZERO) GO TO 60
	  DO 50 J=1,NOP
   50     G(IROW,J)=P(J)
	  G(IROW,I)=P(I)+STEP(I)
	  IROW=IROW+1
   60   CONTINUE
	NP1=NAP+1
	NEVAL=0
	DO 90 I=1,NP1
	  DO 70 J=1,NOP
   70     P(J)=G(I,J)
	  CALL FUNCTN(P,H(I))
	  NEVAL=NEVAL+1
	  IF(IPRINT.LE.0) GO TO 90
	  WRITE(LOUT,1010) NEVAL,H(I),(P(J),J=1,NOP)
 1010     FORMAT(/I4, 2X, G12.5, 2X, 5G12.5, 3(/20X, 5G12.5))
   90   CONTINUE
C
C       START OF MAIN CYCLE.
C
C       FIND MAX. & MIN. VALUES FOR CURRENT SIMPLEX (HMAX & HMIN).
C
  100   LOOP=LOOP+1
	IMAX=1
	IMIN=1
	HMAX=H(1)
	HMIN=H(1)
	DO 120 I=2,NP1
	  IF(H(I).LE.HMAX) GO TO 110
	  IMAX=I
	  HMAX=H(I)
	  GO TO 120
  110     IF(H(I).GE.HMIN) GO TO 120
	  IMIN=I
	  HMIN=H(I)
  120   CONTINUE
C
C       FIND THE CENTROID OF THE VERTICES OTHER THAN P(IMAX)
C
	DO 130 I=1,NOP
  130   PBAR(I)=ZERO
	DO 150 I=1,NP1
	  IF(I.EQ.IMAX) GO TO 150
	  DO 140 J=1,NOP
  140     PBAR(J)=PBAR(J)+G(I,J)
  150   CONTINUE
	DO 160 J=1,NOP
	  FNAP = NAP
  160   PBAR(J)=PBAR(J)/FNAP
C
C       REFLECT MAXIMUM THROUGH PBAR TO PSTAR,
C       HSTAR = FUNCTION VALUE AT PSTAR.
C
	DO 170 I=1,NOP
  170   PSTAR(I)=A*(PBAR(I)-G(IMAX,I))+PBAR(I)
	CALL FUNCTN(PSTAR,HSTAR)
	NEVAL=NEVAL+1
	IF(IPRINT.LE.0) GO TO 180
	IF(MOD(NEVAL,IPRINT).EQ.0) WRITE(LOUT,1010) NEVAL,HSTAR,
     1  (PSTAR(J),J=1,NOP)
C
C       IF HSTAR < HMIN, REFLECT PBAR THROUGH PSTAR,
C       HSTST = FUNCTION VALUE AT PSTST.
C
  180   IF(HSTAR.GE.HMIN) GO TO 220
	DO 190 I=1,NOP
  190   PSTST(I)=C*(PSTAR(I)-PBAR(I))+PBAR(I)
	CALL FUNCTN(PSTST,HSTST)
	NEVAL=NEVAL+1
	IF(IPRINT.LE.0) GO TO 200
	IF(MOD(NEVAL,IPRINT).EQ.0) WRITE(LOUT,1010) NEVAL,HSTST,
     1  (PSTST(J),J=1,NOP)
C
C       IF HSTST < HMIN REPLACE CURRENT MAXIMUM POINT BY PSTST AND
C       HMAX BY HSTST, THEN TEST FOR CONVERGENCE.
C
  200   IF(HSTST.GE.HMIN) GO TO 320
	DO 210 I=1,NOP
	  IF(STEP(I).NE.ZERO) G(IMAX,I)=PSTST(I)
  210   CONTINUE
	H(IMAX)=HSTST
	GO TO 340
C
C       HSTAR IS NOT < HMIN.
C       TEST WHETHER IT IS < FUNCTION VALUE AT SOME POINT OTHER THAN
C       P(IMAX).   IF IT IS REPLACE P(IMAX) BY PSTAR & HMAX BY HSTAR.
C
  220   DO 230 I=1,NP1
	  IF(I.EQ.IMAX) GO TO 230
	  IF(HSTAR.LT.H(I)) GO TO 320
  230   CONTINUE
C
C       HSTAR > ALL FUNCTION VALUES EXCEPT POSSIBLY HMAX.
C       IF HSTAR <= HMAX, REPLACE P(IMAX) BY PSTAR & HMAX BY HSTAR.
C
	IF(HSTAR.GT.HMAX) GO TO 260
	DO 250 I=1,NOP
	  IF(STEP(I).NE.ZERO) G(IMAX,I)=PSTAR(I)
  250   CONTINUE
	HMAX=HSTAR
	H(IMAX)=HSTAR
C
C       CONTRACTED STEP TO THE POINT PSTST,
C       HSTST = FUNCTION VALUE AT PSTST.
C
  260   DO 270 I=1,NOP
  270   PSTST(I)=B*G(IMAX,I) + (1.d0-B)*PBAR(I)
	CALL FUNCTN(PSTST,HSTST)
	NEVAL=NEVAL+1
	IF(IPRINT.LE.0) GO TO 280
	IF(MOD(NEVAL,IPRINT).EQ.0) WRITE(LOUT,1010) NEVAL,HSTST,
     1  (PSTST(J),J=1,NOP)
C
C       IF HSTST < HMAX REPLACE P(IMAX) BY PSTST & HMAX BY HSTST.
C
  280   IF(HSTST.GT.HMAX) GO TO 300
	DO 290 I=1,NOP
	  IF(STEP(I).NE.ZERO) G(IMAX,I)=PSTST(I)
  290   CONTINUE
	H(IMAX)=HSTST
	GO TO 340
C
C       HSTST > HMAX.
C       SHRINK THE SIMPLEX BY REPLACING EACH POINT, OTHER THAN THE CURRENT
C       MINIMUM, BY A POINT MID-WAY BETWEEN ITS CURRENT POSITION AND THE
C       MINIMUM.
C
  300   DO 315 I=1,NP1
	  IF(I.EQ.IMIN) GO TO 315
	  DO 310 J=1,NOP
	    IF(STEP(J).NE.ZERO) G(I,J)=(G(I,J)+G(IMIN,J))*HALF
	    P(J)=G(I,J)
  310     CONTINUE
	  CALL FUNCTN(P,H(I))
	  NEVAL=NEVAL+1
	  IF(IPRINT.LE.0) GO TO 315
	  IF(MOD(NEVAL,IPRINT).EQ.0) WRITE(LOUT,1010) NEVAL,H(I),
     1                (P(J),J=1,NOP)
  315   CONTINUE
	GO TO 340
C
C       REPLACE MAXIMUM POINT BY PSTAR & H(IMAX) BY HSTAR.
C
  320   DO 330 I=1,NOP
	  IF(STEP(I).NE.ZERO) G(IMAX,I)=PSTAR(I)
  330   CONTINUE
	H(IMAX)=HSTAR
C
C       IF LOOP = NLOOP TEST FOR CONVERGENCE, OTHERWISE REPEAT MAIN CYCLE.
C
  340   IF(LOOP.LT.NLOOP) GO TO 100
C
C       CALCULATE MEAN & STANDARD DEVIATION OF FUNCTION VALUES FOR THE
C       CURRENT SIMPLEX.
C
	HSTD=ZERO
	HMEAN=ZERO
	DO 350 I=1,NP1
  350   HMEAN=HMEAN+H(I)
	FNP1 = NP1
	HMEAN=HMEAN/FNP1
	DO 360 I=1,NP1
  360   HSTD=HSTD+(H(I)-HMEAN)**2
	HSTD=SQRT(HSTD/FLOAT(NP1))
C
C       IF THE RMS > STOPCR, SET IFLAG & LOOP TO ZERO AND GO TO THE
C       START OF THE MAIN CYCLE AGAIN.
C
	IF(HSTD.LE.STOPCR.OR.NEVAL.GT.MAX) GO TO 410
	IFLAG=0
	LOOP=0
	GO TO 100
C
C       FIND THE CENTROID OF THE CURRENT SIMPLEX AND THE FUNCTION VALUE THERE.
C
  410   DO 380 I=1,NOP
	  IF(STEP(I).EQ.ZERO) GO TO 380
	  P(I)=ZERO
	  DO 370 J=1,NP1
  370     P(I)=P(I)+G(J,I)
	  FNP1 = NP1
	  P(I)=P(I)/FNP1
  380   CONTINUE
	CALL FUNCTN(P,FUNC)
	NEVAL=NEVAL+1
	IF(IPRINT.LE.0) GO TO 390
	IF(MOD(NEVAL,IPRINT).EQ.0) WRITE(LOUT,1010) NEVAL,FUNC,
     1  (P(J),J=1,NOP)
C
C       TEST WHETHER THE NO. OF FUNCTION VALUES ALLOWED, MAX, HAS BEEN
C       OVERRUN       IF SO, EXIT WITH IFAULT = 1.
C
  390   IF(NEVAL.LE.MAX) GO TO 420
	IFAULT=1
	IF(IPRINT.LT.0) RETURN
	WRITE(LOUT,1020) MAX
 1020   FORMAT(' NO. OF FUNCTION EVALUATIONS EXCEEDS',I5)
	WRITE(LOUT,1030) HSTD
 1030   FORMAT(' RMS OF FUNCTION VALUES OF LAST SIMPLEX =',G14.6)
	WRITE(LOUT,1040)(P(I),I=1,NOP)
 1040   FORMAT(' CENTROID OF LAST SIMPLEX =',4(/1X,6G13.5))
	WRITE(LOUT,1050) FUNC
 1050   FORMAT(' FUNCTION VALUE AT CENTROID =',G14.6)
	RETURN
C
C       CONVERGENCE CRITERION SATISFIED.
C       IF IFLAG = 0, SET IFLAG & SAVE HMEAN.
C       IF IFLAG = 1 & CHANGE IN HMEAN <= STOPCR THEN SEARCH IS COMPLETE.
C
  420   IF(IPRINT.LT.0) GO TO 430
	WRITE(LOUT,1060)
 1060   FORMAT(/' EVIDENCE OF CONVERGENCE')
	WRITE(LOUT,1040)(P(I),I=1,NOP)
	WRITE(LOUT,1050) FUNC
  430   IF(IFLAG.GT.0) GO TO 450
	IFLAG=1
  440   SAVEMN=HMEAN
	LOOP=0
	GO TO 100
  450   IF(ABS(SAVEMN-HMEAN).GE.STOPCR) GO TO 440
	IF(IPRINT.LT.0) GO TO 460
	WRITE(LOUT,1070) NEVAL
 1070   FORMAT(//' MINIMUM FOUND AFTER',I5,' FUNCTION EVALUATIONS')
	WRITE(LOUT,1080)(P(I),I=1,NOP)
 1080   FORMAT(' MINIMUM AT',4(/1X,6G13.6))
	WRITE(LOUT,1090) FUNC
 1090   FORMAT(' FUNCTION VALUE AT MINIMUM =',G14.6)
  460   IF(IQUAD.LE.0) RETURN
C-------------------------------------------------------------------
C
C       QUADRATIC SURFACE FITTING
C
	IF(IPRINT.GE.0) WRITE(LOUT,1110)
 1110   FORMAT(/' QUADRATIC SURFACE FITTING ABOUT SUPPOSED MINIMUM'/)
C
C       EXPAND THE FINAL SIMPLEX, IF NECESSARY, TO OVERCOME ROUNDING
C       ERRORS.
C
	NEVAL=0
	DO 490 I=1,NP1
  470     TEST=ABS(H(I)-FUNC)
	  IF(TEST.GE.SIMP) GO TO 490
	  DO 480 J=1,NOP
	    IF(STEP(J).NE.ZERO) G(I,J)=(G(I,J)-P(J))+G(I,J)
	    PSTST(J)=G(I,J)
  480     CONTINUE
	  CALL FUNCTN(PSTST,H(I))
	  NEVAL=NEVAL+1
	  GO TO 470
  490   CONTINUE
C
C       FUNCTION VALUES ARE CALCULATED AT AN ADDITIONAL NAP POINTS.
C
	DO 510 I=1,NAP
	  I1=I+1
	  DO 500 J=1,NOP
  500     PSTAR(J)=(G(1,J)+G(I1,J))*HALF
	  CALL FUNCTN(PSTAR,AVAL(I))
	  NEVAL=NEVAL+1
  510   CONTINUE
C
C       THE MATRIX OF ESTIMATED SECOND DERIVATIVES IS CALCULATED AND ITS
C       LOWER TRIANGLE STORED IN BMAT.
C
	A0=H(1)
	DO 540 I=1,NAP
	  I1=I-1
	  I2=I+1
	  IF(I1.LT.1) GO TO 540
	  DO 530 J=1,I1
	    J1=J+1
	    DO 520 K=1,NOP
  520       PSTST(K)=(G(I2,K)+G(J1,K))*HALF
	    CALL FUNCTN(PSTST,HSTST)
	    NEVAL=NEVAL+1
	    L=I*(I-1)/2+J
	    BMAT(L)=TWO*(HSTST+A0-AVAL(I)-AVAL(J))
  530     CONTINUE
  540   CONTINUE
	L=0
	DO 550 I=1,NAP
	  I1=I+1
	  L=L+I
	  BMAT(L)=TWO*(H(I1)+A0-TWO*AVAL(I))
  550   CONTINUE
C
C       THE VECTOR OF ESTIMATED FIRST DERIVATIVES IS CALCULATED AND
C       STORED IN AVAL.
C
	DO 560 I=1,NAP
	  I1=I+1
	  AVAL(I)=TWO*AVAL(I)-(H(I1)+THREE*A0)*HALF
  560   CONTINUE
C
C       THE MATRIX Q OF NELDER & MEAD IS CALCULATED AND STORED IN G.
C
	DO 570 I=1,NOP
  570   PMIN(I)=G(1,I)
	DO 580 I=1,NAP
	  I1=I+1
	  DO 580 J=1,NOP
	  G(I1,J)=G(I1,J)-G(1,J)
  580   CONTINUE
	DO 590 I=1,NAP
	  I1=I+1
	  DO 590 J=1,NOP
	    G(I,J)=G(I1,J)
  590   CONTINUE
C
C       INVERT BMAT
C
	CALL SYMINV(BMAT,NAP,BMAT,TEMP,NULLTY,IFAULT,RMAX)
	IF(IFAULT.NE.0) GO TO 600
	IRANK=NAP-NULLTY
	GO TO 610
  600   IF(IPRINT.GE.0) WRITE(LOUT,1120)
 1120   FORMAT(/' MATRIX OF ESTIMATED SECOND DERIVATIVES NOT +VE DEFN.'/
     1  ' MINIMUM PROBABLY NOT FOUND'/)
	IFAULT=2
	RETURN
C
C       BMAT*A/2 IS CALCULATED AND STORED IN H.
C
  610   DO 650 I=1,NAP
	  H(I)=ZERO
	  DO 640 J=1,NAP
	    IF(J.GT.I) GO TO 620
	    L=I*(I-1)/2+J
	    GO TO 630
  620       L=J*(J-1)/2+I
  630       H(I)=H(I)+BMAT(L)*AVAL(J)
  640     CONTINUE
  650   CONTINUE
C
C       FIND THE POSITION, PMIN, & VALUE, YMIN, OF THE MINIMUM OF THE
C       QUADRATIC.
C
	YMIN=ZERO
	DO 660 I=1,NAP
  660   YMIN=YMIN+H(I)*AVAL(I)
	YMIN=A0-YMIN
	DO 670 I=1,NOP
	  PSTST(I)=ZERO
	  DO 670 J=1,NAP
  670   PSTST(I)=PSTST(I)+H(J)*G(J,I)
	DO 680 I=1,NOP
  680   PMIN(I)=PMIN(I)-PSTST(I)
	IF(IPRINT.LT.0) GO TO 682
	WRITE(LOUT,1130) YMIN,(PMIN(I),I=1,NOP)
 1130   FORMAT(' MINIMUM OF QUADRATIC SURFACE =',G14.6,' AT',
     1  4(/1X,6G13.5))
	WRITE(LOUT,1150)
 1150   FORMAT(' IF THIS DIFFERS BY MUCH FROM THE MINIMUM ESTIMATED',
     1  1X,'FROM THE MINIMIZATION,'/
     2  ' THE MINIMUM MAY BE FALSE &/OR THE INFORMATION MATRIX MAY BE',
     3  1X,'INACCURATE'/)
c
c       Calculate true function value at the minimum of the quadratic.
c
  682   neval = neval + 1
	call functn(pmin, hstar)
c
c       If HSTAR < FUNC, replace search minimum with quadratic minimum.
c
	if (hstar .ge. func) go to 690
	func = hstar
	do 684 i = 1, nop
  684   p(i) = pmin(i)
	if(iprint.ge.0)write(lout, 1140) func
 1140   format(' True func. value at minimum of quadratic = ', g14.6/)
C
C       Q*BMAT*Q'/2 IS CALCULATED & ITS LOWER TRIANGLE STORED IN VC
C
  690   DO 760 I=1,NOP
	  DO 730 J=1,NAP
	    H(J)=ZERO
	    DO 720 K=1,NAP
	      IF(K.GT.J) GO TO 700
	      L=J*(J-1)/2+K
	      GO TO 710
  700         L=K*(K-1)/2+J
  710         H(J)=H(J)+BMAT(L)*G(K,I)*HALF
  720       CONTINUE
  730     CONTINUE
	  DO 750 J=I,NOP
	    L=J*(J-1)/2+I
	    VC(L)=ZERO
	    DO 740 K=1,NAP
  740       VC(L)=VC(L)+H(K)*G(K,J)
  750     CONTINUE
  760   CONTINUE
C
C       THE ELEMENTS OF VC ARE COPIED INTO VAR.
C
	J=0
	DO 770 I=1,N2
	  J=J+1
	  VAR(I)=VC(I)
  770   CONTINUE
	IF(IPRINT.LT.0) RETURN
	WRITE(LOUT,1160) IRANK
 1160   FORMAT(' RANK OF INFORMATION MATRIX =',I3/
     1  ' GENERALIZED INVERSE OF INFORMATION MATRIX:-')
	IJK=1
	GO TO 880
  790   CONTINUE
	WRITE(LOUT,1170)
 1170   FORMAT(/' IF THE FUNCTION MINIMIZED WAS -LOG(LIKELIHOOD),'/
     1  ' THIS IS THE COVARIANCE MATRIX OF THE PARAMETERS'/
     2  ' IF THE FUNCTION WAS A SUM OF SQUARES OF RESIDUALS'/
     3  ' THIS MATRIX MUST BE MULTIPLIED BY TWICE THE ESTIMATED',
     4  1X,'RESIDUAL VARIANCE'/' TO OBTAIN THE COVARIANCE MATRIX.'/)
	CALL SYMINV(VC,NAP,BMAT,TEMP,NULLTY,IFAULT,RMAX)
C
C       BMAT NOW CONTAINS THE INFORMATION MATRIX
C
	WRITE(LOUT,1190)
 1190   FORMAT(' INFORMATION MATRIX:-'/)
	IJK=3
	GO TO 880
c
c       Calculate correlations of parameter estimates, put into VC.
c
  800   IJK=2
	II=0
	IJ=0
	DO 840 I=1,NOP
	  II=II+I
	  IF(VC(II).LE.ZERO) GO TO 810
	  VC(II)=ONE/SQRT(VC(II))
	  GO TO 820
  810     VC(II)=ZERO
  820     IF(I.EQ.1) GO TO 840
	  IM1=I-1
	  JJ=0
	  DO 830 J=1,IM1
	    JJ=JJ+J
	    IJ=IJ+1
	    VC(IJ)=VC(IJ)*VC(II)*VC(JJ)
  830     CONTINUE
  840     IJ=IJ+1
	CONTINUE
	WRITE(LOUT,1200)
 1200   FORMAT(/' CORRELATION MATRIX:-')
	II=0
	DO 850 I=1,NOP
	  II=II+I
	  IF(VC(II).NE.ZERO) VC(II)=ONE
  850   CONTINUE
	GO TO 880
  860   WRITE(LOUT,1210) NEVAL
 1210   FORMAT(/' A FURTHER',I4,' FUNCTION EVALUATIONS HAVE BEEN USED'/)
	RETURN
c
c       Pseudo-subroutine to print VC if IJK = 1 or 2, or
c       BMAT if IJK = 3.
c
  880   L=1
  890   IF(L.GT.NOP) GO TO (790,860,800),IJK
	II=L*(L-1)/2
	DO 910 I=L,NOP
	  I1=II+L
	  II=II+I
	  I2=MIN(II,I1+5)
	  IF(IJK.EQ.3) GO TO 900
	  WRITE(LOUT,1230)(VC(J),J=I1,I2)
	  GO TO 910
  900     WRITE(LOUT,1230)(BMAT(J),J=I1,I2)
  910   CONTINUE
 1230   FORMAT(1X,6G13.5)
	WRITE(LOUT,1240)
 1240   FORMAT(/)
	L=L+6
	GO TO 890
	END




	SUBROUTINE SYMINV(A,N,C,W,NULLTY,IFAULT,RMAX)
C
C       ALGORITHM AS7, APPLIED STATISTICS, VOL.17, 1968.
C
C       ARGUMENTS:-
C       A()     = INPUT, THE SYMMETRIC MATRIX TO BE INVERTED, STORED IN
C                 LOWER TRIANGULAR FORM
C       N       = INPUT, ORDER OF THE MATRIX
C       C()     = OUTPUT, THE INVERSE OF A (A GENERALIZED INVERSE IF C IS
C                 SINGULAR), ALSO STORED IN LOWER TRIANGULAR.
C                 C AND A MAY OCCUPY THE SAME LOCATIONS.
C       W()     = WORKSPACE, DIMENSION AT LEAST N.
C       NULLTY  = OUTPUT, THE RANK DEFICIENCY OF A.
C       IFAULT  = OUTPUT, ERROR INDICATOR
C                       = 1 IF N < 1
C                       = 2 IF A IS NOT +VE SEMI-DEFINITE
C                       = 0 OTHERWISE
C       RMAX    = OUTPUT, APPROXIMATE BOUND ON THE ACCURACY OF THE DIAGONAL
C                 ELEMENTS OF C.  E.G. IF RMAX = 1.E-04 THEN THE DIAGONAL
C                 ELEMENTS OF C WILL BE ACCURATE TO ABOUT 4 DEC. DIGITS.
C
C       LATEST REVISION - 18 October 1985
C
C***************************************************************************
C
	implicit double precision (a-h, o-z)
	DIMENSION A(*),C(*),W(N)
	DATA ZERO/0.D0/, ONE/1.D0/
C
	NROW=N
	IFAULT=1
	IF(NROW.LE.0) GO TO 100
	IFAULT=0
C
C       CHOLESKY FACTORIZATION OF A, RESULT IN C
C
	CALL CHOLA(A,NROW,C,NULLTY,IFAULT,RMAX,W)
	IF(IFAULT.NE.0) GO TO 100
C
C       INVERT C & FORM THE PRODUCT (CINV)'*CINV, WHERE CINV IS THE INVERSE
C       OF C, ROW BY ROW STARTING WITH THE LAST ROW.
C       IROW = THE ROW NUMBER, NDIAG = LOCATION OF LAST ELEMENT IN THE ROW.
C
	NN=NROW*(NROW+1)/2
	IROW=NROW
	NDIAG=NN
   10   IF(C(NDIAG).EQ.ZERO) GO TO 60
	L=NDIAG
	DO 20 I=IROW,NROW
	  W(I)=C(L)
	  L=L+I
   20   CONTINUE
	ICOL=NROW
	JCOL=NN
	MDIAG=NN
   30   L=JCOL
	X=ZERO
	IF(ICOL.EQ.IROW) X=ONE/W(IROW)
	K=NROW
   40   IF(K.EQ.IROW) GO TO 50
	X=X-W(K)*C(L)
	K=K-1
	L=L-1
	IF(L.GT.MDIAG) L=L-K+1
	GO TO 40
   50   C(L)=X/W(IROW)
	IF(ICOL.EQ.IROW) GO TO 80
	MDIAG=MDIAG-ICOL
	ICOL=ICOL-1
	JCOL=JCOL-1
	GO TO 30
c
c       Special case, zero diagonal element.
c
   60   L=NDIAG
	DO 70 J=IROW,NROW
	  C(L)=ZERO
	  L=L+J
   70   CONTINUE
c
c       End of row.
c
   80   NDIAG=NDIAG-IROW
	IROW=IROW-1
	IF(IROW.NE.0) GO TO 10
  100   RETURN
	END





	SUBROUTINE CHOLA(A, N, U, NULLTY, IFAULT, RMAX, R)
C
C       ALGORITHM AS6, APPLIED STATISTICS, VOL.17, 1968, WITH
C       MODIFICATIONS BY A.J.MILLER
C
C       ARGUMENTS:-
C       A()     = INPUT, A +VE DEFINITE MATRIX STORED IN LOWER-TRIANGULAR
C                 FORM.
C       N       = INPUT, THE ORDER OF A
C       U()     = OUTPUT, A LOWER TRIANGULAR MATRIX SUCH THAT U*U' = A.
C                 A & U MAY OCCUPY THE SAME LOCATIONS.
C       NULLTY  = OUTPUT, THE RANK DEFICIENCY OF A.
C       IFAULT  = OUTPUT, ERROR INDICATOR
C                       = 1 IF N < 1
C                       = 2 IF A IS NOT +VE SEMI-DEFINITE
C                       = 0 OTHERWISE
C       RMAX    = OUTPUT, AN ESTIMATE OF THE RELATIVE ACCURACY OF THE
C                 DIAGONAL ELEMENTS OF U.
C       R()     = OUTPUT, ARRAY CONTAINING BOUNDS ON THE RELATIVE ACCURACY
C                 OF EACH DIAGONAL ELEMENT OF U.
C
C       LATEST REVISION - 18 October 1985
C
C***************************************************************************
C
	implicit double precision (a-h, o-z)
	DIMENSION A(*),U(*),R(N)
C
C       ETA SHOULD BE SET EQUAL TO THE SMALLEST +VE VALUE SUCH THAT
C       1.0 + ETA IS CALCULATED AS BEING GREATER THAN 1.0 IN THE ACCURACY
C       BEING USED.
C
	DATA ETA/1.D-16/, ZERO/0.D0/, FIVE/5.D0/
C
	IFAULT=1
	IF(N.LE.0) GO TO 100
	IFAULT=2
	NULLTY=0
	RMAX=ETA
	R(1)=ETA
	J=1
	K=0
C
C       FACTORIZE COLUMN BY COLUMN, ICOL = COLUMN NO.
C
	DO 80 ICOL=1,N
	  L=0
C
C       IROW = ROW NUMBER WITHIN COLUMN ICOL
C
	  DO 40 IROW=1,ICOL
	    K=K+1
	    W=A(K)
	    IF(IROW.EQ.ICOL) RSQ=(W*ETA)**2
	    M=J
	    DO 10 I=1,IROW
	      L=L+1
	      IF(I.EQ.IROW) GO TO 20
	      W=W-U(L)*U(M)
	      IF(IROW.EQ.ICOL) RSQ=RSQ+(U(L)**2*R(I))**2
	      M=M+1
   10       CONTINUE
   20       IF(IROW.EQ.ICOL) GO TO 50
	    IF(U(L).EQ.ZERO) GO TO 30
	    U(K)=W/U(L)
	    GO TO 40
   30       U(K)=ZERO
	    IF(ABS(W).GT.ABS(RMAX*A(K))) GO TO 100
   40     CONTINUE
C
C       END OF ROW, ESTIMATE RELATIVE ACCURACY OF DIAGONAL ELEMENT.
C
   50     RSQ=SQRT(RSQ)
	  IF(ABS(W).LE.FIVE*RSQ) GO TO 60
	  IF(W.LT.ZERO) GO TO 100
	  U(K)=SQRT(W)
	  R(I)=RSQ/W
	  IF(R(I).GT.RMAX) RMAX=R(I)
	  GO TO 70
   60     U(K)=ZERO
	  NULLTY=NULLTY+1
   70     J=J+ICOL
   80   CONTINUE
	IFAULT=0
C
  100   RETURN
	END




	SUBROUTINE NPROB(Z,P,Q,PDF)
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C       P, Q = PROBABILITIES TO THE LEFT AND RIGHT OF Z
C       FOR THE STANDARD NORMAL DISTRIBUTION.
C       PDF  = THE PROBABILITY DENSITY FUNCTION
C
C       REFERENCE: ADAMS,A.G. AREAS UNDER THE NORMAL CURVE,
C       ALGORITHM 39, COMPUTER J., VOL. 12, 197-8, 1969.
C
C       LATEST REVISION - 23 JANUARY 1981
C
C********************************************************************
C
	DATA A0,A1,A2,A3,A4,A5,A6,A7/0.5D0, 0.398942280444D0,
     1  0.399903438504D0, 5.75885480458D0, 29.8213557808D0,
     2  2.62433121679D0, 48.6959930692D0, 5.92885724438D0/,
     3  B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11/0.398942280385D0,
     4  3.8052D-8, 1.00000615302D0, 3.98064794D-4, 1.98615381364D0,
     5  0.151679116635D0, 5.29330324926D0, 4.8385912808D0,
     6  15.1508972451D0, 0.742380924027D0, 30.789933034D0,
     7  3.99019417011D0/
C
	ZABS = ABS(Z)
	IF(ZABS.GT.12.7D0) GO TO 20
	Y = A0*Z*Z
	PDF = EXP(-Y)*B0
	IF(ZABS.GT.1.28D0) GO TO 10
C
C       Z BETWEEN -1.28 AND +1.28
C
	Q = A0-ZABS*(A1-A2*Y/(Y+A3-A4/(Y+A5+A6/(Y+A7))))
	IF(Z.LT.0.D0) GO TO 30
	P = 1.D0-Q
	RETURN
C
C       ZABS BETWEEN 1.28 AND 12.7
C
   10   Q = PDF/(ZABS-B1+B2/(ZABS+B3+B4/(ZABS-B5+B6/(ZABS+B7-B8/
     1  (ZABS+B9+B10/(ZABS+B11))))))
	IF(Z.LT.0.D0) GO TO 30
	P = 1.D0-Q
	RETURN
C
C       Z FAR OUT IN TAIL
C
   20   Q = 0.D0
	PDF = 0.D0
	IF(Z.LT.0.D0) GO TO 30
	P = 1.D0
	RETURN
C
C       NEGATIVE Z, INTERCHANGE P AND Q
C
   30   P = Q
	Q = 1.D0-P
	RETURN
	END
