啊哈磊_编程从这里起步
标题:
求助!关于cholesky矩阵分解的代码,但是无法输出结果,望帮忙解决,十分感谢!!!
[打印本页]
作者:
我爱GM
时间:
2020-10-6 02:27
标题:
求助!关于cholesky矩阵分解的代码,但是无法输出结果,望帮忙解决,十分感谢!!!
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
double**MonAlloc(int n,int m)
{
int i,j;
double *vect,**mat;
vect=(double*)malloc((n*m)*sizeof(double));
mat=(double**)malloc(n*sizeof(double));
for(i=0;i<n;i++)
{
mat[i]=vect;
vect=vect+m;
}
return(mat);
}
double* sol(int N,double** A,double* b)
{
int i,j;
double sum,*x;
x=(double*)malloc(N*sizeof(double));
for(i=0;i<N;i++)
{
sum=0.0;
for(j=0;j<i;j++)
{
sum=sum + A[i][j]*x[j];
}
x[i]=(b[i]-sum)/A[i][i];
}
return(x);
}
double* Cholesky(int N,double** A,double* b)
{
int i,j,k;
double **L,**temp,*x,*y,sum;
L=MonAlloc(N,N);
temp=MonAlloc(N,N);
L[0][0]=sqrt(A[0][0]);
for(k=0;k<N;k++)
{
sum=0;
for(i=0;i<k;i++)
sum+=sum+L[k][j]*L[k][i];
sum=A[k][k] - sum;
L[k][k] = sqrt(sum>0? sum : 0);
for(i=k+1;i<N;i++)
{
sum=0;
for(j=0;j<k;j++)
sum += L[i][j]*L[k][j];
L[i][k]=(A[i][k]-sum)/L[k][k];
}
for(int j=0; j<k; j++)
L[j][k]=0;
}
y=sol(N,L,b);
for(i = 0;i < N;i++)
for(j = 0;j < i;j++)
{
temp[i][j] = L[i][j];
L[i][j] = L[j][i];
L[j][i] = temp[i][j];
}
x=sol(N,L,y);
return(x);
}
void main()
{
int i,N;
double **A,*b,*x;
N=4;
A=MonAlloc(N,N);
x=(double*)malloc(N*sizeof(double));
b=(double*)malloc(N*sizeof(double));
A[0][0]=10; A[0][1]=7; A[0][2]=8; A[0][3]=7;
A[1][0]=7; A[1][1]=5; A[1][2]=6; A[1][3]=5;
A[2][0]=8; A[2][1]=6; A[2][2]=10; A[2][3]=9;
A[3][0]=7; A[3][1]=5; A[3][2]=9; A[3][3]=10;
b[0]=32; b[0]=23; b[0]=33; b[0]=31;
Cholesky(N,A,b);
}
复制代码
欢迎光临 啊哈磊_编程从这里起步 (https://bbs.codeaha.com/)
Powered by Discuz! X3.2