/*
 * to solve a positive definite sysmetric system ax=b
 *   a[1->n*(n+1)/2]        row-wise
 *   g[1->n*(n+1)/2]        cholesky storage
 *   b(1->n]                
 *   x[1->n]                answer vector
 *   n                      system dimension
 *   nono > 0 is the level at which p.d. failed
 *   (a,g) and (b,y,x) may be equivalenced.
 */

#include <math.h>

void choles(a,b,x,n,nono)
float *a,*b,*x;
int    n,*nono;
{
float *y,*g,*farray1();
int i,j,k,kj,kz,jmin,jmax,kmax;
float sg,gkz;
     
     y=farray1(1,n);
     g=a;
     *nono = 0;
     if(a[1]<0.0){
        *nono = 1;
        return;
     }
     
     g[1]=sqrt(a[1]);
     y[1]=b[1]/g[1];
     
     for(i=2;i<=n;i++){
         kz=(i*(i-1))/2;
         g[kz+1]=a[kz+1]/g[1];
         sg=g[kz+1]*g[kz+1];
         y[i]=b[i]-g[kz+1]*y[1];
         if(i>2){
            jmax=i-1;
            for(j=2;j<=jmax;j++){
                gkz=a[kz+j];
                kj=(j*(j-1))/2;
                kmax=j-1;
      
                for(k=1;k<=kmax;k++)
                    gkz-=g[kz+k]*g[kj+k];
          
                g[kz+j]=gkz/g[kj+j];
                y[i]-=g[kz+j]*y[j];
                sg+=g[kz+j]*g[kz+j];
            }
         }
      
         gkz=a[kz+i]-sg;
      
         if(gkz<=0.0){
            *nono=i;
            return;
         }
      
         g[kz+i]=sqrt(gkz);
         y[i]/=g[kz+i];
      }
      
      kz=(n*(n-1))/2;
      x[n]=y[n]/g[kz+n];
      if(n<=1) return;
      
      for(k=2;k<=n;k++){
          i=n+1-k;
          x[i]=y[i];
          jmin=i+1;
      
          for(j=jmin;j<=n;j++){
              kj=(j*(j-1))/2;
              x[i]-=g[kj+i]*x[j];
          }

          kz=(i*(i+1))/2;
          x[i]/=g[kz];
      }
}