lu.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_SOLVER_LU_H
#define pRC_ALGORITHMS_SOLVER_LU_H
 
#include <prc/algorithms/lu.hpp>
#include <prc/algorithms/solve.hpp>
#include <prc/algorithms/solver/backward_substitution.hpp>
#include <prc/algorithms/solver/forward_substitution.hpp>
#include <prc/core/functors/eval.hpp>
#include <prc/core/tensor/declarations.hpp>
#include <prc/core/tensor/operator/hint.hpp>
#include <prc/core/tensor/operator/restrict.hpp>
#include <prc/core/tensor/operator/transform.hpp>
 
namespace pRC::Solver
{
    class LU
    {
    public:
        template<Operator::Transform OT = Operator::Transform::None,
            Operator::Restrict OR = Operator::Restrict::None,
            Operator::Hint OH = Operator::Hint::None, class XA,
            class RA = RemoveConstReference<XA>, class XB,
            class RB = RemoveReference<XB>, class XX = decltype(zero<RB>()),
            class RX = RemoveReference<XX>>
            requires IsInvocable<Apply<OT, OR, OH>, XA, XX>
        inline constexpr auto operator()(XA &&A, XB &&b) const
        {
            constexpr auto N = RA::size(0);
 
            auto A_ = lu<OH>(transform<OT>(forward<XA>(A)));
 
            Tensor Upper = upperTriangular(A_);
            Tensor Lower = strictlyLowerTriangular(A_);
            extractDiagonal(Lower) = unit();
 
            auto y = pRC::solve<pRC::Solver::ForwardSubstitution>(Lower, b);
            auto x = pRC::solve<pRC::Solver::BackwardSubstitution>(Upper, y);
 
            auto const z = eval(apply<OT, OR, OH>(A, x));
 
            if(!isApprox(z, b))
            {
                Logging::error("LU solver failed");
            }
 
            return x;
        }
    };
}
#endif // pRC_ALGORITHMS_SOLVER_LU_H