adaptors.hpp

Go to the documentation of this file.

#ifndef TDLS_TFEL_ADAPTORS_HPP
#define TDLS_TFEL_ADAPTORS_HPP
 
 
 
 
 
 
#include <cstddef>
#include <type_traits>
 
#include <tdls/core/macros.hpp>
#include <tdls/solvers/tiled_lupp/solver_dynamic.hpp>
#include <tdls/solvers/tiled_lupp/solver_static.hpp>
 
 
 
namespace tdls {
 
 
 
namespace detail {
 
template<typename T, typename = void>
struct has_data : std::false_type {};
template<typename T>
struct has_data<T, std::void_t<decltype(std::declval<const T&>().data())>>
    : std::is_pointer<decltype(std::declval<const T&>().data())> {};
 
template<typename T, typename = void>
struct has_pair_data : std::false_type {};
template<typename T>
struct has_pair_data<T, std::void_t<decltype(std::declval<const T&>().data().first),
                                    decltype(std::declval<const T&>().data().second)>>
    : std::bool_constant<std::is_pointer_v<decltype(std::declval<const T&>().data().first)> &&
                         std::is_integral_v<decltype(std::declval<const T&>().data().second)>> {};
 
template<typename T, bool = has_pair_data<T>::value>
struct const_data_pointer {
    using type = decltype(std::declval<const T&>().data());
};
template<typename T>
struct const_data_pointer<T, true> {
    using type = decltype(std::declval<const T&>().data().first);
};
 
template<typename T, bool = has_pair_data<T>::value>
struct mutable_data_pointer {
    using type = decltype(std::declval<T&>().data());
};
template<typename T>
struct mutable_data_pointer<T, true> {
    using type = decltype(std::declval<T&>().data().first);
};
 
template<typename T, typename = void>
struct has_stride : std::false_type {};
template<typename T>
struct has_stride<T, std::void_t<decltype(std::declval<const T&>().stride())>> : std::true_type {};
 
template<typename T, typename = void>
struct has_get_stride : std::false_type {};
template<typename T>
struct has_get_stride<T, std::void_t<decltype(std::declval<const T&>().getStride())>>
    : std::true_type {};
 
template<typename T, typename = void>
struct has_indexing_policy : std::false_type {};
template<typename T>
struct has_indexing_policy<T, std::void_t<typename T::indexing_policy>> : std::true_type {};
 
template<typename T>
inline constexpr bool is_dense_v =
    (has_data<T>::value || has_pair_data<T>::value) && has_indexing_policy<T>::value;
 
} // namespace detail
 
 
 
template<typename DenseType, typename = void>
struct storage_traits;
 
template<typename DenseType>
struct storage_traits<DenseType, std::enable_if_t<detail::is_dense_v<DenseType>>> {
    using indexing_policy = typename DenseType::indexing_policy;
    using policy_size_type = typename indexing_policy::size_type;
    using value_type = std::remove_cv_t<
        std::remove_pointer_t<typename detail::const_data_pointer<DenseType>::type>>;
    static constexpr bool is_mutable = !std::is_const_v<
        std::remove_pointer_t<typename detail::mutable_data_pointer<DenseType>::type>>;
    static constexpr int arity = static_cast<int>(indexing_policy::arity);
    static_assert(arity == 1 || arity == 2, "tdls adaptors: only vector-like (arity 1) and "
                                            "matrix-like (arity 2) objects are supported");
    static constexpr int extent0 = static_cast<int>(indexing_policy{}.size(0));
    static constexpr int extent1 = (arity == 2) ? static_cast<int>(indexing_policy{}.size(1)) : 1;
 
    static constexpr int compute_policy_stride() {
        if constexpr (arity == 1) {
            return static_cast<int>(indexing_policy{}.getIndex(policy_size_type(1)));
        } else {
            return static_cast<int>(
                indexing_policy{}.getIndex(policy_size_type(0), policy_size_type(1)));
        }
    }
    static constexpr int policy_stride = compute_policy_stride();
 
    [[nodiscard]] static constexpr bool has_uniform_rows() {
        if constexpr (arity == 1) {
            return true;
        } else {
            return static_cast<int>(indexing_policy{}.getIndex(
                       policy_size_type(1), policy_size_type(0))) == extent1 * policy_stride;
        }
    }
    // Single-stride addressing requires uniform rows. Row-strided matrix
    // views (sub-matrix views) violate this and cannot be passed to the
    // solvers; solve the full matrix or use an element-strided
    // (coalesced) view instead.
    static_assert(has_uniform_rows(),
                  "tdls adaptors: row-strided matrix views (sub-matrix views) cannot "
                  "be expressed by the single-stride addressing of the TiledLUpp solvers");
 
