C中的Cholesky分解？答案

【问题标题】：Cholesky Factorization in C?C中的Cholesky分解？
【发布时间】：2014-02-16 16:59:51
【问题描述】：

我正在 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;
}

我在这里做错了吗？

【问题讨论】：

你试过调试器吗？使用调试器，您至少可以确定有问题的行，这可能会让您自己清楚问题所在。这也可以告诉您错误是什么：它是否只是跳过函数的执行并因此终止？它会因分段错误而被杀死吗？浮点异常？
一个明确的问题是您假设 C 中的数组从 1 开始索引。
退出是什么意思？在什么时候？叫什么？
当我以前在 C 中编写这样的代码时，我广泛使用 Valgrind 来检查与内存相关的错误和错误。
你做了很多“错误”，但主要的是，对于现代 C，根本没有理由用指针表来模拟 2D 矩阵。它只会给你带来麻烦并且性能更差。

标签： c matrix linear-algebra factorization linear-equation

【解决方案1】：

double** cholesky(double **A, int N)

在此函数中，您假设数组长度为N。这意味着数组的最后一个索引位于N-1 而不是N。把代码改成：

for ( j = 0; j < N; ++j)

其他类似。

【讨论】：