c***********************************************************************
c    Finds maximum likelihood estimates of parameters of symmetric, 
c  mean-reverting, correlated diffusions from discrete time 
c  observations.  State variables assumed to follow processes dx_i(t):
c
c     dx_i = k ( u - x_i ) dt + s ( p^1/2 dz_0 + (1-p)^1/2 dz_i )
c
c  where k, u are mean reversion parameters, s is common total volatilty
c  and p is common pairwise instantaneous correlation coefficient.
c  Processes z_0 and z_i, i=1..n, are independent standard Weiner
c  processes.
c     State variables are assumed observed at IT discrete intervals of
c  length H.
c
c  Inputs:     N(IT)   vector no. of correlated diffusions each t
c              IT      no. of observation dates
c              H       length of time between observations
c              X(*)    state variable levels at start of each interval
c                      concatenated: t=1 vector, t=2 vector, etc.
c              XH(*)   state variable levels time H later of same
c              IFLAG   integer flag selecting options: set to 
c                      IMr + 2 * ICov + 4 * IGenZ + 8 * INsame  
c                                      where
c                      IMr    = 0 (no drift)  or 1 (mean-reverting)
c                      ICov   = 0 (skip COV)  or 1 (return COV)
c                      IGenZ  = 0 (skip ZHAT) or 1 (return ZHAT)
c                      INsame = 0 (N's vary)  or 1 (N(1) applies to all) 
c
c  Outputs:    K       estimated mean reversion coefficient
c              U       estimated mean reversion target
c              SIG     estimated instantaneous volatility
c              RHO     estimated common correlation coeffient
c              LMAX    maximized 2 x log-likelihood value
c              ZHAT()  estimated common factor series
c              COV     lower triangular part of est. cov. matrix
c
c  Parameters: MAXT    max number of observation dates (change to suit)
c
c              Otherwise no limitations on size of data arrays
c
c  Author:     R. Jones  4 June 1999
c              add estim common factor series as output 24/8/99
c              add estim cov matrix of parameters 24/9/99
c              change to BHH covariance matrix estimate 9/10/99
c              change to permit incomplete data series 13/10/99
c***********************************************************************

      subroutine 
     & MAXLIKE( N, IT, H, X, XH, IFLAG, K, U, SIG, RHO, LMAX, ZHAT, COV)
      
      implicit double precision (A-H, K-L, O-Z)
      logical    MREVERT, GENCOV, GENZHAT, NSAME
      external   LAMBDA
      parameter (MAXT = 500)                     ! max no. of time points

      parameter (EPS = 1d-2)
      dimension  X(*), XH(*), ZHAT(IT), COV(10), N(IT)
      dimension  xe(MAXT), xhe(MAXT), xx(MAXT), xxh(MAXT), 
     &           xhxh(MAXT), g(MAXT), NN(MAXT)
      common / LLCOM / XFLAG, DAT(18), xe, xhe, xx, xxh, xhxh, g, NN
      equivalence ( DAT( 1), XT      ),
     &            ( DAT( 2), XN      ), 
     &            ( DAT( 3), Sxe     ), 
     &            ( DAT( 4), Sxx     ), 
     &            ( DAT( 5), Sxhe    ), 
     &            ( DAT( 6), Sxhxh   ), 
     &            ( DAT( 7), Sxxh    ), 
     &            ( DAT( 8), Skxe    ),
     &            ( DAT( 9), Skxhe   ), 
     &            ( DAT(10), Skxeex  ),
     &            ( DAT(11), Skxeexh ),
     &            ( DAT(12), Skxheexh),
     &            ( DAT(13), Skn     )
      
      if ( IT .gt. MAXT ) then
         print*,'MAXLIKE: Observation dates exceeds MAXT'
         return
      endif

c     Determine which options apply
      if ( IFLAG .gt. 15 .or. IFLAG .lt. 0 ) then
        print*, 'MAXLIKE.FOR: Inadmissible Iflag = ',IFLAG
        return
      endif
      II = IFLAG
      if ( mod(II,2) .eq. 1 ) then
         MREVERT = .true.
         else
         MREVERT = .false.
      endif
      II = II / 2
      if ( mod(II,2) .eq. 1 ) then
         GENCOV  = .true.
         else
         GENCOV  = .false.
      endif
      II = II / 2
      if ( mod(II,2) .eq. 1 ) then
         GENZHAT = .true.
         else
         GENZHAT = .false.
      endif
      II = II / 2
      if ( mod(II,2) .eq. 1 ) then
         NSAME   = .true.
         else
         NSAME   = .false.
      endif

c     Copy some stuff to COMMON
      XT = dble(IT)

      if ( NSAME ) then
         do I = 1, IT
           NN(I) = N(1)
         enddo
        else
         do I = 1, IT
           NN(I) = N(I)
         enddo
      endif
      
c     Accumulate sample moments of the data
      do I = 2, 13
         DAT(I) = 0d0
      enddo

c     For each observation date
      IX = 0                           ! position in X() vector
      XN = 0                           ! total no. of observations
      do J = 1, IT
         IN      = NN(J)
         xx(J)   = 0d0
         xe(J)   = 0d0
         xxh(J)  = 0d0
         xhe(J)  = 0d0
         xhxh(J) = 0d0
         XN = XN + IN
         do I = 1, IN  
            y    = X (IX+I)
            yh   = XH(IX+I)
            xx(J)  = xx(J) + y * y
            xxh(J) = xxh(J) + y * yh
            xhxh(J)  = xhxh(J) + yh * yh
            xe(J)  = xe(J) + y
            xhe(J) = xhe(J) + yh
         enddo
         IX = IX + IN  
      enddo

c     Accumulate total moments not depending on RHO
      Sxx   = 0d0
      Sxxh  = 0d0
      Sxhxh = 0d0
      Sxe   = 0d0
      Sxhe  = 0d0
      do J = 1, IT
         Sxx = Sxx + xx(J)
         Sxxh = Sxxh + xxh(J)
         Sxhxh = Sxhxh + xhxh(J)
         Sxe = Sxe + xe(J)
         Sxhe = Sxhe + xhe(J)
      enddo

c     Determine case
      if ( MREVERT ) then                        ! mean-reverting
         XFLAG = 1d0
      else                                       ! zero drift
         A     = 1d0
         B     = 0d0
         S     = (Sxhxh - 2d0 * Sxxh + Sxx) / XN
         XFLAG = 0d0
      endif

c     Find rho maximizing concentrated likelihood function 
      ax =  .00d0                                ! initial guess range
      bx = +.03d0
      call MNBRAK(ax,bx,cx,fa,fb,fc,LAMBDA)      ! bracket maximum
      tol = 1d-8
      LMAX = - BRENT(ax,bx,cx,LAMBDA,tol,RHOD)   ! find maximum

c     Compute remaining parameters a,b,s of discrete-time model:
      if ( MREVERT ) then
         RHO = RHOD
         TOP = Skn * Sxxh - RHO * Skn * Skxeexh - (1d0-RHO)*Skxe*Skxhe
         BOT = Skn * Sxx  - RHO * Skn * Skxeex  - (1d0-RHO)*Skxe**2   
         A   = TOP / BOT
         B   = (Skxhe - A * Skxe) / Skn
         TOP = Sxhxh + A*A*Sxx + B*B*XN - 2d0*(A*Sxxh + B*Sxhe -A*B*Sxe)
         S   = TOP / XN 
      endif

c     Convert to continuous-time model parameters:
      if ( .not. MREVERT ) then
        K = 0d0
        U = 0d0
        SIG = sqrt( S / H )
      elseif ( MREVERT ) then
        if ( A .gt. 0d0 ) then
           K = - log(A) / H
        else
           K = 1d1 / H
           A = exp( - K * H )
           print*,' Maxlike: warning A < 0 ... kappa = 10 output'
        endif
        U = B / ( 1d0 - A ) 
        SIG = sqrt( 2d0 * S * log(A) / ( H * (A**2 - 1d0 ) ) )
      endif
      RHO = RHOD

c     Get true 2 * log-likelihood function value:
      LMAX =  LMAX - XN * log(4d0 * datan(1d0)) 

c     Generate estimated common factor movements:
      if ( GENZHAT ) then
        if (  S * RHO  .lt. 1d-10 ) then         ! indeterminate or RHO < 0
          do J = 1, IT
             ZHAT(J) = 0d0
          enddo
          Sze  = 0d0
          Sze2 = 0d0
        else                                     ! positive RHO   
          Sze   = 0d0
          Sze2  = 0d0
          do J = 1, IT
             IN = NN(J)
             scale = 1d0 / ( IN *  sqrt(S * RHO / H ) )
             ze = ( xhe(J) - A * xe(J) - IN * B ) * scale
             ZHAT(J) = ze         
             Sze = Sze + ze
             Sze2 = Sze2 + ze * ze
          enddo
        endif
      endif 

c     Generate parameter cov estimate using Berndt/Hall/Hausman:
      if ( GENCOV ) then
         do I = 1, 10
            COV(I) = 0d0
         enddo
         do J = 1, IT
c           Generate numerical first partials of log Likelihood
            D = EPS * SIG  + 1d-5   
            DLDS = (LOGL(K,U,SIG+D,RHO,H,J)
     &            - LOGL(K,U,SIG-D,RHO,H,J)) / (2d0 * D)
            D = EPS * RHO  + 1d-5   
            DLDR = (LOGL(K,U,SIG,RHO+D,H,J)
     &            - LOGL(K,U,SIG,RHO-D,H,J)) / (2d0 * D)
          
            if ( MREVERT ) then
               D = EPS * K    + 1d-5   
               DLDK = (LOGL(K+D,U,SIG,RHO,H,J)
     &               - LOGL(K-D,U,SIG,RHO,H,J)) / (2d0 * D)
               D = EPS * U    + 1d-5   
               DLDU = (LOGL(K,U+D,SIG,RHO,H,J)
     &               - LOGL(K,U-D,SIG,RHO,H,J)) / (2d0 * D)
            endif

c           Accumulate cross product of partials into COV order S,R,K,U

            COV(1) = COV(1)  + DLDS * DLDS
            COV(2) = COV(2)  + DLDS * DLDR
            COV(3) = COV(3)  + DLDR * DLDR
      
            if ( MREVERT ) then
               COV(4) = COV(4)  + DLDS * DLDK
               COV(5) = COV(5)  + DLDR * DLDK
               COV(6) = COV(6)  + DLDK * DLDK
               COV(7) = COV(7)  + DLDS * DLDU
               COV(8) = COV(8)  + DLDR * DLDU
               COV(9) = COV(9)  + DLDK * DLDU
               COV(10)= COV(10) + DLDU * DLDU
            endif
         enddo

c        Invert information matrix estimate to get COV matrix
         if ( .not. MREVERT ) then
            call RSINV(COV,2,IER)
         elseif ( MREVERT ) then
            call RSINV(COV,4,IER)
         endif

         if (IER .eq. -1) then                   ! return 0 if not + def
           do i = 1, 10
              COV(i) = 0d0
           enddo
         endif     
      endif                                      ! end of GENCOV block

      return

      end
      

c***********************************************************************
c  Function calculates SINGLE DAY log-likelihood function for correlated 
c  diffusions with mean reversion from data moments and parameters.
c
c  Arguments:  KAP, U, SIG, R   values of model parameters
c              H                time interval between observations
c              I                observation number
c
c  Common   :  DAT   a vector containing the data moments as follows:
c                 1  N   number of jointly correlated diffusions
c                 2  T   number of discrete time observations
c                 3  Sum_t x'e
c                 4  Sum_t x'x
c                 5  Sum_t xh'e
c                 6  Sum_t xh'xh
c                 7  Sum_t x'ee'x
c                 8  Sum_t xh'ee'xh
c                 9  Sum_t x'xh
c                10  Sum_t x'ee'xh
c  Author:     R. Jones  15/10/99
c***********************************************************************

      double precision function LOGL( KAP, U, SIG, R, H, I )
      
      implicit double precision (A-Z)
      integer    NN, MAXT, I
      parameter (MAXT = 500)
      dimension  xe(MAXT), xhe(MAXT), xx(MAXT), xxh(MAXT), 
     &           xhxh(MAXT), g(MAXT), NN(MAXT)
      common / LLCOM / XFLAG, DAT(18), xe, xhe, xx, xxh, xhxh, g, NN
      
c     Transform to discrete time model parameters
      A = exp( - H * KAP )
      B = ( 1d0 - A ) * U
      if ( KAP .eq. 0d0 ) then
         S = SIG ** 2 * H
        else
         S = SIG ** 2 * ( 1d0 - A*A ) / ( 2d0 * KAP )
      endif
      R1 = 1d0 - R

      N = dble( NN(I) )
      k = 1d0 / ( 1d0 + ( N - 1d0 ) * R )
      PI2 = 4d0 * datan(1d0)

      Sxx = xx(I)
      Sxe = xe(I)
      Sxhe = xhe(I)
      Sxxh = xxh(I)
      Sxhxh = xhxh(I)      

      L = - N * log(PI2) - N * log(S) - (N-1d0) * log(R1) + log(k)
      JNK = Sxhxh + A*A * Sxx - 2d0 * A * Sxxh + B*B * R1 * k * N
     &    - R * k * ( Sxhe - A * Sxe ) ** 2
     &    - 2d0 * B * R1 * k * ( Sxhe - A * Sxe )
     
      L = L - JNK / (R1 * S)
      LOGL = L / 2d0

      return
      end


c***********************************************************************
c  Subroutine calculating CONCENTRATED likelihood function for 
c  estimation of symmetric correlated diffusions with mean reversion.
c  I.e., 2 log Likelihood maximized with respect to
c  mean reversion parameters a, b, s, and is a function just of rho
c  and the moments of the data.
c
c  Arguments:  RHO   a value for the common correlation coefficient
c              DATX  a vector containing the data moments as follows:
c                 1  N   number of jointly correlated diffusions
c                 2  T   number of discrete time observations
c                 3  Sum_t x'e
c                 4  Sum_t x'x
c                 5  Sum_t xh'e
c                 6  Sum_t xh'xh
c                 7  Sum_t x'ee'x
c                 8  Sum_t xh'ee'xh
c                 9  Sum_t x'xh
c                10  Sum_t x'ee'xh
c
c  Notes   : x(t)  = column n-vector of state values at time t
c            xh(t) = column n-vector of state values at time t+h
c            e     = column n-vector of 1's
c  ***      -log L is returned so Numerical Recipes minimization
c            routines can be used
c
c  Author  : 14 October 1999  R. Jones                          
c***********************************************************************

      double precision function LAMBDA( RHO ) 
      
      implicit double precision (A-Z)
      integer IT, NN, MAXT
      parameter (MAXT = 500)
      dimension  xe(MAXT), xhe(MAXT), xx(MAXT), xxh(MAXT), 
     &           xhxh(MAXT), k(MAXT), NN(MAXT)
      common / LLCOM / XFLAG, DAT(18), xe, xhe, xx, xxh, xhxh, k, NN
      equivalence ( DAT( 1), T       ),
     &            ( DAT( 2), N       ), 
     &            ( DAT( 3), Sxe     ), 
     &            ( DAT( 4), Sxx     ), 
     &            ( DAT( 5), Sxhe    ), 
     &            ( DAT( 6), Sxhxh   ), 
     &            ( DAT( 7), Sxxh    ), 
     &            ( DAT( 8), Skxe    ),
     &            ( DAT( 9), Skxhe   ), 
     &            ( DAT(10), Skxeex  ),
     &            ( DAT(11), Skxeexh ),
     &            ( DAT(12), Skxheexh),
     &            ( DAT(13), Skn     )      

c    (Check next condition for applicability with varying NN(t))

c     Constrain RHO to be in range [-1,1] by returning bad value:
      if  ( RHO .lt. (1d-10 - 1d0/(NN(1)-1)) 
     & .or. RHO .gt. (1d0 - 1d-10) ) then
         LAMBDA = 1d20
         return
      endif

      Slogz = 0d0
      do IT = 1, nint(T)                         ! needed for both
         g = 1d0 + ( NN(IT) - 1 ) * RHO
         Slogz = Slogz + log(g)
         k(IT) = 1d0 / g
      enddo

      if ( XFLAG .eq. 0d0 ) then                 ! no drift case
         Skyeey = 0d0
         do IT = 1, nint(T)
            Skyeey = Skyeey + k(IT) * (xhe(IT) - xe(IT))**2
         enddo
         S = (Sxhxh - 2d0 * Sxxh + Sxx) / N
         L = (N-T) * log(1d0 - RHO) + N * log(S) + Slogz 
     &     + N / (1d0 - RHO) - RHO * Skyeey / ((1d0 - RHO) * S)
         LAMBDA = L
         return
      elseif ( XFLAG .eq. 1d0 ) then             ! mean reverting
         Skn      = 0d0
         Skxe     = 0d0
         Skxhe    = 0d0
         Skxeexh  = 0d0
         Skxeex   = 0d0
         Skxheexh = 0d0
         do IT = 1, nint(T)
           kt   = k(IT)
           xet  = xe(IT)
           xhet = xhe(IT)
           Skn      = Skn + kt * NN(IT)
           Skxe     = Skxe + kt * xet
           Skxhe    = Skxhe + kt * xhet
           Skxeexh  = Skxeexh + kt * xet * xhet
           Skxeex   = Skxeex + kt * xet ** 2
           Skxheexh = Skxheexh + kt * xhet ** 2
         enddo

         TOP = Skn * Sxxh - RHO * Skn * Skxeexh - (1d0-RHO)*Skxe*Skxhe
         BOT = Skn * Sxx  - RHO * Skn * Skxeex  - (1d0-RHO)*Skxe**2   
         A   = TOP / BOT
         B   = (Skxhe - A * Skxe) / Skn
         TOP = Sxhxh + A*A*Sxx + B*B*N - 2d0*(A*Sxxh + B*Sxhe -A*B*Sxe)
         S   = TOP / N 
         JNK = Sxhxh + A*A*Sxx - 2d0*A*Sxxh + B*B*(1d0-RHO)*Skn
     &       - RHO * (Skxheexh - 2d0 * A * Skxeexh + A*A * Skxeex)
     &       - 2d0 * B * (1d0-RHO) * (Skxhe - A * Skxe)
         L   = (N-T) * log(1d0 - RHO) + N * log(S) + Slogz 
     &       + JNK / (1d0-RHO) / S
         LAMBDA = L
         return
      else                                       ! flag set wrong
         LAMBDA = 0d0
         return
      endif
      end
            
c***********************************************************************
      SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func)
      DOUBLE PRECISION ax,bx,cx,fa,fb,fc,func,GOLD,GLIMIT,TINY
      EXTERNAL func
      PARAMETER (GOLD=1.618034d0, GLIMIT=100.d0, TINY=1.d-20)
      DOUBLE PRECISION dum,fu,q,r,u,ulim
      fa=func(ax)
      fb=func(bx)
      if(fb.gt.fa)then
        dum=ax
        ax=bx
        bx=dum
        dum=fb
        fb=fa
        fa=dum
      endif
      cx=bx+GOLD*(bx-ax)
      fc=func(cx)
