c - C中的Cholesky分解失败

我对c还是很陌生，想根据Wikipedia的伪代码实现cholesky分解。
要求ram是动态分配的。

我使用以下示例矩阵尝试了代码：

 4.000  2.000  0.000  0.000
 2.000  5.000  2.000  0.000
 0.000  2.000 10.000  3.000
 0.000  0.000  3.000  2.000

这应该导致：

2.000 0.000 0.000 0.000
1.000 2.000 0.000 0.000
0.000 1.000 3.000 0.000
0.000 0.000 1.000 1.000

而是让我回来。

4.000 2.000 0.000  0.000
0.000 5.000 2.000  0.000
0.000 0.000 10.000 3.000
0.000 0.000 0.000  2.000

我想我在使用指针时有些误解。我尝试根据this link动态分配。

谁能告诉我，为什么正确的值未写入矩阵？
这是我的代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int Cholesky(int n, double **A){
    double sum;
    sum = 0.0f;
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < i; j++)
        {
            for(int k = 0; k < j-1; k++)
            {
                sum = sum - A[i][k]*A[j][k];
            }
            if(i > j)
            {
                A[i][j] = sum / A[j][j];
            } else
            {
                if(sum > 0)
                {
                    A[i][i] = sqrt(sum);
                } else {
                    printf("Die Matrix ist nicht symetrisch positiv\n");
                    return -1;
                }
            }
        }
    }

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++)
            printf("%.5f ", A[i][j]);
            printf("\n");
    }

}

int main(){
    int n = 4;
    double ** matrix;

    double test[4][4] = {{4.0f,2.0f,0.0f,0.0f},{2.0f,5.0f,2.0f,0.0f},{0.0f,2.0f,10.0f, 3.0f},{0.0f,0.0f,3.0f,2.0f}};

    /* Speicher reservieren für die int-Zeiger (=zeile) */
    matrix = malloc(n * sizeof(double *));
    if(NULL == matrix) {
       printf("Kein virtueller RAM mehr vorhanden ... !");
       return -1;
    }
    /* jetzt noch Speicher reservieren für die einzelnen Spalten
     * der i-ten Zeile */
    for(int i = 0; i < n; i++) {
        matrix[i] = malloc(n * sizeof(double));
        if(NULL == matrix[i]) {
            printf("Kein Speicher mehr fuer Zeile %d\n",i);
            return -1;
        }
    }
    /* mit beliebigen Werten initialisieren */
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            matrix[i][j] = test[i][j];

    /* Inhalt der Matrix entsprechend ausgeben */
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++)
            printf("%.5f ", matrix[i][j]);
            printf("\n");
    }

    Cholesky(n, matrix);

    /* Spalten der i-ten Zeile freigeben */
    for(int i = 0; i < n; i++)
        free(matrix[i]);
    /* Jetzt können die leeren Zeilen freigegeben werden. */
    free(matrix);
    int x;
    scanf("%d", x);
    return 0;
}

最佳答案

您动态分配的矩阵是正确的，并且还使用了指针。

正确的值未写入矩阵，因为您用于cholesky分解的算法是错误的。您可以在此处找到正确的算法，例如：

http://www2.denizyuret.com/bib/press/www.library.cornell.edu/nr/bookcpdf/c2-9.pdf

您的函数仅按以下说明修改矩阵：

if (i > j) A[i][j] = sum / A[j][j];

这意味着对角线下方的元素将被清零，如您的帖子中所示。

4.000 2.000 0.000  0.000
0.000 5.000 2.000  0.000
0.000 0.000 10.000 3.000
0.000 0.000 0.000  2.000