SUBROUTINE MINIMIZE_2D(X0,FCN,X,F,Info)	
	IMPLICIT NONE

!	Arguments

	REAL*8,  INTENT(INOUT) :: X0(2),X(2),F
	INTEGER, INTENT(INOUT) :: Info
	EXTERNAL FCN

!	Local Variables

	INTEGER Nop, MaxF, IPrint, NLoop,IQuad
	REAL*8 Step(2),StopCryterion,Simplex,Var(2)

!	Unpack and call the minimizer

	X(1) = X0(1)
	X(2) = X0(2)
	Step(1) = 1d-2
	Step(2) = 1d-2
	Nop = 2
	MaxF = 400
	IPrint = -1
	StopCryterion = 1.d-10
	NLoop = 4
	IQuad = 1
	Simplex = 1.d-4

	CALL minim(X, Step, Nop, F, MaxF, IPrint, StopCryterion, NLoop,   &
              IQuad, Simplex, Var, FCN, Info)

END SUBROUTINE MINIMIZE_2D

SUBROUTINE MINIMIZE_1D(X0,FCN,X,F,Info)
	IMPLICIT NONE

!	Arguments

	REAL*8,  INTENT(INOUT) :: X0,X,F
	INTEGER, INTENT(INOUT) ::  Info
	EXTERNAL FCN

!	Local Variables

	INTEGER Nop, MaxF, IPrint, NLoop,IQuad
	REAL*8 Step,StopCryterion,Simplex,Var

!	Unpack and call the minimizer

	X = X0
	Step = 0.4D0
	Nop = 1
	MaxF = 400
	IPrint = -1
	StopCryterion = 1.d-10
	NLoop = 4
	IQuad = 0
	Simplex = 1.d-5

	CALL minim(X, Step, Nop, F, MaxF, IPrint, StopCryterion, NLoop,   &
              IQuad, Simplex, Var, FCN, Info)

END SUBROUTINE MINIMIZE_1D


SUBROUTINE minim(p, step, nop, func, maxfn, iprint, stopcr, nloop, iquad,  &
                 simp, var, functn, ifault)

!     A PROGRAM FOR FUNCTION MINIMIZATION USING THE SIMPLEX METHOD.

!     FOR DETAILS, SEE NELDER & MEAD, THE COMPUTER JOURNAL, JANUARY 1965

!     PROGRAMMED BY D.E.SHAW,
!     CSIRO, DIVISION OF MATHEMATICS & STATISTICS
!     P.O. BOX 218, LINDFIELD, N.S.W. 2070

!     WITH AMENDMENTS BY R.W.M.WEDDERBURN
!     ROTHAMSTED EXPERIMENTAL STATION
!     HARPENDEN, HERTFORDSHIRE, ENGLAND

!     Further amended by Alan Miller
!     CSIRO Division of Mathematics & Statistics
!     Private Bag 10, CLAYTON, VIC. 3169

!     Fortran 90 conversion by Alan Miller, June 1995
!     Latest revision - 14 September 1995

