我正在用 C 实现 Cholesky 方法,但程序在到达此点时退出。
答案之后:现在它可以工作了,感谢 (devnull & piotruś) 的答案,但它没有给我正确的答案
/* Ax=b
*This algorithm does:
* A = U * U'
* with
* U := lower left triangle matrix
* U' := the transposed form of U.
*/
double** cholesky(double **A, int N) //this gives me the U (edited)
{
int i, j, k;
double sum, **p, **U;
U=(double**)malloc(N*sizeof(double*));
for(p=U; p<U+N; p++)
*p=(double*)malloc(N*sizeof(double));
for (j=0; j<N; j++) {
sum = A[j][j];
for (k=0; k<(j-1); k++) sum -= U[k][j]*U[k][j];
U[j][j] = sqrt(sum);
for (i=j; i<N; i++) {
sum = A[j][i];
for (k=0; k<(j-1); k++)
sum -= U[k][j]*U[k][i];
U[j][i] = sum/U[j][j];
}
}
return U;
}
我在这里做错了什么吗?
最佳答案
double** cholesky(double **A, int N)
在此函数中,您假设数组长度为N
。这意味着数组的最后一个索引位于 N-1
处,而不是 N
处。将代码改为:
for ( j = 0; j < N; ++j)
其余类似。
关于C语言中的乔列斯基因式分解?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21814399/