C++，Lapack乔列斯基分解实现不准确的结果

Question

我试图在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 

虽然结果看起来近，它们之间有细微的差别。

R23和R33略有不同。我检查了结果。从八度的结果是正确的，我的一个不是：

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 

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

Answer 1

在最里面的循环应该是

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

而不是

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