!     ARGUMENTS:-
!     P()     = INPUT, STARTING VALUES OF PARAMETERS
!               OUTPUT, FINAL VALUES OF PARAMETERS
!     STEP()  = INPUT, INITIAL STEP SIZES
!     NOP     = INPUT, NO. OF PARAMETERS, INCL. ANY TO BE HELD FIXED
!     FUNC    = OUTPUT, THE FUNCTION VALUE CORRESPONDING TO THE FINAL
!                 PARAMETER VALUES.
!     maxfn     = INPUT, THE MAXIMUM NO. OF FUNCTION EVALUATIONS ALLOWED.
!               Say, 20 times the number of parameters, NOP.
!     IPRINT  = INPUT, PRINT CONTROL PARAMETER
!                 < 0 NO PRINTING
!                 = 0 PRINTING OF PARAMETER VALUES AND THE FUNCTION
!                     VALUE AFTER INITIAL EVIDENCE OF CONVERGENCE.
!                 > 0 AS FOR IPRINT = 0 PLUS PROGRESS REPORTS AFTER
!                     EVERY IPRINT EVALUATIONS, PLUS PRINTING FOR THE
!                     INITIAL SIMPLEX.
!     STOPCR  = INPUT, STOPPING CRITERION.
!               The criterion is applied to the standard deviation of
!               the values of FUNC at the points of the simplex.
!     NLOOP   = INPUT, THE STOPPING RULE IS APPLIED AFTER EVERY NLOOP
!               FUNCTION EVALUATIONS.   Normally NLOOP should be slightly
!               greater than NOP, say NLOOP = 2*NOP.
!     IQUAD   = INPUT, = 1 IF FITTING OF A QUADRATIC SURFACE IS REQUIRED
!                      = 0 IF NOT
!               N.B. The fitting of a quadratic surface is strongly
!               recommended, provided that the fitted function is
!               continuous in the vicinity of the minimum.   It is often
!               a good indicator of whether a premature termination of
!               the search has occurred.
!     SIMP    = INPUT, CRITERION FOR EXPANDING THE SIMPLEX TO OVERCOME
!               ROUNDING ERRORS BEFORE FITTING THE QUADRATIC SURFACE.
!               The simplex is expanded so that the function values at
!               the points of the simplex exceed those at the supposed
!               minimum by at least an amount SIMP.
!     VAR()   = OUTPUT, CONTAINS THE DIAGONAL ELEMENTS OF THE INVERSE OF
!               THE INFORMATION MATRIX.
!     FUNCTN  = INPUT, NAME OF THE USER'S SUBROUTINE - ARGUMENTS (P,FUNC)
!               WHICH RETURNS THE FUNCTION VALUE FOR A GIVEN SET OF
!               PARAMETER VALUES IN ARRAY P.
!****     FUNCTN MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM.
!     IFAULT  = OUTPUT, = 0 FOR SUCCESSFUL TERMINATION
!                 = 1 IF MAXIMUM NO. OF FUNCTION EVALUATIONS EXCEEDED
!                 = 2 IF INFORMATION MATRIX IS NOT +VE SEMI-DEFINITE
!                 = 3 IF NOP < 1
!                 = 4 IF NLOOP < 1

!     N.B. P, STEP AND VAR (IF IQUAD = 1) MUST HAVE DIMENSION AT LEAST NOP
!          IN THE CALLING PROGRAM.

!*****************************************************************************

IMPLICIT DOUBLE PRECISION (a-h,o-z)

INTEGER, INTENT(IN)             :: nop, maxfn, iprint, nloop, iquad
INTEGER, INTENT(OUT)            :: ifault
DOUBLE PRECISION, INTENT(IN)    :: stopcr, simp
DOUBLE PRECISION, INTENT(INOUT) :: p(nop), step(nop)
DOUBLE PRECISION, INTENT(OUT)   :: var(nop), func
EXTERNAL functn

!     Local variables

DOUBLE PRECISION   :: g(nop+1,nop), h(nop+1), pbar(nop), pstar(nop),      &
                      pstst(nop), aval(nop), pmin(nop), temp(nop),        &
                      bmat(nop*(nop+1)/2), vc(nop*(nop+1)/2),             &
                      zero = 0.d0, half = 0.5d0, one = 1.d0, two = 2.d0

!     A = REFLECTION COEFFICIENT, B = CONTRACTION COEFFICIENT, AND
!     C = EXPANSION COEFFICIENT.

DOUBLE PRECISION, PARAMETER :: a = 1.d0, b = 0.5d0, c = 2.d0

!     SET LOUT = LOGICAL UNIT NO. FOR OUTPUT

INTEGER, PARAMETER :: lout = 6

!     IF PROGRESS REPORTS HAVE BEEN REQUESTED, PRINT HEADING

IF (iprint.GT.0) WRITE (lout,5000) iprint

!     CHECK INPUT ARGUMENTS

ifault = 0
IF (nop.LE.0) ifault = 3
IF (nloop.LE.0) ifault = 4
IF (ifault.NE.0) RETURN

