#include <stdio.h>
#include <math.h>
/* cholesky.c : 1999 3 14 by Oguni */
#define N 3
main(){
  double A[N][N]={3.,2.,1.,2.,3.,2.,1.,2.,3.},L[N][N];
  double tol=0.000001;
  int n=3, i, j;
  cholesky(A,L,tol,n);
  printf("Matrix L\n");
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      printf("  %f",L[i][j]); }
    printf("\n");}
}
/* function */
cholesky(A,L,tol,n)
double A[N][N],L[N][N],tol;
int n;
{
  double s[N],ss;
  int i,j,k;
  for (j=0; j<n; j++){
    if (j==0) {
      for (i=j; i<n; i++){ s[i]=A[i][j]; } 
    }
    if (j!=0) {
      for (i=j; i<n; i++){
        ss=0.;
        for (k=0; k<=j-1; k++){ ss+=L[i][k]*L[j][k];}
        s[i]=A[i][j]-ss; }
    }
    if (fabs(s[j])<tol) {
      printf("Too small pivot\n"); return; }
    for (i=j; i<n; i++){ L[i][j]=s[i]/sqrt(s[j]); }
  }
  return;
}