    static constexpr bool has_runtime_stride = detail::has_pair_data<DenseType>::value ||
                                               detail::has_stride<DenseType>::value ||
                                               detail::has_get_stride<DenseType>::value;
    static constexpr bool is_internal = (policy_stride == 1) && !has_runtime_stride;
 
    static constexpr bool has_runtime_extents = (extent0 == 0);
 
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr int
    runtime_extent0(const DenseType& o) noexcept {
        return static_cast<int>(o.getIndexingPolicy().size(0));
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr int
    runtime_extent1(const DenseType& o) noexcept {
        return static_cast<int>(o.getIndexingPolicy().size(1));
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr int
    runtime_row_stride(const DenseType& o) noexcept {
        return static_cast<int>(
            o.getIndexingPolicy().getIndex(policy_size_type(1), policy_size_type(0)));
    }
 
    [[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr int
    runtime_col_stride(const DenseType& o) noexcept {
        return static_cast<int>(
            o.getIndexingPolicy().getIndex(policy_size_type(0), policy_size_type(1)));
    }
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr value_type*
    pointer(const DenseType& o) noexcept {
        if constexpr (detail::has_pair_data<DenseType>::value) {
            return const_cast<value_type*>(static_cast<const value_type*>(o.data().first));
        } else {
            return const_cast<value_type*>(static_cast<const value_type*>(o.data()));
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr int
    runtime_view_stride(const DenseType& o) noexcept {
        if constexpr (detail::has_pair_data<DenseType>::value) {
            return static_cast<int>(o.data().second);
        } else if constexpr (detail::has_stride<DenseType>::value) {
            return static_cast<int>(o.stride());
        } else if constexpr (detail::has_get_stride<DenseType>::value) {
            return static_cast<int>(o.getStride());
        } else {
            return 1;
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr int stride(const DenseType& o) noexcept {
        return policy_stride * runtime_view_stride(o);
    }
 
    // The three accessors below describe a matrix-like block as the
    // _multirhs entry points of the raw API see it: a column count, a
    // total distance between row-consecutive elements and a total
    // distance between column-consecutive elements. Each reads the
    // policy type when the extents are compile-time and the policy
    // instance of the object otherwise, so a fixed-size block needs no
    // getIndexingPolicy() member, as the structural contract states.
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr int columns(const DenseType& o) noexcept {
        if constexpr (has_runtime_extents) {
            return runtime_extent1(o);
        } else {
            return extent1;
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr int row_stride(const DenseType& o) noexcept {
        if constexpr (has_runtime_extents) {
            return runtime_row_stride(o) * runtime_view_stride(o);
        } else {
            return extent1 * stride(o);
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr int col_stride(const DenseType& o) noexcept {
        if constexpr (has_runtime_extents) {
            return runtime_col_stride(o) * runtime_view_stride(o);
        } else {
            return stride(o);
        }
    }
};
 
 
 
namespace detail {
 
template<typename T, typename = void>
struct has_storage_traits : std::false_type {};
template<typename T>
struct has_storage_traits<T, std::void_t<typename storage_traits<T>::value_type>> : std::true_type {
};
 
template<typename MatrixType>
struct matrix_scalar_of {
    static_assert(has_storage_traits<std::remove_cv_t<MatrixType>>::value,
                  "tdls adaptors: A must be a dense object exposing data() and an indexing_policy "
                  "type (a gather view holding one pointer per element is not one)");
    using type = typename storage_traits<std::remove_cv_t<MatrixType>>::value_type;
};
 
template<typename MatrixType>
using matrix_scalar = typename matrix_scalar_of<MatrixType>::type;
 
template<typename ConfigType>
struct is_solver_config : std::false_type {};
template<typename Scalar>
struct is_solver_config<TiledLUppConfig<Scalar>> : std::true_type {};
 
template<typename ConfigType>
concept solver_config = is_solver_config<std::remove_cvref_t<ConfigType>>::value;
 
template<typename MatrixType, auto UserConfig>
consteval auto checked_config() {
    static_assert(std::is_same_v<std::remove_cvref_t<decltype(UserConfig)>,
                                 TiledLUppConfig<matrix_scalar<MatrixType>>>,
                  "tdls adaptors: the config scalar type must match the matrix element type");
    return UserConfig;
}
 
template<typename MatrixType, TiledLUppConfig<matrix_scalar<MatrixType>> Config>
struct adaptor_context {
    using mtraits = storage_traits<std::remove_cv_t<MatrixType>>;
    using scalar = typename mtraits::value_type;
    static_assert(std::is_floating_point_v<scalar>,
                  "tdls adaptors: the element type must be float, double or long double");
    static_assert(mtraits::arity == 2, "tdls adaptors: A must be matrix-like");
    static constexpr bool runtime_sized = mtraits::has_runtime_extents;
    static_assert(runtime_sized || mtraits::extent0 == mtraits::extent1,
                  "tdls adaptors: A must be square");
    static_assert(Config.layout == MatrixLayout::RowMajor,
                  "tdls adaptors: dense objects are addressed row-major, the TFEL convention; "
                  "a column-major configuration cannot be used through the adaptors");
    static constexpr int N = mtraits::extent0;
    using solver = TiledLUppSolverStatic<scalar, N, Config>;
    using dynamic_solver = TiledLUppSolverDynamic<scalar, Config>;
};
 
template<typename VectorType, typename Scalar, int N>
TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr void check_vector() {
    using vtraits = storage_traits<std::remove_cv_t<VectorType>>;
    static_assert(vtraits::arity == 1, "tdls adaptors: expected a vector-like object");
    static_assert(N == 0 || vtraits::has_runtime_extents || vtraits::extent0 == N,
                  "tdls adaptors: vector extent does not match the system dimension");
    static_assert(std::is_same_v<typename vtraits::value_type, Scalar>,
                  "tdls adaptors: mixed scalar types");
}
 
template<typename RhsType, typename SolutionType>
TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr void check_rhs_pair() {
    using bt = storage_traits<std::remove_cv_t<RhsType>>;
    using xt = storage_traits<std::remove_cv_t<SolutionType>>;
    static_assert(bt::is_internal == xt::is_internal && bt::policy_stride == xt::policy_stride &&
                      bt::has_runtime_stride == xt::has_runtime_stride,
                  "tdls adaptors: b and x must share the same residency and "
                  "stride (the raw API carries a single rhs_stride for both)");
    static_assert(bt::has_runtime_extents || xt::has_runtime_extents || bt::extent0 == xt::extent0,
                  "tdls adaptors: b and x extents do not match");
}
 
template<typename RhsType, typename SolutionType, typename Scalar, int N>
TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr void check_multirhs_pair() {
    using bt = storage_traits<std::remove_cv_t<RhsType>>;
    using xt = storage_traits<std::remove_cv_t<SolutionType>>;
    static_assert(xt::arity == 2, "tdls adaptors: a matrix-like right-hand side requires a "
                                  "matrix-like solution");
    static_assert(std::is_same_v<typename bt::value_type, Scalar> &&
                      std::is_same_v<typename xt::value_type, Scalar>,
                  "tdls adaptors: mixed scalar types");
    static_assert(bt::has_runtime_extents == xt::has_runtime_extents,
                  "tdls adaptors: B and X must both be fixed-size or both runtime-sized");
    static_assert(N == 0 || !bt::has_runtime_extents,
                  "tdls adaptors: a runtime-sized right-hand-side block requires a "
                  "runtime-sized matrix A");
    static_assert(bt::has_runtime_extents ||
                      (bt::extent0 == xt::extent0 && bt::extent1 == xt::extent1 &&
                       (N == 0 || bt::extent0 == N)),
                  "tdls adaptors: B and X must have N rows and the same column count");
    static_assert(bt::policy_stride == xt::policy_stride &&
                      bt::has_runtime_stride == xt::has_runtime_stride,
                  "tdls adaptors: B and X must share the same layout (the raw API "
                  "carries a single stride pair for both)");
}
 
template<typename PivotType, int N>
struct pivot_access {
    static constexpr bool is_raw =
        std::is_pointer_v<std::decay_t<PivotType>> &&
        std::is_same_v<std::remove_cv_t<std::remove_pointer_t<std::decay_t<PivotType>>>, int>;
    static constexpr bool is_dense = has_storage_traits<std::remove_cv_t<PivotType>>::value;
    static_assert(is_raw || is_dense, "tdls adaptors: the pivot must be an int pointer/array or "
                                      "a dense int object");
 
    static constexpr bool has_matching_shape = [] {
        if constexpr (is_raw) {
            return true;
        } else {
            using pt = storage_traits<std::remove_cv_t<PivotType>>;
            return pt::arity == 1 && (N == 0 || pt::has_runtime_extents || pt::extent0 == N);
        }
    }();
    static_assert(has_matching_shape, "tdls adaptors: the pivot must be vector-like, of the "
                                      "system dimension when both are fixed-size");
 
    static constexpr bool is_mutable = [] {
        if constexpr (is_raw) {
            return !std::is_const_v<std::remove_pointer_t<std::decay_t<PivotType>>>;
        } else {
            return !std::is_const_v<PivotType> &&
                   storage_traits<std::remove_cv_t<PivotType>>::is_mutable;
        }
    }();
 
    static constexpr bool is_internal = [] {
        if constexpr (is_raw) {
            return true;
        } else {
            return storage_traits<std::remove_cv_t<PivotType>>::is_internal;
        }
    }();
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr int* pointer(const PivotType& p) noexcept {
        if constexpr (is_raw) {
            return const_cast<int*>(static_cast<const int*>(p));
        } else {
            using pt = storage_traits<std::remove_cv_t<PivotType>>;
            static_assert(std::is_same_v<typename pt::value_type, int>,
                          "tdls adaptors: the pivot element type must be int");
            return pt::pointer(p);
        }
    }
 
    TDLS_HOST_DEVICE TDLS_FORCEINLINE static constexpr int stride(const PivotType& p) noexcept {
        if constexpr (is_raw) {
            return 1;
        } else {
            return storage_traits<std::remove_cv_t<PivotType>>::stride(p);
        }
    }
};
 
template<typename Ctx, int pass_width, bool inplace, typename MatrixType, typename PivotType,
         typename RhsType, typename SolutionType>
TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr void
substitute_multirhs_dispatch(const MatrixType& A, const PivotType& piv, const RhsType& b,
                             SolutionType& x) {
    using mt = typename Ctx::mtraits;
    using pa = pivot_access<PivotType, Ctx::N>;
    using bt = storage_traits<std::remove_cv_t<RhsType>>;
    using xt = storage_traits<std::remove_cv_t<SolutionType>>;
    static_assert(pass_width >= 0, "tdls adaptors: pass_width must not be negative");
    if constexpr (Ctx::runtime_sized) {
        const int n    = mt::runtime_extent0(A);
        const int m    = bt::columns(b);
        const int rs   = xt::row_stride(x);
        const int xcol = xt::col_stride(x);
        if constexpr (inplace) {
            Ctx::dynamic_solver::template substitute_inplace_multirhs<pass_width>(
                n, m, mt::pointer(A), mt::stride(A), pa::pointer(piv), pa::stride(piv),
                xt::pointer(x), rs, xcol);
        } else {
            Ctx::dynamic_solver::template substitute_multirhs<pass_width>(
                n, m, mt::pointer(A), mt::stride(A), pa::pointer(piv), pa::stride(piv),
                bt::pointer(b), xt::pointer(x), rs, xcol);
        }
    } else {
        constexpr int M = bt::extent1;
        static_assert(M >= 1, "tdls adaptors: B must have at least one column");
        const int rs   = xt::row_stride(x);
        const int xcol = xt::col_stride(x);
        if constexpr (inplace) {
            Ctx::solver::template substitute_inplace_multirhs<M, false, pa::is_internal,
                                                              mt::is_internal, pass_width>(
                mt::pointer(A), mt::stride(A), pa::pointer(piv), pa::stride(piv), xt::pointer(x),
                rs, xcol);
        } else {
            Ctx::solver::template substitute_multirhs<M, false, pa::is_internal, mt::is_internal,
                                                      pass_width>(
                mt::pointer(A), mt::stride(A), pa::pointer(piv), pa::stride(piv), bt::pointer(b),
                xt::pointer(x), rs, xcol);
        }
    }
}
 
template<auto UserConfig, bool oot_diagnostics, typename MatrixType, typename PivotType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
factorize_dispatch(MatrixType& A, PivotType& piv, int& oot_count) {
    using ctx = adaptor_context<MatrixType, checked_config<MatrixType, UserConfig>()>;
    using mt  = typename ctx::mtraits;
    using pa  = pivot_access<PivotType, ctx::N>;
    static_assert(mt::is_mutable, "tdls adaptors: factorize writes into A");
    static_assert(!std::is_const_v<MatrixType>, "tdls adaptors: A must not be const here "
                                                "(factorize writes it)");
    static_assert(pa::is_mutable, "tdls adaptors: the pivot must be mutable here "
                                  "(factorize writes it)");
    if constexpr (ctx::runtime_sized) {
        return ctx::dynamic_solver::template factorize<oot_diagnostics>(
            mt::runtime_extent0(A), mt::pointer(A), mt::stride(A), pa::pointer(piv),
            pa::stride(piv), oot_count);
    } else {
        return ctx::solver::template factorize<pa::is_internal, mt::is_internal, oot_diagnostics>(
            mt::pointer(A), mt::stride(A), pa::pointer(piv), pa::stride(piv), oot_count);
    }
}
 
template<auto UserConfig, int pass_width, bool oot_diagnostics, typename MatrixType,
         typename PivotType, typename RhsType, typename SolutionType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
solve_dispatch(MatrixType& A, PivotType& piv, const RhsType& b, SolutionType& x, int& oot_count) {
    using ctx = adaptor_context<MatrixType, checked_config<MatrixType, UserConfig>()>;
    using mt  = typename ctx::mtraits;
    using pa  = pivot_access<PivotType, ctx::N>;
    using bt  = storage_traits<std::remove_cv_t<RhsType>>;
    using xt  = storage_traits<std::remove_cv_t<SolutionType>>;
    static_assert(mt::is_mutable, "tdls adaptors: solve factors A in place");
    static_assert(xt::is_mutable, "tdls adaptors: solve writes into x");
    static_assert(!std::is_const_v<MatrixType> && !std::is_const_v<SolutionType>,
                  "tdls adaptors: A and x must not be const here (solve writes them)");
    static_assert(pa::is_mutable, "tdls adaptors: the pivot must be mutable here "
                                  "(solve writes it)");
    if constexpr (bt::arity == 2) {
        check_multirhs_pair<RhsType, SolutionType, typename ctx::scalar, ctx::N>();
        if (!factorize_dispatch<UserConfig, oot_diagnostics>(A, piv, oot_count)) return false;
        substitute_multirhs_dispatch<ctx, pass_width, false>(A, piv, b, x);
        return true;
    } else {
        static_assert(pass_width == 0,
                      "tdls adaptors: pass_width only applies to matrix-like right-hand sides");
        check_vector<RhsType, typename ctx::scalar, ctx::N>();
        check_vector<SolutionType, typename ctx::scalar, ctx::N>();
        check_rhs_pair<RhsType, SolutionType>();
        if constexpr (ctx::runtime_sized) {
            return ctx::dynamic_solver::template solve<oot_diagnostics>(
                mt::runtime_extent0(A), mt::pointer(A), mt::stride(A), pa::pointer(piv),
                pa::stride(piv), bt::pointer(b), xt::pointer(x), xt::stride(x), oot_count);
        } else {
            return ctx::solver::template solve<xt::is_internal, pa::is_internal, mt::is_internal,
                                               oot_diagnostics>(
                mt::pointer(A), mt::stride(A), pa::pointer(piv), pa::stride(piv), bt::pointer(b),
                xt::pointer(x), xt::stride(x), oot_count);
        }
    }
}
 
template<auto UserConfig, int pass_width, bool oot_diagnostics, typename MatrixType,
         typename PivotType, typename VectorType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
solve_inplace_dispatch(MatrixType& A, PivotType& piv, VectorType& y, int& oot_count) {
    using ctx = adaptor_context<MatrixType, checked_config<MatrixType, UserConfig>()>;
    using mt  = typename ctx::mtraits;
    using pa  = pivot_access<PivotType, ctx::N>;
    using yt  = storage_traits<std::remove_cv_t<VectorType>>;
    static_assert(mt::is_mutable, "tdls adaptors: solve_inplace factors A in place");
    static_assert(yt::is_mutable, "tdls adaptors: solve_inplace writes into y");
    static_assert(!std::is_const_v<MatrixType> && !std::is_const_v<VectorType>,
                  "tdls adaptors: A and y must not be const here (solve_inplace writes them)");
    static_assert(pa::is_mutable, "tdls adaptors: the pivot must be mutable here "
                                  "(solve_inplace writes it)");
    if constexpr (yt::arity == 2) {
        check_multirhs_pair<VectorType, VectorType, typename ctx::scalar, ctx::N>();
        if (!factorize_dispatch<UserConfig, oot_diagnostics>(A, piv, oot_count)) return false;
        substitute_multirhs_dispatch<ctx, pass_width, true>(A, piv, y, y);
        return true;
    } else {
        static_assert(pass_width == 0,
                      "tdls adaptors: pass_width only applies to matrix-like right-hand sides");
        check_vector<VectorType, typename ctx::scalar, ctx::N>();
        if constexpr (ctx::runtime_sized) {
            return ctx::dynamic_solver::template solve_inplace<oot_diagnostics>(
                mt::runtime_extent0(A), mt::pointer(A), mt::stride(A), pa::pointer(piv),
                pa::stride(piv), yt::pointer(y), yt::stride(y), oot_count);
        } else {
            return ctx::solver::template solve_inplace<yt::is_internal, pa::is_internal,
                                                       mt::is_internal, oot_diagnostics>(
                mt::pointer(A), mt::stride(A), pa::pointer(piv), pa::stride(piv), yt::pointer(y),
                yt::stride(y), oot_count);
        }
    }
}
 
} // namespace detail
 
 
 
template<detail::solver_config auto UserConfig, typename MatrixType, typename PivotType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool factorize(MatrixType& A,
                                                                         PivotType& piv) {
    int unused = 0;
    return detail::factorize_dispatch<UserConfig, false>(A, piv, unused);
}
 
template<typename MatrixType, typename PivotType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool factorize(MatrixType& A,
                                                                         PivotType& piv) {
    return factorize<TiledLUppConfig<detail::matrix_scalar<MatrixType>>{}>(A, piv);
}
 
template<detail::solver_config auto UserConfig, typename MatrixType, typename PivotType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
factorize(MatrixType& A, PivotType& piv, int& oot_count) {
    return detail::factorize_dispatch<UserConfig, true>(A, piv, oot_count);
}
 
template<typename MatrixType, typename PivotType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
factorize(MatrixType& A, PivotType& piv, int& oot_count) {
    return factorize<TiledLUppConfig<detail::matrix_scalar<MatrixType>>{}>(A, piv, oot_count);
}
 
template<detail::solver_config auto UserConfig, int pass_width = 0, typename MatrixType,
         typename PivotType, typename RhsType, typename SolutionType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
solve(MatrixType& A, PivotType& piv, const RhsType& b, SolutionType& x) {
    int unused = 0;
    return detail::solve_dispatch<UserConfig, pass_width, false>(A, piv, b, x, unused);
}
 
template<int pass_width = 0, typename MatrixType, typename PivotType, typename RhsType,
         typename SolutionType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
solve(MatrixType& A, PivotType& piv, const RhsType& b, SolutionType& x) {
    return solve<TiledLUppConfig<detail::matrix_scalar<MatrixType>>{}, pass_width>(A, piv, b, x);
}
 
template<detail::solver_config auto UserConfig, int pass_width = 0, typename MatrixType,
         typename PivotType, typename RhsType, typename SolutionType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
solve(MatrixType& A, PivotType& piv, const RhsType& b, SolutionType& x, int& oot_count) {
    return detail::solve_dispatch<UserConfig, pass_width, true>(A, piv, b, x, oot_count);
}
 
template<int pass_width = 0, typename MatrixType, typename PivotType, typename RhsType,
         typename SolutionType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
solve(MatrixType& A, PivotType& piv, const RhsType& b, SolutionType& x, int& oot_count) {
    return solve<TiledLUppConfig<detail::matrix_scalar<MatrixType>>{}, pass_width>(A, piv, b, x,
                                                                                   oot_count);
}
 
template<detail::solver_config auto UserConfig, int pass_width = 0, typename MatrixType,
         typename PivotType, typename VectorType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
solve_inplace(MatrixType& A, PivotType& piv, VectorType& y) {
    int unused = 0;
    return detail::solve_inplace_dispatch<UserConfig, pass_width, false>(A, piv, y, unused);
}
 
template<int pass_width = 0, typename MatrixType, typename PivotType, typename VectorType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
solve_inplace(MatrixType& A, PivotType& piv, VectorType& y) {
    return solve_inplace<TiledLUppConfig<detail::matrix_scalar<MatrixType>>{}, pass_width>(A, piv,
                                                                                           y);
}
 
template<detail::solver_config auto UserConfig, int pass_width = 0, typename MatrixType,
         typename PivotType, typename VectorType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
solve_inplace(MatrixType& A, PivotType& piv, VectorType& y, int& oot_count) {
    return detail::solve_inplace_dispatch<UserConfig, pass_width, true>(A, piv, y, oot_count);
}
 
template<int pass_width = 0, typename MatrixType, typename PivotType, typename VectorType>
[[nodiscard]] TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr bool
solve_inplace(MatrixType& A, PivotType& piv, VectorType& y, int& oot_count) {
    return solve_inplace<TiledLUppConfig<detail::matrix_scalar<MatrixType>>{}, pass_width>(
        A, piv, y, oot_count);
}
 
template<detail::solver_config auto UserConfig, int pass_width = 0, typename MatrixType,
         typename PivotType, typename RhsType, typename SolutionType>
TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr void
substitute(const MatrixType& A, const PivotType& piv, const RhsType& b, SolutionType& x) {
    using ctx =
        detail::adaptor_context<MatrixType, detail::checked_config<MatrixType, UserConfig>()>;
    using mt = typename ctx::mtraits;
    using pa = detail::pivot_access<PivotType, ctx::N>;
    using bt = storage_traits<std::remove_cv_t<RhsType>>;
    using xt = storage_traits<std::remove_cv_t<SolutionType>>;
    static_assert(xt::is_mutable, "tdls adaptors: substitute writes into x");
    static_assert(!std::is_const_v<SolutionType>,
                  "tdls adaptors: x must not be const here (substitute writes it)");
    if constexpr (bt::arity == 2) {
        detail::check_multirhs_pair<RhsType, SolutionType, typename ctx::scalar, ctx::N>();
        detail::substitute_multirhs_dispatch<ctx, pass_width, false>(A, piv, b, x);
    } else {
        static_assert(pass_width == 0,
                      "tdls adaptors: pass_width only applies to matrix-like right-hand sides");
        detail::check_vector<RhsType, typename ctx::scalar, ctx::N>();
        detail::check_vector<SolutionType, typename ctx::scalar, ctx::N>();
        detail::check_rhs_pair<RhsType, SolutionType>();
        if constexpr (ctx::runtime_sized) {
            ctx::dynamic_solver::substitute(mt::runtime_extent0(A), mt::pointer(A), mt::stride(A),
                                            pa::pointer(piv), pa::stride(piv), bt::pointer(b),
                                            xt::pointer(x), xt::stride(x));
        } else {
            ctx::solver::template substitute<xt::is_internal, pa::is_internal, mt::is_internal>(
                mt::pointer(A), mt::stride(A), pa::pointer(piv), pa::stride(piv), bt::pointer(b),
                xt::pointer(x), xt::stride(x));
        }
    }
}
 
template<int pass_width = 0, typename MatrixType, typename PivotType, typename RhsType,
         typename SolutionType>
TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr void
substitute(const MatrixType& A, const PivotType& piv, const RhsType& b, SolutionType& x) {
    substitute<TiledLUppConfig<detail::matrix_scalar<MatrixType>>{}, pass_width>(A, piv, b, x);
}
 
template<detail::solver_config auto UserConfig, int pass_width = 0, typename MatrixType,
         typename PivotType, typename SolutionType>
TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr void
substitute_inplace(const MatrixType& A, const PivotType& piv, SolutionType& x) {
    using ctx =
        detail::adaptor_context<MatrixType, detail::checked_config<MatrixType, UserConfig>()>;
    using mt = typename ctx::mtraits;
    using pa = detail::pivot_access<PivotType, ctx::N>;
    using xt = storage_traits<std::remove_cv_t<SolutionType>>;
    static_assert(xt::is_mutable, "tdls adaptors: substitute_inplace writes into x");
    static_assert(!std::is_const_v<SolutionType>,
                  "tdls adaptors: x must not be const here (substitute_inplace writes it)");
    if constexpr (xt::arity == 2) {
        detail::check_multirhs_pair<SolutionType, SolutionType, typename ctx::scalar, ctx::N>();
        detail::substitute_multirhs_dispatch<ctx, pass_width, true>(A, piv, x, x);
    } else {
        static_assert(pass_width == 0,
                      "tdls adaptors: pass_width only applies to matrix-like right-hand sides");
        detail::check_vector<SolutionType, typename ctx::scalar, ctx::N>();
        if constexpr (ctx::runtime_sized) {
            ctx::dynamic_solver::substitute_inplace(mt::runtime_extent0(A), mt::pointer(A),
                                                    mt::stride(A), pa::pointer(piv),
                                                    pa::stride(piv), xt::pointer(x), xt::stride(x));
        } else {
            ctx::solver::template substitute_inplace<xt::is_internal, pa::is_internal,
                                                     mt::is_internal>(
                mt::pointer(A), mt::stride(A), pa::pointer(piv), pa::stride(piv), xt::pointer(x),
                xt::stride(x));
        }
    }
}
 
template<int pass_width = 0, typename MatrixType, typename PivotType, typename SolutionType>
TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr void
substitute_inplace(const MatrixType& A, const PivotType& piv, SolutionType& x) {
    substitute_inplace<TiledLUppConfig<detail::matrix_scalar<MatrixType>>{}, pass_width>(A, piv, x);
}
 
template<detail::solver_config auto UserConfig, int pass_width = 0, typename MatrixType,
         typename PivotType, typename SolutionType>
TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr void
substitute_canonical(const MatrixType& A, const PivotType& piv, const int col, SolutionType& x) {
    using ctx =
        detail::adaptor_context<MatrixType, detail::checked_config<MatrixType, UserConfig>()>;
    using mt = typename ctx::mtraits;
    using pa = detail::pivot_access<PivotType, ctx::N>;
    using xt = storage_traits<std::remove_cv_t<SolutionType>>;
    static_assert(xt::is_mutable, "tdls adaptors: substitute_canonical writes into x");
    static_assert(!std::is_const_v<SolutionType>,
                  "tdls adaptors: x must not be const here (substitute_canonical writes it)");
    if constexpr (xt::arity == 2) {
        detail::check_multirhs_pair<SolutionType, SolutionType, typename ctx::scalar, ctx::N>();
        static_assert(pass_width >= 0, "tdls adaptors: pass_width must not be negative");
        if constexpr (ctx::runtime_sized) {
            ctx::dynamic_solver::template substitute_canonical_multirhs<pass_width>(
                mt::runtime_extent0(A), xt::columns(x), mt::pointer(A), mt::stride(A),
                pa::pointer(piv), pa::stride(piv), col, xt::pointer(x), xt::row_stride(x),
                xt::col_stride(x));
        } else {
            constexpr int M = xt::extent1;
            static_assert(M >= 1, "tdls adaptors: x must have at least one column");
            ctx::solver::template substitute_canonical_multirhs<M, false, pa::is_internal,
                                                                mt::is_internal, pass_width>(
                mt::pointer(A), mt::stride(A), pa::pointer(piv), pa::stride(piv), col,
                xt::pointer(x), xt::row_stride(x), xt::col_stride(x));
        }
    } else {
        static_assert(pass_width == 0,
                      "tdls adaptors: pass_width only applies to matrix-like right-hand sides");
        detail::check_vector<SolutionType, typename ctx::scalar, ctx::N>();
        if constexpr (ctx::runtime_sized) {
            ctx::dynamic_solver::substitute_canonical(
                mt::runtime_extent0(A), mt::pointer(A), mt::stride(A), pa::pointer(piv),
                pa::stride(piv), col, xt::pointer(x), xt::stride(x));
        } else {
            ctx::solver::template substitute_canonical<xt::is_internal, pa::is_internal,
                                                       mt::is_internal>(
                mt::pointer(A), mt::stride(A), pa::pointer(piv), pa::stride(piv), col,
                xt::pointer(x), xt::stride(x));
        }
    }
}
 
template<int pass_width = 0, typename MatrixType, typename PivotType, typename SolutionType>
TDLS_HOST_DEVICE TDLS_FORCEINLINE constexpr void
substitute_canonical(const MatrixType& A, const PivotType& piv, const int col, SolutionType& x) {
    substitute_canonical<TiledLUppConfig<detail::matrix_scalar<MatrixType>>{}, pass_width>(A, piv,
                                                                                           col, x);
}
 
 
 
} // namespace tdls
 
 
 
#endif // TDLS_TFEL_ADAPTORS_HPP