!     SET NAP = NO. OF PARAMETERS TO BE VARIED, I.E. WITH STEP.NE.0

nap = COUNT(step.NE.zero)
neval = 0
loop = 0
iflag = 0

!     IF NAP = 0 EVALUATE FUNCTION AT THE STARTING POINT AND RETURN

IF (nap.LE.0) THEN
  CALL functn(p,func)
  RETURN
END IF

!     SET UP THE INITIAL SIMPLEX

20 g(1,:) = p
irow = 2
DO i = 1, nop
  IF (step(i).NE.zero) THEN
    g(irow,:) = p
    g(irow,i) = p(i) + step(i)
    irow = irow + 1
  END IF
END DO

np1 = nap + 1
DO i = 1, np1
  p = g(i,:)
  CALL functn(p,h(i))
  neval = neval + 1
  IF (iprint.GT.0) THEN
    WRITE (lout,5100) neval, h(i), p
  END IF
END DO

!     START OF MAIN CYCLE.

!     FIND MAX. & MIN. VALUES FOR CURRENT SIMPLEX (HMAX & HMIN).

Main_loop: DO
  loop = loop + 1
  imax = 1
  imin = 1
  hmax = h(1)
  hmin = h(1)
  DO i = 2, np1
    IF (h(i).GT.hmax) THEN
      imax = i
      hmax = h(i)
    ELSE
      IF (h(i).LT.hmin) THEN
        imin = i
        hmin = h(i)
      END IF
    END IF
  END DO

!     FIND THE CENTROID OF THE VERTICES OTHER THAN P(IMAX)

  pbar = zero
  DO i = 1, np1
    IF (i.NE.imax) THEN
      pbar = pbar + g(i,:)
    END IF
  END DO
  pbar = pbar / FLOAT(nap)

!     REFLECT MAXIMUM THROUGH PBAR TO PSTAR,
!     HSTAR = FUNCTION VALUE AT PSTAR.

  pstar = a * (pbar - g(imax,:)) + pbar
  CALL functn(pstar,hstar)
  neval = neval + 1
  IF (iprint.GT.0) THEN
    IF (MOD(neval,iprint).EQ.0) WRITE (lout,5100) neval, hstar, pstar
  END IF

!     IF HSTAR < HMIN, REFLECT PBAR THROUGH PSTAR,
!     HSTST = FUNCTION VALUE AT PSTST.

  IF (hstar.LT.hmin) THEN
    pstst = c * (pstar - pbar) + pbar
    CALL functn(pstst,hstst)
    neval = neval + 1
    IF (iprint.GT.0) THEN
      IF (MOD(neval,iprint).EQ.0) WRITE (lout,5100) neval, hstst, pstst
    END IF

!     IF HSTST < HMIN REPLACE CURRENT MAXIMUM POINT BY PSTST AND
!     HMAX BY HSTST, THEN TEST FOR CONVERGENCE.

    IF (hstst.GE.hmin) THEN   ! REPLACE MAXIMUM POINT BY PSTAR & H(IMAX) BY HSTAR.
      g(imax,:) = pstar
      h(imax) = hstar
    ELSE
      g(imax,:) = pstst
      h(imax) = hstst
    END IF
    GO TO 250
  END IF

!     HSTAR IS NOT < HMIN.
!     TEST WHETHER IT IS < FUNCTION VALUE AT SOME POINT OTHER THAN
!     P(IMAX).   IF IT IS REPLACE P(IMAX) BY PSTAR & HMAX BY HSTAR.

  DO i = 1, np1
    IF (i.NE.imax) THEN
      IF (hstar.LT.h(i)) THEN  ! REPLACE MAXIMUM POINT BY PSTAR & H(IMAX) BY HSTAR.
        g(imax,:) = pstar
        h(imax) = hstar
        GO TO 250
      END IF
    END IF
  END DO

!     HSTAR > ALL FUNCTION VALUES EXCEPT POSSIBLY HMAX.
!     IF HSTAR <= HMAX, REPLACE P(IMAX) BY PSTAR & HMAX BY HSTAR.

  IF (hstar.LE.hmax) THEN
    g(imax,:) = pstar
    hmax = hstar
    h(imax) = hstar
  END IF