1     if(fb.ge.fc)then
        r=(bx-ax)*(fb-fc)
        q=(bx-cx)*(fb-fa)
        u=bx-((bx-cx)*q-(bx-ax)*r)/(2.d0*sign(max(abs(q-r),TINY),q-r))
        ulim=bx+GLIMIT*(cx-bx)
        if((bx-u)*(u-cx).gt.0.d0)then
          fu=func(u)
          if(fu.lt.fc)then
            ax=bx
            fa=fb
            bx=u
            fb=fu
            return
          else if(fu.gt.fb)then
            cx=u
            fc=fu
            return
          endif
          u=cx+GOLD*(cx-bx)
          fu=func(u)
        else if((cx-u)*(u-ulim).gt.0.d0)then
          fu=func(u)
          if(fu.lt.fc)then
            bx=cx
            cx=u
            u=cx+GOLD*(cx-bx)
            fb=fc
            fc=fu
            fu=func(u)
          endif
        else if((u-ulim)*(ulim-cx).ge.0.d0)then
          u=ulim
          fu=func(u)
        else
          u=cx+GOLD*(cx-bx)
          fu=func(u)
        endif
        ax=bx
        bx=cx
        cx=u
        fa=fb
        fb=fc
        fc=fu
        goto 1
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software *1`3.d0
      
c***********************************************************************
      FUNCTION brent(ax,bx,cx,f,tol,xmin)
      INTEGER ITMAX
      DOUBLE PRECISION brent,ax,bx,cx,tol,xmin,f,CGOLD,ZEPS
      EXTERNAL f
      PARAMETER (ITMAX=100,CGOLD=.3819660d0,ZEPS=1.d-18)
      INTEGER iter
      DOUBLE PRECISION a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x
     *,xm
      a=min(ax,cx)
      b=max(ax,cx)
      v=bx
      w=v
      x=v
      e=0.d0
      fx=f(x)
      fv=fx
      fw=fx
      do 11 iter=1,ITMAX
        xm=0.5d0*(a+b)
        tol1=tol*abs(x)+ZEPS
        tol2=2.d0*tol1
        if(abs(x-xm).le.(tol2-.5d0*(b-a))) goto 3
        if(abs(e).gt.tol1) then
          r=(x-w)*(fx-fv)
          q=(x-v)*(fx-fw)
          p=(x-v)*q-(x-w)*r
          q=2.d0*(q-r)
          if(q.gt.0.d0) p=-p
          q=abs(q)
          etemp=e
          e=d
          if(abs(p).ge.abs(.5d0*q*etemp).or.p.le.q*(a-x).or.p.ge.q*(b-x)
     *) 
     *goto 1
          d=p/q
          u=x+d
          if(u-a.lt.tol2 .or. b-u.lt.tol2) d=sign(tol1,xm-x)
          goto 2
        endif
1       if(x.ge.xm) then
          e=a-x
        else
          e=b-x
        endif
        d=CGOLD*e
2       if(abs(d).ge.tol1) then
          u=x+d
        else
          u=x+sign(tol1,d)
        endif
        fu=f(u)
        if(fu.le.fx) then
          if(u.ge.x) then
            a=x
          else
            b=x
          endif
          v=w
          fv=fw
          w=x
          fw=fx
          x=u
          fx=fu
        else
          if(u.lt.x) then
            a=u
          else
            b=u
          endif
          if(fu.le.fw .or. w.eq.x) then
            v=w
            fv=fw
            w=u
            fw=fu
          else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then
            v=u
            fv=fu
          endif
        endif
11    continue
      pause 'brent exceed maximum iterations'
3     xmin=x
      brent=fx
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software *1`3.d0
      

