- 62

- 2

I have to calculate the Cholesky decomposition of a symmetric matrix and this is the C ++ code I wrote:

C++:

```
boost::numeric::ublas::matrix<double> Math::cholesky(const boost::numeric::ublas::matrix<double> &MatrixA)
{
int dim = MatrixA.size1();
boost::numeric::ublas::matrix<double> MatrixL (dim,dim);
double sum;
for (int i = 0; i < dim; ++i)
{
for (int j = 0; j <= i; ++j)
{
sum = 0;
if (j == i)
{
for (int k = 0; k < j - 1; ++k)
{
sum += (MatrixL(j,k) * MatrixL(j,k));
MatrixL(j,j) = sqrt(MatrixA(j,j) - sum);
}
}
else
{
for (int k = 0; k < j; ++k)
{
sum += (MatrixL(i,k) * MatrixL(j,k));
MatrixL(i,j) = (MatrixA(i,j) - sum) / MatrixL(j,j);
}
}
}
}
return MatrixL;
}
```

My output is wrong:

MatrixA:

```
25 15 -5
15 18 0
-5 0 11
```

MatrixL:

```
4.6396e-310 0 0
0 3.2387e-319 2.122e-314
0 0 3.31662
```

thank you!