!     CONTRACTED STEP TO THE POINT PSTST,
!     HSTST = FUNCTION VALUE AT PSTST.

  pstst = b * g(imax,:) + (one-b) * pbar
  CALL functn(pstst,hstst)
  neval = neval + 1
  IF (iprint.GT.0) THEN
    IF (MOD(neval,iprint).EQ.0) WRITE (lout,5100) neval, hstst, pstst
  END IF

!     IF HSTST < HMAX REPLACE P(IMAX) BY PSTST & HMAX BY HSTST.

  IF (hstst.LE.hmax) THEN
    g(imax,:) = pstst
    h(imax) = hstst
    GO TO 250
  END IF

!     HSTST > HMAX.
!     SHRINK THE SIMPLEX BY REPLACING EACH POINT, OTHER THAN THE CURRENT
!     MINIMUM, BY A POINT MID-WAY BETWEEN ITS CURRENT POSITION AND THE
!     MINIMUM.

  DO i = 1, np1
    IF (i.NE.imin) THEN
      DO j = 1, nop
        IF (step(j).NE.zero) g(i,j) = (g(i,j) + g(imin,j)) * half
        p(j) = g(i,j)
      END DO
      CALL functn(p,h(i))
      neval = neval + 1
      IF (iprint.GT.0) THEN
        IF (MOD(neval,iprint).EQ.0) WRITE (lout,5100) neval, h(i), p
      END IF
    END IF
  END DO

!     IF LOOP = NLOOP TEST FOR CONVERGENCE, OTHERWISE REPEAT MAIN CYCLE.

  250 IF (loop.LT.nloop) CYCLE Main_loop

!     CALCULATE MEAN & STANDARD DEVIATION OF FUNCTION VALUES FOR THE
!     CURRENT SIMPLEX.

  hmean = SUM( h(1:np1) ) / FLOAT(np1)
  hstd = SUM( (h(1:np1) - hmean) ** 2 )
  hstd = SQRT(hstd / FLOAT(np1))

!     IF THE RMS > STOPCR, SET IFLAG & LOOP TO ZERO AND GO TO THE
!     START OF THE MAIN CYCLE AGAIN.

  IF (hstd.GT.stopcr .AND. neval.LE.maxfn) THEN
    iflag = 0
    loop = 0
    CYCLE Main_loop
  END IF

!     FIND THE CENTROID OF THE CURRENT SIMPLEX AND THE FUNCTION VALUE THERE.

  DO i = 1, nop
    IF (step(i).NE.zero) THEN
      p(i) = SUM( g(1:np1,i) ) / FLOAT(np1)
    END IF
  END DO
  CALL functn(p,func)
  neval = neval + 1
  IF (iprint.GT.0) THEN
    IF (MOD(neval,iprint).EQ.0) WRITE (lout,5100) neval, func, p
  END IF

!     TEST WHETHER THE NO. OF FUNCTION VALUES ALLOWED, maxfn, HAS BEEN
!     OVERRUN; IF SO, EXIT WITH IFAULT = 1.

  IF (neval.GT.maxfn) THEN
    ifault = 1
    IF (iprint.LT.0) RETURN
    WRITE (lout,5200) maxfn
    WRITE (lout,5300) hstd
    WRITE (lout,5400) p
    WRITE (lout,5500) func
    RETURN
  END IF

!     CONVERGENCE CRITERION SATISFIED.
!     IF IFLAG = 0, SET IFLAG & SAVE HMEAN.
!     IF IFLAG = 1 & CHANGE IN HMEAN <= STOPCR THEN SEARCH IS COMPLETE.

  IF (iprint.GE.0) THEN
    WRITE (lout,5600)
    WRITE (lout,5400) p
    WRITE (lout,5500) func
  END IF

  IF (iflag.EQ.0 .OR. ABS(savemn-hmean).GE.stopcr) THEN
    iflag = 1
    savemn = hmean
    loop = 0
  ELSE
    EXIT Main_loop
  END IF

