Cholesky分解法的思想及C语言编程

设线性方程组Axb （1） 式中A为对称、正定的矩阵。对于对称、正定的矩阵A，可进行分解

ALDLT （2）

式中L是下单位三角阵，D是对角线矩阵。 右端项列向量（列阵）也作相应的分解

%bLDb （3）

将式（2）和式（3）代入方程（1），得到上三角方程组

%LTxb

再按诸如高斯消元法的回代过程就可解出x

1213x112x2354162下面，通过编写一个完整的C语言程序求，完整的C语言程序代码如下： 1551x3183417x420#include<> #include<> #include<> #define N 4

void Cholesky(int n,double A[N][N],double x[N],double b[N]) {

int i,j,k;

double L[N][N],D[N][N],b2[N]; i=1;

D[i-1][i-1]=A[i-1][i-1]; for(i=2;i<=n;i++) {

j=1;

L[i-1][j-1]=A[i-1][j-1]/D[j-1][j-1]; if(j==i-1) {

D[i-1][i-1]=0; for(k=1;k<=i-1;k++)

D[i-1][i-1]+=pow(L[i-1][k-1],2)*D[k-1][k-1];

D[i-1][i-1]=A[i-1][i-1]-D[i-1][i-1]; continue; } else {

for(j=2;j<=i-1;j++)

{

L[i-1][j-1]=0; for(k=1;k<=j-1;k++)

}

L[i-1][j-1]+=L[i-1][k-1]*L[j-1][k-1]*D[k-1][k-1]; L[i-1][j-1]=(A[i-1][j-1]-L[i-1][j-1])/D[j-1][j-1]; if(j==i-1) {

D[i-1][i-1]=0;

for(k=1;k<=i-1;k++) }

D[i-1][i-1]+=pow(L[i-1][k-1],2)*D[k-1][k-1]; D[i-1][i-1]=A[i-1][i-1]-D[i-1][i-1]; continue; } }

i=1;

b2[i-1]=b[i-1]/D[i-1][i-1];

for(i=2;i<=n;i++)

{

b2[i-1]=0;

for(k=1;k<=i-1;k++)

b2[i-1]+=L[i-1][k-1]*D[k-1][k-1]*b2[k-1]; b2[i-1]=(b[i-1]-b2[i-1])/D[i-1][i-1];

}

x[n-1]=b2[n-1] ;

for(i=n-1;i>=1;i--)

{

x[i-1]=;

for(k=i+1;k<=n;k++)

x[i-1]+=L[k-1][i-1]*x[k-1];

x[i-1]=b2[i-1]-x[i-1]; }

for(i=0;i<=n-1;i++)

}

printf(\"%f\\n\

void main() {

double A[4][4]={{1,2,1,3},{2,3,-5,4},{1,-5,5,-1},{3,4,-1,7}}; double x[4]={0};

double b[4]={12,16,18,20}; Cholesky(4,A,x,b); }

程序运行结果：