c语言实现cholesky分解

/* file: cholesky.c */


/* 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;
}