END DO Main_loop

IF (iprint.GE.0) THEN
  WRITE (lout,5700) neval
  WRITE (lout,5800) p
  WRITE (lout,5900) func
END IF
IF (iquad.LE.0) RETURN

!------------------------------------------------------------------

!     QUADRATIC SURFACE FITTING

IF (iprint.GE.0) WRITE (lout,6000)

!     EXPAND THE FINAL SIMPLEX, IF NECESSARY, TO OVERCOME ROUNDING
!     ERRORS.

hmin = func
nmore = 0
DO i = 1, np1
  DO
    test = ABS(h(i)-func)
    IF (test.LT.simp) THEN
      DO j = 1, nop
        IF (step(j).NE.zero) g(i,j) = (g(i,j)-p(j)) + g(i,j)
        pstst(j) = g(i,j)
      END DO
      CALL functn(pstst,h(i))
      nmore = nmore + 1
      neval = neval + 1
      IF (h(i).GE.hmin) CYCLE
      hmin = h(i)
      IF (iprint.GE.0) WRITE (lout,5100) neval, hmin, pstst
    ELSE
      EXIT
    END IF
  END DO
END DO

!     FUNCTION VALUES ARE CALCULATED AT AN ADDITIONAL NAP POINTS.

DO i = 1, nap
  i1 = i + 1
  pstar = (g(1,:) + g(i1,:)) * half
  CALL functn(pstar,aval(i))
  nmore = nmore + 1
  neval = neval + 1
END DO

!     THE MATRIX OF ESTIMATED SECOND DERIVATIVES IS CALCULATED AND ITS
!     LOWER TRIANGLE STORED IN BMAT.

a0 = h(1)
DO i = 1, nap
  i1 = i - 1
  i2 = i + 1
  DO j = 1, i1
    j1 = j + 1
    pstst = (g(i2,:) + g(j1,:)) * half
    CALL functn(pstst,hstst)
    nmore = nmore + 1
    neval = neval + 1
    l = i * (i-1) / 2 + j
    bmat(l) = two * (hstst + a0 - aval(i) - aval(j))
  END DO
END DO

l = 0
DO i = 1, nap
  i1 = i + 1
  l = l + i
  bmat(l) = two * (h(i1) + a0 - two*aval(i))
END DO

!     THE VECTOR OF ESTIMATED FIRST DERIVATIVES IS CALCULATED AND
!     STORED IN AVAL.

DO i = 1, nap
  i1 = i + 1
  aval(i) = two * aval(i) - (h(i1) + 3.d0*a0) * half
END DO

!     THE MATRIX Q OF NELDER & MEAD IS CALCULATED AND STORED IN G.

pmin = g(1,:)
DO i = 1, nap
  i1 = i + 1
  g(i1,:) = g(i1,:) - g(1,:)
END DO

DO i = 1, nap
  i1 = i + 1
  g(i,:) = g(i1,:)
END DO

!     INVERT BMAT

CALL syminv(bmat, nap, bmat, temp, nullty, ifault, rmax)
IF (ifault.EQ.0) THEN
  irank = nap - nullty
ELSE                                 ! BMAT not +ve definite
                                     ! Resume search for the minimum
  IF (iprint.GE.0) WRITE (lout,6100)
  ifault = 2
  IF (neval.GT.maxfn) RETURN
  IF (iprint.GE.0) WRITE (lout,6200)
  step = half * step
  GO TO 20
END IF

!     BMAT*A/2 IS CALCULATED AND STORED IN H.

DO i = 1, nap
  h(i) = zero
  DO j = 1, nap
    IF (j.LE.i) THEN
      l = i * (i-1) / 2 + j
    ELSE
      l = j * (j-1) / 2 + i
    END IF
    h(i) = h(i) + bmat(l) * aval(j)
  END DO
END DO

!     FIND THE POSITION, PMIN, & VALUE, YMIN, OF THE MINIMUM OF THE
!     QUADRATIC.

