c 中的 Cholesky 分解失败

Question

我对 c 相当陌生，想根据维基百科的伪代码实现 cholesky 分解。需要动态分配内存。

我用以下示例矩阵尝试了我的代码:

 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;
}

Answer 1

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

正确的值没有写入您的矩阵，因为您用于胆汁分解的算法是错误的。您可以在这里找到正确的算法:

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