source: CIVL/examples/compare/CholeskyDecomposition/cholesky.c

/* file: Choesky decomposition */

/* Take the cholesky decomposition in the manner described in FA Graybill
   (1976).
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>


int cholesky(double **orig, int n, double **aug, int mcol,double **chol, double **cholaug, int ofs)
     /* 
        Do the augmented cholesky decomposition as described in FA Graybill
        (1976) Theory and Application of the Linear Model. The original matrix
        must be symmetric positive definite. The augmentation matrix, or
        series of column vectors, are multiplied by C^-t, where C is the
        upper triangular cholesky matrix, ie C^t * C = M and M is the original
        matrix. Returns with a value of 0 if M is a non-positive definite 
        matrix. Returns with a value of 1 with succesful completion.

        Arguments:

        orig (input) double n x n array. The matrix to take the Cholesky
              decomposition of.
        n    (input) integer. Number of rows and columns in orig.
        aug  (input) double n x mcol array. The matrix for the augmented
              part of the decomposition.
        mcol (input) integer. Number of columns in aug.
        chol (output) double n x n array. Holds the upper triangular matrix
              C on output. The lower triangular portion remains unchanged.
              This maybe the same as orig, in which case the upper triangular
              portion of orig is overwritten.
        cholaug (output) double n x mcol array. Holds the product C^-t * aug.
                 May be the same as aug, in which case aug is over written.
        ofs (input) integer. The index of the first element in the matrices.
             Normally this is 0, but commonly is 1 (but may be any integer).
     */
{
   int i, j, k, l;
   int retval = 1;

   for (i=ofs; i<n+ofs; i++) {
      chol[i][i] = orig[i][i];
      for (k=ofs; k<i; k++)
         chol[i][i] -= chol[k][i]*chol[k][i];
      if (chol[i][i] <= 0) {
         fprintf(stderr,"\nERROR: non-positive definite matrix!\n");
         printf("\nproblem from %d %f\n",i,chol[i][i]);
         retval = 0;
         return retval;
      }
      chol[i][i] = sqrt(chol[i][i]);

      /*This portion multiplies the extra matrix by C^-t */
      for (l=ofs; l<mcol+ofs; l++) {
         cholaug[i][l] = aug[i][l];
         for (k=ofs; k<i; k++) {
            cholaug[i][l] -= cholaug[k][l]*chol[k][i];
         }
         cholaug[i][l] /= chol[i][i];
      }

      for (j=i+1; j<n+ofs; j++) {
         chol[i][j] = orig[i][j];
         for (k=ofs; k<i; k++)
            chol[i][j] -= chol[k][i]*chol[k][j];
         chol[i][j] /= chol[i][i];
      }
   }

   return retval;
}