编一个程序，用C++实现矩阵的cholesky分解

给，下凯空迟面是Cholesky分解法的C++经典算法：

//-------------------------------------------------------------------

// Cholesky分解法

//-------------------------------------------------------------------

template <class T>

int cholesky(Matrix<T>& mat, double epsilon=EPSILON) {

size_t i, j, k

for (i=0 i<mat.Rows() ++i) {

// 计算第 i 轮主元

for (k=0 k<i ++k) mat[i][i] -= mat[k][i]*mat[k][i]

mat[i][i] = sqrt(mat[i][i])

// 计算第 i 轮主元结束

if (fabs(mat[i][i])<epsilon) break

//计算第 i 列

// for (j=i+1 j<mat.Rows() ++j) {

// for (k=0 k<i ++k) mat[j][i] -= mat[j][k]*mat[i][k]

// }

// 计算第 i 列结束

/盯李/ 计算第 i 行

for (j=i+1 j<mat.Cols() ++j) {

for (k=0 k<i ++k) mat[i][j] -= mat[k][i]*mat[k][j]

mat[i][j] /= mat[i][i]

}

// 计算第 i 行结束

}

return (i==mat.Rows())

}

下面这是另一种实现方法，输入输出语句自己根据需要写一下吧：

#include<malloc.h>

#include<math.h>

void cholesky(double **a,double *b,int n,double *x)

{

int i,j,m,k

double **L

L=(double **)malloc(n*sizeof(double))

for(i=0i<ni++)

L[i]=(double *)malloc(n*sizeof(double))

for(i=0i<ni++)

for(j=0j<nj++)

L[i][j]=0

for(k=0k<nk++)

{

L[k][k]=a[k][k]

for(m=0m<亏判km++)

L[k][k]-=L[k][m]*L[k][m]

L[k][k]=sqrt(L[k][k])

for(i=k+1i<ni++)

{

L[i][k]=a[i][k]

for(m=0m<km++)

L[i][k]-=L[i][m]*L[k][m]

L[i][k]/=L[k][k]

}

}

for(i=0i<ni++)

{

x[i]=b[i]

for(m=0m<im++)

x[i]-=L[i][m]*x[m]

x[i]/=L[i][i]

}

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

{

for(m=i+1m<nm++)

x[i]-=L[m][i]*x[m]

x[i]/=L[i][i]

}

}

void main()

{

int i,j,n=3

double **a,*b,*x

a=(double **)malloc(n*sizeof(double))

for(i=0i<ni++)

a[i]=(double *)malloc(n*sizeof(double))

for(i=0i<ni++)

for(j=0j<nj++)

a[i][j]=1

for(i=0i<ni++)

a[i][i]=100

b=(double *)malloc(n*sizeof(double))

x=(double *)malloc(n*sizeof(double))

for(i=0i<ni++)

b[i]=3

cholesky(a,b,n,x)

}

以前写的一段代码，看一祥埋下能用不。

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

double** alloc(int m, int n)

{

int i

double** a = (double**)malloc(m * sizeof(double*))

for (i = 0i <mi++)

a[i] = (double*)malloc(n * sizeof(double))

return a

}

void dealloc(double** a, int m)

{

int i

for (i = 0i <mi++)

free(a[i])

free(a)

}

void cholesky(double** a, int m, double* x)

{

int i, j, k

double s

a[0][0] = sqrt(a[0][0])

for (i = 0i <mi++)

a[i][0] = a[i][0] / a[0][0]

for (i = 1i <mi++) {

for (j = 0j <= ij++) {

for (s = 0, k = 0k <jk++)

s += a[i][k] * a[j][k]

if (i >j)

a[i][j] = (a[i][j] - s) / a[j][j]

else

a[i][j] = sqrt(a[i][j] - s)

}

}

a[0][m] = a[0][m] / a[0][0]

for (i = 1i <mi++) {

for (s = 0, k = 0k <谨耐蚂 ik++)

s += a[i][k] * a[k][m]

a[i][m] = (a[i][m] - s) / a[i][i]

}

x[m - 1] = a[m - 1][m] / a[m - 1][m - 1]

for (i = m - 2i >= 0i--) {

for (s = 0, k = i + 1k <mk++)

s += a[k][i] * x[k]

x[i] = (a[i][m] - s) / a[i][i]

}

}

int main()

{

int i, j, n

double** a, *x

FILE* fp = freopen("cholesky.txt", "r", stdin)

scanf("%d", &n)

a = alloc(n, n + 1)

for (i = 0i <ni++) {

for (j = 0j <n + 1j++) {

scanf("%lf", &a[i][j])

}

}

x = (double*)malloc(n * sizeof(double))

cholesky(a, n, x)

for (i = 0i <ni++)

printf("x%d = %f\亩稿n", i, x[i])

dealloc(a, n)

free(x)

return 0

}

---