backward_substitution.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_ALGORITHMS_SOLVER_BACKWARD_SUBSTITUTION_H
#define pRC_ALGORITHMS_SOLVER_BACKWARD_SUBSTITUTION_H
 
#include <prc/config.hpp>
#include <prc/core/basic/functions/copy.hpp>
#include <prc/core/basic/range.hpp>
#include <prc/core/tensor/functions/reshape.hpp>
#include <prc/core/tensor/operator/functions/is_upper_triangular.hpp>
#include <prc/core/tensor/operator/functions/restrict.hpp>
#include <prc/core/tensor/operator/functions/transform.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/limits.hpp>
 
namespace pRC::Solver
{
    class BackwardSubstitution
    {
    public:
        template<Operator::Transform OT = Operator::Transform::None,
            Operator::Restrict OR = Operator::Restrict::None,
            Operator::Hint OH = Operator::Hint::None, class XA,
            class RA = RemoveReference<XA>, class TA = typename RA::Type,
            class VA = typename TA::Value, class XB,
            class RB = RemoveReference<XB>, class TB = typename RB::Type,
            class VB = typename TB::Value,
            If<All<IsFloat<VA>, IsFloat<VB>>> = 0, If<IsTensorish<RA>> = 0,
            If<IsTensorish<RB>> = 0,
            If<IsSatisfied<(typename RA::Dimension() == 2)>> = 0,
            If<IsSatisfied<(typename RB::Dimension() == 1 ||
                typename RB::Dimension() == 2)>> = 0,
            If<IsSatisfied<(RA::size(0) == RA::size(1))>> = 0,
            If<IsSatisfied<(RA::size(0) == RB::size(0))>> = 0>
        inline constexpr auto operator()(XA &&A, XB &&b) const
        {
            auto const U = transform<OT>(restrict<OR>(forward<XA>(A)));
 
            if constexpr(cDebugLevel >= DebugLevel::High)
            {
                if(!isUpperTriangular(U))
                {
                    Logging::error("U is not upper triangular");
                }
            }
 
            decltype(auto) x = copy<!(!IsReference<XB>() && !IsConst<RB>() &&
                cDebugLevel < DebugLevel::High)>(eval(b));
 
            range<typename RB::Sizes>(
                [&U, &x](auto const p, auto const... indices)
                {
                    Index const k = RB::size(0) - 1 - p;
                    if(x(k, indices...) != zero())
                    {
                        x(k, indices...) /= U(k, k);
                        for(Index i = 0; i < k; ++i)
                        {
                            x(i, indices...) -= x(k, indices...) * U(i, k);
                        }
                    }
                });
 
            if constexpr(cDebugLevel >= DebugLevel::High)
            {
                if(!isApprox(U * x, b))
                {
                    Logging::error("Backward substitution failed.");
                }
            }
 
            if constexpr(IsReference<decltype(x)>())
            {
                return forward<XB>(b);
            }
            else
            {
                return x;
            }
        }
    };
}
#endif // pRC_ALGORITHMS_SOLVER_BACKWARD_SUBSTITUTION_H