C++，Lapack-Cholesky分解实现结果不准确

C++, Lapack Cholesky decomposition implementation inaccurate result

本文关键字：结果 不准确 实现 分解 Lapack-Cholesky C++      更新时间：2023-10-16

我正在尝试在C++中实现Cholesky分解，它以前是在lapack dpotrf_中实现的。

Cholesky分解：R' * R = A

代码：

#include <iostream>
#include <armadillo>


long my_chol(
    arma::mat &R,
    const arma::mat A,
    long lda
    )
{
    arma::arma_debug_check( (A.is_square() == false), "chol(): given matrix must be square sized" );
    arma::arma_debug_check( (lda<std::max(1l,(long)A.n_rows)), "chol(): LDA must be equal to or greater than max(1,N)" );
    double sum;
    long i, j, k;
    long n = A.n_rows;
    if(lda==(long)A.n_rows)
        R=A;
    else
    {
        R.zeros(lda,A.n_rows);
        R.submat(0,lda-1,0,A.n_rows-1)=A;
    }
    for(i=0; i<n; ++i)
        for(j=i+1; j<n; ++j)
            R(j,i)=0;

    for( i=0; i<n; ++i )
    {
        /* j == i */
        sum = R(i,i);
        for( k=(i-1); k>=0; --k )
            sum -= R(k,i)*R(k,i);
        if ( sum > 0.0 )
            R(i,i) = sqrt( sum );
        else
        {
            R(0) = sum; /* tunnel negative diagonal element to caller */
            return (long)i+1;
        }
        for( j=(i+1); j<n; ++j )
        {
            sum = R(i,j);
            for( k=(i-1); k>=0; --k )
                sum -= R(k,i) * R(k,i);
            R(i,j) = sum / R(i,i);
        }
    }
    return 0;
}

我使用以下代码测试此功能：

int main()
{
    arma::mat A={{10, 3, 5},{3, 60, 7},{5, 7, 9}};
    arma::mat B;
    my_chol(B,A,3);
    std::cout<<"---------------------------n";
    std::cout<<"A:n";
    A.print();
    std::cout<<"B:n";
    B.print();
    std::cout<<"---------------------------n";
    return 0;
}

这就是结果：

---------------------------
A:
   10.0000    3.0000    5.0000
    3.0000   60.0000    7.0000
    5.0000    7.0000    9.0000
B:
   3.1623   0.9487   1.5811
        0   7.6877   0.7935
        0        0   2.4229
---------------------------

但在倍频程中测试相同的矩阵会给我一个不同的结果：

A=[10,3,5;3,60,7;5,7,9];
chol(A)
   3.16228   0.94868   1.58114
   0.00000   7.68765   0.71543
   0.00000   0.00000   2.44707

尽管结果看起来很接近，但它们之间有细微的差异。

CCD_ 2和CCD_。我检查了结果。八度音阶的结果是正确的，而我的结果不是：

R=[3.1623,0.9487,1.5811;0,7.6877,0.7935;0,0,2.4229];
R'*R
ans =
   10.0001    3.0001    4.9999
    3.0001   60.0008    7.6002
    4.9999    7.6002    9.0000

为什么我的代码给出了错误的结果？

在最里面的循环中，它应该是

sum -= R(k,i) * R(k,j);

代替

sum -= R(k,i) * R(k,i);
相关文章：