cholesky.h

Go to the documentation of this file.

#ifndef FUSED_ALGORITHMS_CHOLESKY_H
#define FUSED_ALGORITHMS_CHOLESKY_H
 
#include "config.h"
#include "utilities/expected.h"
#include "matrix_traits.h"
#include "math/math_utils.h"    // math::sqrt
#include "simple_type_traits.h" // for move
 
namespace matrix_algorithms
{
 
    using matrix_traits::MatrixStatus;
 
    template <typename MatrixType>
    Expected<MatrixType, MatrixStatus> cholesky(
        const MatrixType &A,
        typename MatrixType::value_type tol = typename MatrixType::value_type(PRECISION_TOLERANCE))
    {
        static_assert(is_floating_point_v<typename MatrixType::value_type>,
                      "Cholesky requires a floating-point scalar type");
        if (!A.isSymmetric())
        {
            return Unexpected{MatrixStatus::NotSymmetric};
        }
 
        MatrixType L(0);
 
        for (my_size_t i = 0; i < A.getDim(0); ++i)
        {
            for (my_size_t j = 0; j <= i; ++j)
            {
                typename MatrixType::value_type sum = 0;
 
                for (my_size_t k = 0; k < j; ++k)
                {
                    sum += L(i, k) * L(j, k);
                }
 
                if (i == j)
                {
                    typename MatrixType::value_type diag = A(i, i) - sum;
 
                    if (diag <= tol)
                    {
                        return Unexpected{MatrixStatus::NotPositiveDefinite};
                    }
 
                    L(i, j) = math::sqrt(diag);
                }
                else
                {
                    L(i, j) = (A(i, j) - sum) / L(j, j);
                }
            }
        }
 
        return move(L);
    }
 
    template <typename MatrixType>
    MatrixType cholesky_or_die(const MatrixType &A)
    {
        auto result = cholesky(A);
 
        if (!result.has_value())
        {
            MyErrorHandler::error("cholesky decomposition failed");
        }
 
        return move(result.value());
    }
 
} // namespace matrix_algorithms
 
#endif // FUSED_ALGORITHMS_CHOLESKY_H