//$$ cholesky.cxx                     cholesky decomposition

// Copyright (C) 1991: R B Davies and DSIR

#define WANT_MATH

#include "include.hxx"

#include "newmat.hxx"


/********* Cholesky decomposition of a positive definite matrix *************/

// Suppose S is symmetrix and positive definite. Then there exists a unique
// lower triangular matrix L such that L L.t() = S;

ReturnMatrix Cholesky(const SymmetricMatrix& S)
{
   int nr = S.Nrows();
   LowerTriangularMatrix T(nr);
   real* s = S.Store(); real* t = T.Store(); real* tk = t;

   for (int i=0; i<nr; i++)
   {
      real* tj = t; real* ti = tk;
      for (int j=0; j<=i; j++)
      {
	 tk = ti; real sum = 0.0; int k = j;
	 while (k--) { sum += *tj++ * *tk++; }
	 if (i==j) *tk++ = sqrt(*s++ - sum);
         else *tk++ = (*s++ - sum) / *tj++;
      }
   }
   T.Release();              // not wanted if routine avoids return init
   return T;
}

