forward_substitution.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_SOLVER_FORWARD_SUBSTITUTION_H
#define pRC_ALGORITHMS_SOLVER_FORWARD_SUBSTITUTION_H
 
#include <prc/core/functors/apply.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 ForwardSubstitution
    {
    public:
        template<Operator::Transform OT = Operator::Transform::None,
            Operator::Restrict OR = Operator::Restrict::None,
            Operator::Hint OH = Operator::Hint::None, class XA,
            IsTensorish RA = RemoveReference<XA>, class XB,
            IsTensorish RB = RemoveReference<XB>>
            requires(RA::Dimension == 2 &&
                        (RB::Dimension == 1 || RB::Dimension == 2)) &&
            IsInvocable<Apply<OT, OR, OH>, XA, XB>
        inline constexpr auto operator()(XA &&A, XB &&b) const
        {
            auto const L = transform<OT>(restrict<OR>(forward<XA>(A)));
 
            if constexpr(cDebugLevel >= DebugLevel::High)
            {
                if(!isLowerTriangular(L))
                {
                    Logging::error("L is not lower triangular");
                }
            }
 
            decltype(auto) x = copy<!(!IsReference<XB> && !IsConst<RB> &&
                cDebugLevel < DebugLevel::High)>(eval(b));
 
            range<typename RB::Sizes>(
                [&L, &x](auto const k, auto const... indices)
                {
                    if(x(k, indices...) != zero())
                    {
                        x(k, indices...) /= L(k, k);
                        for(Index i = k + 1; i < RB::size(0); ++i)
                        {
                            x(i, indices...) -= x(k, indices...) * L(i, k);
                        }
                    }
                });
 
            if constexpr(cDebugLevel >= DebugLevel::High)
            {
                if(!isApprox(b, L * x))
                {
                    Logging::error("Forward substitution failed.");
                }
            }
 
            if constexpr(IsReference<decltype(x)>)
            {
                return forward<XB>(b);
            }
            else
            {
                return x;
            }
        }
    };
}
#endif // pRC_ALGORITHMS_SOLVER_FORWARD_SUBSTITUTION_H