cholesky.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_CHOLESKY_H
#define pRC_ALGORITHMS_CHOLESKY_H
 
#include <prc/config.hpp>
#include <prc/core/basic/functions/copy.hpp>
#include <prc/core/basic/range.hpp>
#include <prc/core/tensor/functions/inner_product.hpp>
#include <prc/core/tensor/operator/functions/is_self_adjoint.hpp>
#include <prc/core/tensor/operator/functions/is_upper_triangular.hpp>
#include <prc/core/tensor/operator/hint.hpp>
#include <prc/core/tensor/type_traits.hpp>
#include <prc/core/value/functions/is_approx.hpp>
#include <prc/core/value/functions/norm.hpp>
 
namespace pRC
{
    template<Operator::Hint H = Operator::Hint::None, class X,
        class R = RemoveReference<X>, If<IsTensorish<R>> = 0,
        If<IsFloat<typename R::Value>> = 0,
        If<IsSatisfied<(typename R::Dimension() == 2)>> = 0,
        If<IsSatisfied<(R::size(0) == R::size(1))>> = 0>
    static inline constexpr auto cholesky(X &&a)
    {
        constexpr auto N = R::size(0);
 
        decltype(auto) u = [&a]() -> decltype(auto)
        {
            if constexpr(H == Operator::Hint::UpperTriangular)
            {
                if constexpr(cDebugLevel >= DebugLevel::High)
                {
                    if(!isUpperTriangular(a))
                    {
                        Logging::error(
                            "Cholesky decomposition failed: Input is not upper "
                            "triangular part of self-adjoint matrix.");
                    }
 
                    return copy(eval(a));
                }
                else
                {
                    return copy<!(!IsReference<X>() && !IsConst<R>())>(eval(a));
                }
            }
            else if constexpr(H == Operator::Hint::LowerTriangular)
            {
                if constexpr(cDebugLevel >= DebugLevel::High)
                {
                    if(!isLowerTriangular(a))
                    {
                        Logging::error(
                            "Cholesky decomposition failed: Input is not lower "
                            "triangular part of self-adjoint matrix.");
                    }
                }
 
                return eval(transpose(a));
            }
            else if constexpr(H == Operator::Hint::None)
            {
                if constexpr(cDebugLevel >= DebugLevel::High)
                {
                    if(!isSelfAdjoint(a))
                    {
                        Logging::error(
                            "Cholesky decomposition failed: Input is not "
                            "self-adjoint.");
                    }
                }
 
                return eval(upperTriangular(a));
            }
            else
            {
                static_assert(H != H,
                    "Unsupported cholesky decomposition hint.");
            }
        }();
 
        range<N>(
            [&u](auto const i)
            {
                auto x = real(u(i, i));
                for(Index k = 0; k < i; ++k)
                {
                    x -= norm<2, 1>(u(k, i));
                }
                x = sqrt(x);
                u(i, i) = x;
                for(Index j = i + 1; j < N; ++j)
                {
                    for(Index k = 0; k < i; ++k)
                    {
                        u(i, j) -= innerProduct(u(k, i), u(k, j));
                    }
                    if(u(i, j) != zero())
                    {
                        u(i, j) /= x;
                    }
                }
            });
 
        if constexpr(cDebugLevel >= DebugLevel::High)
        {
            if(!isUpperTriangular(u))
            {
                Logging::error(
                    "Cholesky decomposition failed: U is not upper "
                    "triangular.");
            }
 
            if(!isSelfAdjoint(adjoint(u) * u))
            {
                Logging::error(
                    "Cholesky decomposition failed: Product is not self "
                    "adjoint.");
            }
 
            if(!isApprox(upperTriangular(a), upperTriangular(adjoint(u) * u)))
            {
                Logging::error("Cholesky decomposition failed.");
            }
        }
 
        if constexpr(IsReference<decltype(u)>())
        {
            return forward<X>(a);
        }
        else
        {
            return u;
        }
    }
}
#endif // pRC_ALGORITHMS_CHOLESKY_H