cholesky.c

用乔里斯基(cholesky)算法求解对称正定方程组。

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
void cholesky(double a[],double b[],int n)
{
	int i,j,k,m;
	double *d,*y,*xl,eps;
	d=(double *)malloc(n*sizeof(double));
	y=(double *)malloc(n*sizeof(double));
	xl=(double *)malloc(n*n*sizeof(double));
	eps=1.0e-15;
	m=0;
	d[0]=a[m];
	for(i=1;i<n;i++)
	{
		for(j=0;j<i;j++)
		{
			m=m+1;
		    xl[i*n+j]=a[m]/d[j];
		    if(j==0)continue;
		    for(k=0;k<j;k++)
			{xl[i*n+j]=xl[i*n+k]*xl[j*n+k]*d[k]/d[j];}
		}
		m=m+1;
	    d[i]=a[m];
	    for(k=0;k<i;k++)
		{d[i]=d[i]-d[k]*xl[i*n+k]*xl[i*n+k];}
		if(fabs(d[i])<eps)
		{
			printf("\nill-conditioned! \n");
			return;
		}
	}
	y[0]=b[0];
	for(k=1;k<n;k++)
	{
		y[k]=b[k];
		for(j=0;j<k;j++)
		{y[k]=y[k]-xl[k*n+j]*y[j];}
	}
	b[n-1]=y[n-1]/d[n-1];
	for(k=(n-2);k>=0;k--)
	{
		b[k]=y[k]/d[k];
		for(j=(k+1);j<n;j++)
		{b[k]=b[k]-xl[j*n+k]*b[j];}
	}
	free(d);
	free(y);
	free(xl);
}