Функция факторизации Cholesky от Matlab до C

Question

Я попытался преобразовать цикл Matlab, выполняющий cholesky-факторизацию в C. Я получил правильные результаты только для первого столбца.

Код Matlab:

M = [5 1.2 0.3 -0.6; 1.2 6 -0.4 0.9; 0.3 -0.4 8 1.7; -0.6 0.9 1.7 10]; 
n = length(M); 
L = zeros(n, n); 
for i=1:n 
    L(i, i) = sqrt(M(i, i) - L(i, :)*L(i, :)'); 

    for j=(i + 1):n 
     L(j, i) = (M(j, i) - L(i, :)*L(j, :)')/L(i, i); 
    end 
end 

Мой код C:

int main() 
{ 

    double M[4*4]={5, 1.2, 0.3, -0.6, 1.2, 6, -0.4, 0.9, 0.3, -0.4, 8, 1.7, -0.6, 0.9, 1.7, 10}; 
    int n=4; 
    double L[4*4]={0}; 
    double sum; 
    for (int i=0;i<n;i++){ 
     sum=0; 
     for (int j=0;j<n;j++){ 
      sum+=L[i*n+j]*L[j*n+i]; 
      } 
    L[i*n+i]=sqrt(M[i*n+i]-sum); 

    for (int j=i+1;j<n;j++){ 
     sum=0; 
     for(int k=0;k<n;k++){ 
      sum+=L[i*n+k]*L[k*n+j]; 
      } 
     L[j*n+i]=(M[j*n+i]-sum)/L[i*n+i]; 

     } 

} 
return 0; 

} 

Результат в Matlab является:

L = 

    2.2361   0   0   0 
    0.5367 2.3900   0   0 
    0.1342 -0.1975 2.8183   0 
    -0.2683 0.4368 0.6466 3.0527 

Но в C результат:

L_c = 

    2.2361   0   0   0 
    0.5367 2.4495   0   0 
    0.1342 -0.1633 2.8284   0 
    -0.2683 0.3674 0.6010 3.1623 

Где ошибка?

Answer 1

Проверьте линии

sum+=L[i*n+j]*L[j*n+i]; 

и

sum+=L[i*n+k]*L[k*n+j]; 

Я думаю, что это должно быть

sum+=L[i*n+j]*L[i*n+j]; 

и

sum+=L[i*n+k]*L[j*n+k]; 

соответственно. Транспонирование в коде Matlab просто изменяет направление вектора, но не переключает аргументы или что-то в этом роде.

Я советую вам сделать функцию, которая вычисляет точечные продукты, чтобы вы могли повторно использовать код.

Answer 2

Смотрите код:

#include <iostream> 
#include <vector>  
#include <iostream> 
#include <math.h> 
using namespace std; 

double M[4][4]={5, 1.2, 0.3, -0.6, 1.2, 6, -0.4, 0.9, 0.3, -0.4, 8, 1.7, -0.6, 0.9, 1.7, 10}; 
vector <vector<double> >L; 
int n=4; 
double calculateFunc1(int i){ 
double temp =0; 
for (int t = 0; t < n; t ++) 
    temp += L[i][t]*L[i][t]; 
return sqrt(M[i][i] - temp); 
} 

double calculateFunc2(int i,int j){ 
double temp =0; 
for (int t = 0; t < n; t ++) 
    temp += L[i][t]*L[j][t]; 
return (M[j][i] - temp)/L[i][i]; 
} 

int main() 
{ 
vector <double>temp; 
for (int i =0; i < n; i ++){ 
    for(int j =0; j < n; j ++){ 
     temp.push_back(0.0); 
    } 
    L.push_back(temp); 
    temp.clear(); 
} 

for(int i =0; i < n; i ++){ 
    L[i][i] = calculateFunc1(i); 

    for(int j= i +1 ; j < n; j++){ 
     L[j][i] = calculateFunc2(i,j); 
    } 
} 


for (int i=0;i<n;i++){ 

    for (int j=0;j<n;j++){ 
     cout<< L[i][j]<<" "; 
     } 
    cout << "\n"; 
} 
return 0; 
}