C**** SUBROUTINE TO INVERT REAL SYMMETRIC MATRIX: IBM SUBROUTINE LIBRARY

      SUBROUTINE RSINV(AMAT,N,IER) 
      IMPLICIT REAL*8(A-H,O-Z)  
      DIMENSION  AMAT(1) 
      CALL RMFSD(AMAT,N,IER) 
      IF(IER) 9,1,1
    1 IPIV = N*(N+1)/2
      IND = IPIV  
      DO 6 I = 1,N
      DIN = 1.D0/AMAT(IPIV)  
      AMAT(IPIV) = DIN
      MIN = N  
      KEND = I - 1
      LANF = N - KEND 
      IF(KEND) 5,5,2  
    2 J = IND  
      DO 4 K = 1,KEND 
      WORK = 0.D0 
      MIN = MIN - 1
      LHOR = IPIV 
      LVER = J 
      DO 3  L = LANF,MIN 
      LVER = LVER + 1 
      LHOR = LHOR + L 
    3 WORK = WORK + AMAT(LVER)*AMAT(LHOR) 
      AMAT(J) = -WORK*DIN
    4 J = J - MIN 
    5 IPIV = IPIV - MIN  
    6 IND = IND - 1
      DO 8 I = 1,N
      IPIV = IPIV + I 
      J = IPIV 
      DO 8 K = I,N
      WORK = 0.D0 
      LHOR = J 
      DO 7 L = K,N
      LVER = LHOR + K - I
      WORK = WORK + AMAT(LHOR)*AMAT(LVER) 
    7 LHOR = LHOR + L 
      AMAT(J) = WORK  
    8 J = J + K
    9 RETURN
      END

      SUBROUTINE RMFSD(AMAT,N,IER)                                      
      IMPLICIT REAL*8(A-D,F-H,O-S,U-Z)                                  
      DIMENSION AMAT(1)                                                 
      EPS = 0.0                                                         
      IF(N-1) 12,1,1                                                    
    1 IER = 0                                                           
      KPIV = 0                                                          
      DO 11 K = 1,N                                                     
      KPIV = KPIV + K                                                   
      IND = KPIV                                                        
      LEND = K - 1                                                      
      TOL = ABS(EPS*SNGL(AMAT(KPIV)))                                   
      DO 11 I = K,N                                                     
      DSUM = 0.D0                                                       
      IF(LEND) 2,4,2                                                    
    2 DO 3 L = 1,LEND                                                   
      LANF = KPIV - L                                                   
      LIND = IND - L                                                    
    3 DSUM = DSUM +AMAT(LANF)*AMAT(LIND)                                
    4 DSUM = AMAT(IND) - DSUM                                           
      IF(I-K) 10,5,10                                                   
    5 IF(SNGL(DSUM)-TOL) 6,6,9                                          
    6 IF(DSUM) 12,12,7                                                  
    7 IF(IER) 8,8,9                                                     
    8 IER = K - 1                                                       
    9 DPIV = DSQRT(DSUM)                                                
      AMAT(KPIV) = DPIV                                                 
      DPIV = 1.D0/DPIV                                                  
      GO TO 11                                                          
   10 AMAT(IND) =DSUM*DPIV                                              
   11 IND  = IND+ I                                                     
      RETURN                                                            
   12 IER = -1                                                          
      RETURN                                                            
      END                                                               
 