ymin = DOT_PRODUCT( h(1:nap), aval(1:nap) )
ymin = a0 - ymin
DO i = 1, nop
  pstst(i) = DOT_PRODUCT( h(1:nap), g(1:nap,i) )
END DO
pmin = pmin - pstst
IF (iprint.GE.0) THEN
  WRITE (lout,6300) ymin, pmin
  WRITE (lout,6400)
END IF

!     Q*BMAT*Q'/2 IS CALCULATED & ITS LOWER TRIANGLE STORED IN VC

DO i = 1, nop
  DO j = 1, nap
    h(j) = zero
    DO k = 1, nap
      IF (k.LE.j) THEN
        l = j * (j-1) / 2 + k
      ELSE
        l = k * (k-1) / 2 + j
      END IF
      h(j) = h(j) + bmat(l) * g(k,i) * half
    END DO
  END DO

  DO j = i, nop
    l = j * (j-1) / 2 + i
    vc(l) = DOT_PRODUCT( h(1:nap), g(1:nap,j) )
  END DO
END DO

!     THE DIAGONAL ELEMENTS OF VC ARE COPIED INTO VAR.

j = 0
DO i = 1, nop
  j = j + i
  var(i) = vc(j)
END DO
IF (iprint.LT.0) RETURN
WRITE (lout,6500) irank
CALL print_tri_matrix(vc, nop, lout)

WRITE (lout,6600)
CALL syminv(vc, nap, bmat, temp, nullty, ifault, rmax)

!     BMAT NOW CONTAINS THE INFORMATION MATRIX

WRITE (lout,6700)
CALL print_tri_matrix(bmat, nop, lout)

ii = 0
ij = 0
DO i = 1, nop
  ii = ii + i
  IF (vc(ii).GT.zero) THEN
    vc(ii) = one / SQRT(vc(ii))
  ELSE
    vc(ii) = zero
  END IF
  jj = 0
  DO j = 1, i - 1
    jj = jj + j
    ij = ij + 1
    vc(ij) = vc(ij) * vc(ii) * vc(jj)
  END DO
  ij = ij + 1
END DO

WRITE (lout,6800)
ii = 0
DO i = 1, nop
  ii = ii + i
  IF (vc(ii).NE.zero) vc(ii) = one
END DO
CALL print_tri_matrix(vc, nop, lout)

!     Exit, on successful termination.

650 WRITE (lout,6900) nmore
RETURN

5000 FORMAT (' Progress Report every',i4,' function evaluations'/, &
             ' EVAL.   FUNC.VALUE.',10X,'PARAMETER VALUES')
5100 FORMAT (/1X,i4,2X,g12.5,2X,5G11.4,3(/21X,5G11.4))
5200 FORMAT (' No. of function evaluations > ',i5)
5300 FORMAT (' RMS of function values of last simplex =',g14.6)
5400 FORMAT (' Centroid of last simplex =',4(/1X,6G13.5))
5500 FORMAT (' Function value at centroid =',g14.6)
5600 FORMAT (/' EVIDENCE OF CONVERGENCE')
5700 FORMAT (/' Minimum found after',i5,' function evaluations')
5800 FORMAT (' Minimum at',4(/1X,6G13.6))
5900 FORMAT (' Function value at minimum =',g14.6)
6000 FORMAT (/' Fitting quadratic surface about supposed minimum'/)
6100 FORMAT (/' MATRIX OF ESTIMATED SECOND DERIVATIVES NOT +VE DEFN.'/ &
             ' MINIMUM PROBABLY NOT FOUND'/)
6200 FORMAT (/10X,'Search restarting'/)
6300 FORMAT (' Minimum of quadratic surface =',g14.6,' at',4(/1X,6G13.5))
6400 FORMAT (' IF THIS DIFFERS BY MUCH FROM THE MINIMUM ESTIMATED',1X,       &
             'FROM THE MINIMIZATION,'/' THE MINIMUM MAY BE FALSE &/OR THE '  &
             'INFORMATION MATRIX MAY BE',1X,'INACCURATE'/)
6500 FORMAT (' Rank of information matrix =',i3/ &
             ' Inverse of information matrix:-')
6600 FORMAT (/' If the function minimized was -LOG(LIKELIHOOD),'/         &
             ' this is the covariance matrix of the parameters.'/         &
             ' If the function was a sum of squares of residuals,'/       &
             ' this matrix must be multiplied by twice the estimated ',   &
             'residual variance'/' to obtain the covariance matrix.'/)
6700 FORMAT (' INFORMATION MATRIX:-'/)
6800 FORMAT (/' CORRELATION MATRIX:-')
6900 FORMAT (/' A further',i4,' function evaluations have been used'/)

END SUBROUTINE minim




SUBROUTINE syminv(a, n, c, w, nullty, ifault, rmax)

!     ALGORITHM AS7, APPLIED STATISTICS, VOL.17, 1968.

!     ARGUMENTS:-
!     A()    = INPUT, THE SYMMETRIC MATRIX TO BE INVERTED, STORED IN
!                LOWER TRIANGULAR FORM
!     N      = INPUT, ORDER OF THE MATRIX
!     C()    = OUTPUT, THE INVERSE OF A (A GENERALIZED INVERSE IF C IS
!                SINGULAR), ALSO STORED IN LOWER TRIANGULAR.
!                C AND A MAY OCCUPY THE SAME LOCATIONS.
!     W()    = WORKSPACE, DIMENSION AT LEAST N.
!     NULLTY = OUTPUT, THE RANK DEFICIENCY OF A.
!     IFAULT = OUTPUT, ERROR INDICATOR
!                 = 1 IF N < 1
!                 = 2 IF A IS NOT +VE SEMI-DEFINITE
!                 = 0 OTHERWISE
!     RMAX   = OUTPUT, APPROXIMATE BOUND ON THE ACCURACY OF THE DIAGONAL
!                ELEMENTS OF C.  E.G. IF RMAX = 1.E-04 THEN THE DIAGONAL
!                ELEMENTS OF C WILL BE ACCURATE TO ABOUT 4 DEC. DIGITS.

!     LATEST REVISION - 1 April 1985

!***************************************************************************

IMPLICIT DOUBLE PRECISION (a-h,o-z)
DIMENSION a(*), c(*), w(n)

DOUBLE PRECISION :: zero = 0.d0, one = 1.d0

nrow = n
ifault = 1
IF (nrow.GT.0) THEN
  ifault = 0

!     CHOLESKY FACTORIZATION OF A, RESULT IN C

  CALL chola(a, nrow, c, nullty, ifault, rmax, w)
  IF (ifault.EQ.0) THEN

!     INVERT C & FORM THE PRODUCT (CINV)'*CINV, WHERE CINV IS THE INVERSE
!     OF C, ROW BY ROW STARTING WITH THE LAST ROW.
!     IROW = THE ROW NUMBER, NDIAG = LOCATION OF LAST ELEMENT IN THE ROW.

    nn = nrow * (nrow+1) / 2
    irow = nrow
    ndiag = nn
    10 IF (c(ndiag).NE.zero) THEN
      l = ndiag
      DO i = irow, nrow
        w(i) = c(l)
        l = l + i
      END DO
      icol = nrow
      jcol = nn
      mdiag = nn

      30 l = jcol
      x = zero
      IF (icol.EQ.irow) x = one / w(irow)
      k = nrow
      40 IF (k.NE.irow) THEN
        x = x - w(k) * c(l)
        k = k - 1
        l = l - 1
        IF (l.GT.mdiag) l = l - k + 1
        GO TO 40
      END IF

      c(l) = x / w(irow)
      IF (icol.EQ.irow) GO TO 60
      mdiag = mdiag - icol
      icol = icol - 1
      jcol = jcol - 1
      GO TO 30
    END IF ! (c(ndiag).NE.zero)

    l = ndiag
    DO j = irow, nrow
      c(l) = zero
      l = l + j
    END DO

    60 ndiag = ndiag - irow
    irow = irow - 1
    IF (irow.NE.0) GO TO 10
  END IF
END IF
RETURN

END SUBROUTINE syminv





SUBROUTINE chola(a, n, u, nullty, ifault, rmax, r)

!     ALGORITHM AS6, APPLIED STATISTICS, VOL.17, 1968, WITH
!     MODIFICATIONS BY A.J.MILLER

!     ARGUMENTS:-
!     A()    = INPUT, A +VE DEFINITE MATRIX STORED IN LOWER-TRIANGULAR
!                FORM.
!     N      = INPUT, THE ORDER OF A
!     U()    = OUTPUT, A LOWER TRIANGULAR MATRIX SUCH THAT U*U' = A.
!                A & U MAY OCCUPY THE SAME LOCATIONS.
!     NULLTY = OUTPUT, THE RANK DEFICIENCY OF A.
!     IFAULT = OUTPUT, ERROR INDICATOR
!                 = 1 IF N < 1
!                 = 2 IF A IS NOT +VE SEMI-DEFINITE
!                 = 0 OTHERWISE
!     RMAX   = OUTPUT, AN ESTIMATE OF THE RELATIVE ACCURACY OF THE
!                DIAGONAL ELEMENTS OF U.
!     R()    = OUTPUT, ARRAY CONTAINING BOUNDS ON THE RELATIVE ACCURACY
!                OF EACH DIAGONAL ELEMENT OF U.

!     LATEST REVISION - 1 April 1985

!***************************************************************************

IMPLICIT DOUBLE PRECISION (a-h,o-z)
DIMENSION a(*), u(*), r(n)

!     ETA SHOULD BE SET EQUAL TO THE SMALLEST +VE VALUE SUCH THAT
!     1.D0 + ETA IS CALCULATED AS BEING GREATER THAN 1.D0 IN THE ACCURACY
!     BEING USED.

DOUBLE PRECISION :: eta = 1.d-16, zero = 0.d0

ifault = 1
IF (n.GT.0) THEN
  ifault = 2
  nullty = 0
  rmax = eta
  r(1) = eta
  j = 1
  k = 0

!     FACTORIZE COLUMN BY COLUMN, ICOL = COLUMN NO.

  DO 50 icol = 1, n
    l = 0

!     IROW = ROW NUMBER WITHIN COLUMN ICOL

    DO 30 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) EXIT
        w = w - u(l) * u(m)
        IF (irow.EQ.icol) rsq = rsq + (u(l)**2*r(i)) ** 2
        m = m + 1
      10 CONTINUE

      IF (irow.EQ.icol) GO TO 40
      IF (u(l).NE.zero) THEN
        u(k) = w / u(l)
      ELSE
        u(k) = zero
        IF (ABS(w).GT.ABS(rmax*a(k))) GO TO 60
      END IF
    30 CONTINUE

!     END OF ROW, ESTIMATE RELATIVE ACCURACY OF DIAGONAL ELEMENT.

    40 rsq = SQRT(rsq)
    IF (ABS(w).GT.5.*rsq) THEN
      IF (w.LT.zero) RETURN
      u(k) = SQRT(w)
      r(i) = rsq / w
      IF (r(i).GT.rmax) rmax = r(i)
    ELSE
      u(k) = zero
      nullty = nullty + 1
    END IF

    j = j + icol
  50 CONTINUE
  ifault = zero
END IF
60 RETURN

END SUBROUTINE chola



SUBROUTINE print_tri_matrix(a, n, lout)

IMPLICIT NONE
INTEGER, INTENT(IN)          :: n, lout
DOUBLE PRECISION, INTENT(IN) :: a(*)

!     Local variables
INTEGER  :: i, ii, i1, i2, l

l = 1
DO l = 1, n, 6
  ii = l * (l-1) / 2
  DO i = l, n
    i1 = ii + l
    ii = ii + i
    i2 = MIN(ii,i1+5)
    WRITE (lout,'(1X,6G13.5)') a(i1:i2)
  END DO
END DO
RETURN

END SUBROUTINE print_tri_matrix