From d0387fb171c56d360cbd1d1c599581b8b08ea5b6 Mon Sep 17 00:00:00 2001 From: alonso130r <131227801+alonso130r@users.noreply.github.com> Date: Wed, 9 Sep 2026 15:24:33 -0400 Subject: [PATCH 1/3] added cholesky factorization/solve, need to reformat to escape karl's wrath --- xpu/linear_algebra.hpp | 593 ++++++++++++++++++++++++----------------- 1 file changed, 348 insertions(+), 245 deletions(-) diff --git a/xpu/linear_algebra.hpp b/xpu/linear_algebra.hpp index 38a918a..3e4dacf 100644 --- a/xpu/linear_algebra.hpp +++ b/xpu/linear_algebra.hpp @@ -7,10 +7,10 @@ #include #if defined(XPU_CUDA) - #include - #include +#include +#include #else - #include +#include #endif #include @@ -22,32 +22,24 @@ namespace xpu { namespace linalg { -enum class status { - success, - singular -}; +enum class status { success, singular, not_pd }; namespace detail { [[noreturn]] -inline auto linalg_error(const char* message) noexcept -> void { - std::fprintf( - stderr, - "xpu: %s\n", - message - ); +inline auto linalg_error(const char *message) noexcept -> void { + std::fprintf(stderr, "xpu: %s\n", message); std::abort(); } #if defined(XPU_CUDA) -template -struct build_identity { - T* matrix; +template struct build_identity { + T *matrix; std::size_t stride; DEVICE_ONLY - auto operator()(const xpu::array& index) const -> void { + auto operator()(const xpu::array &index) const -> void { const auto row{index[0]}; const auto column{index[1]}; @@ -59,13 +51,12 @@ struct build_identity { } }; -template -struct transpose_square { - T* matrix; +template struct transpose_square { + T *matrix; std::size_t stride; DEVICE_ONLY - auto operator()(const xpu::array& index) const -> void { + auto operator()(const xpu::array &index) const -> void { const auto row{index[0]}; const auto column{index[1]}; const auto skip_pair{row >= column}; @@ -91,11 +82,8 @@ inline auto create_cusolver_handle() -> cusolverDnHandle_t { } template -inline auto getrf_workspace_size( - cusolverDnHandle_t handle, - std::size_t order, - std::size_t stride -) -> std::size_t { +inline auto getrf_workspace_size(cusolverDnHandle_t handle, std::size_t order, + std::size_t stride) -> std::size_t { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; @@ -103,18 +91,10 @@ inline auto getrf_workspace_size( if constexpr (std::same_as) { xpu::cu_check(cusolverDnSgetrf_bufferSize( - handle, - vendor_order, vendor_order, - nullptr, vendor_stride, - &size - )); + handle, vendor_order, vendor_order, nullptr, vendor_stride, &size)); } else { xpu::cu_check(cusolverDnDgetrf_bufferSize( - handle, - vendor_order, vendor_order, - nullptr, vendor_stride, - &size - )); + handle, vendor_order, vendor_order, nullptr, vendor_stride, &size)); } const auto workspace_size{xpu::detail::checked_cast(size)}; @@ -123,170 +103,202 @@ inline auto getrf_workspace_size( } template -inline auto cusolver_getrf( - cusolverDnHandle_t handle, - std::size_t order, - std::size_t stride, - T* matrix, - T* workspace, - int* pivot, - int* info -) -> void { +inline auto potrf_workspace_size(cusolverDnHandle_t handle, std::size_t order, + std::size_t stride) -> std::size_t { + const auto vendor_order{xpu::detail::checked_cast(order)}; + const auto vendor_stride{xpu::detail::checked_cast(stride)}; + + auto size{0}; + + if constexpr (std::same_as) { + xpu::cu_check(cusolverDnSpotrf_bufferSize(handle, CUBLAS_FILL_MODE_UPPER, + vendor_order, nullptr, + vendor_stride, &size)); + } else { + xpu::cu_check(cusolverDnDpotrf_bufferSize(handle, CUBLAS_FILL_MODE_UPPER, + vendor_order, nullptr, + vendor_stride, &size)); + } + + const auto workspace_size{xpu::detail::checked_cast(size)}; + + return workspace_size; +} + +template +inline auto cusolver_getrf(cusolverDnHandle_t handle, std::size_t order, + std::size_t stride, T *matrix, T *workspace, + int *pivot, int *info) -> void { + const auto vendor_order{xpu::detail::checked_cast(order)}; + const auto vendor_stride{xpu::detail::checked_cast(stride)}; + + if constexpr (std::same_as) { + xpu::cu_check(cusolverDnSgetrf(handle, vendor_order, vendor_order, matrix, + vendor_stride, workspace, pivot, info)); + } else { + xpu::cu_check(cusolverDnDgetrf(handle, vendor_order, vendor_order, matrix, + vendor_stride, workspace, pivot, info)); + } +} + +template +inline auto cusolver_potrf(cusolverDnHandle_t handle, std::size_t order, + std::size_t stride, T *matrix, T *workspace, + std::size_t workspace_size, int *info) -> void { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; + const auto vendor_workspace_size{ + xpu::detail::checked_cast(workspace_size)}; if constexpr (std::same_as) { - xpu::cu_check(cusolverDnSgetrf( - handle, - vendor_order, vendor_order, - matrix, vendor_stride, - workspace, pivot, info - )); + xpu::cu_check(cusolverDnSpotrf(handle, CUBLAS_FILL_MODE_UPPER, vendor_order, + matrix, vendor_stride, workspace, + vendor_workspace_size, info)); } else { - xpu::cu_check(cusolverDnDgetrf( - handle, - vendor_order, vendor_order, - matrix, vendor_stride, - workspace, pivot, info - )); + xpu::cu_check(cusolverDnDpotrf(handle, CUBLAS_FILL_MODE_UPPER, vendor_order, + matrix, vendor_stride, workspace, + vendor_workspace_size, info)); } } template -inline auto cusolver_getrs( - cusolverDnHandle_t handle, - cublasOperation_t operation, - std::size_t order, - std::size_t right_hand_sides, - const T* lower_upper, - std::size_t lower_upper_stride, - const int* pivot, - T* solution, - std::size_t solution_stride, - int* info -) -> void { +inline auto cusolver_getrs(cusolverDnHandle_t handle, + cublasOperation_t operation, std::size_t order, + std::size_t right_hand_sides, const T *lower_upper, + std::size_t lower_upper_stride, const int *pivot, + T *solution, std::size_t solution_stride, int *info) + -> void { const auto vendor_order{xpu::detail::checked_cast(order)}; - const auto vendor_right_hand_sides{xpu::detail::checked_cast(right_hand_sides)}; - const auto vendor_lower_upper_stride{xpu::detail::checked_cast(lower_upper_stride)}; - const auto vendor_solution_stride{xpu::detail::checked_cast(solution_stride)}; + const auto vendor_right_hand_sides{ + xpu::detail::checked_cast(right_hand_sides)}; + const auto vendor_lower_upper_stride{ + xpu::detail::checked_cast(lower_upper_stride)}; + const auto vendor_solution_stride{ + xpu::detail::checked_cast(solution_stride)}; if constexpr (std::same_as) { - xpu::cu_check(cusolverDnSgetrs( - handle, operation, - vendor_order, vendor_right_hand_sides, - lower_upper, vendor_lower_upper_stride, - pivot, - solution, vendor_solution_stride, - info - )); + xpu::cu_check(cusolverDnSgetrs(handle, operation, vendor_order, + vendor_right_hand_sides, lower_upper, + vendor_lower_upper_stride, pivot, solution, + vendor_solution_stride, info)); } else { - xpu::cu_check(cusolverDnDgetrs( - handle, operation, - vendor_order, vendor_right_hand_sides, - lower_upper, vendor_lower_upper_stride, - pivot, - solution, vendor_solution_stride, - info - )); + xpu::cu_check(cusolverDnDgetrs(handle, operation, vendor_order, + vendor_right_hand_sides, lower_upper, + vendor_lower_upper_stride, pivot, solution, + vendor_solution_stride, info)); + } +} + +template +inline auto cusolver_potrs(cusolverDnHandle_t handle, std::size_t order, + std::size_t right_hand_sides, const T *factor, + std::size_t factor_stride, T *solution, + std::size_t solution_stride, int *info) -> void { + const auto vendor_order{xpu::detail::checked_cast(order)}; + const auto vendor_right_hand_sides{ + xpu::detail::checked_cast(right_hand_sides)}; + const auto vendor_factor_stride{ + xpu::detail::checked_cast(factor_stride)}; + const auto vendor_solution_stride{ + xpu::detail::checked_cast(solution_stride)}; + + if constexpr (std::same_as) { + xpu::cu_check(cusolverDnSpotrs( + handle, CUBLAS_FILL_MODE_UPPER, vendor_order, vendor_right_hand_sides, + factor, vendor_factor_stride, solution, vendor_solution_stride, info)); + } else { + xpu::cu_check(cusolverDnDpotrs( + handle, CUBLAS_FILL_MODE_UPPER, vendor_order, vendor_right_hand_sides, + factor, vendor_factor_stride, solution, vendor_solution_stride, info)); } } #else template -inline auto lapacke_getrf( - T* RESTRICT matrix, - lapack_int* RESTRICT pivot, - std::size_t order, - std::size_t stride -) -> lapack_int { +inline auto lapacke_getrf(T *RESTRICT matrix, lapack_int *RESTRICT pivot, + std::size_t order, std::size_t stride) -> lapack_int { + const auto vendor_order{xpu::detail::checked_cast(order)}; + const auto vendor_stride{xpu::detail::checked_cast(stride)}; + + if constexpr (std::same_as) { + return LAPACKE_sgetrf(LAPACK_ROW_MAJOR, vendor_order, vendor_order, matrix, + vendor_stride, pivot); + } else { + return LAPACKE_dgetrf(LAPACK_ROW_MAJOR, vendor_order, vendor_order, matrix, + vendor_stride, pivot); + } +} + +template +inline auto lapacke_potrf(T *RESTRICT matrix, std::size_t order, + std::size_t stride) -> lapack_int { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; if constexpr (std::same_as) { - return LAPACKE_sgetrf( - LAPACK_ROW_MAJOR, - vendor_order, vendor_order, - matrix, vendor_stride, - pivot - ); + return LAPACKE_spotrf(LAPACK_ROW_MAJOR, 'L', vendor_order, matrix, + vendor_stride); // use lower triangular } else { - return LAPACKE_dgetrf( - LAPACK_ROW_MAJOR, - vendor_order, vendor_order, - matrix, vendor_stride, - pivot - ); + return LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'L', vendor_order, matrix, + vendor_stride); } } template -inline auto lapacke_getri( - T* RESTRICT inverse, - const lapack_int* RESTRICT pivot, - std::size_t order, - std::size_t stride -) -> lapack_int { +inline auto lapacke_getri(T *RESTRICT inverse, const lapack_int *RESTRICT pivot, + std::size_t order, std::size_t stride) -> lapack_int { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; if constexpr (std::same_as) { - return LAPACKE_sgetri( - LAPACK_ROW_MAJOR, - vendor_order, - inverse, vendor_stride, - pivot - ); + return LAPACKE_sgetri(LAPACK_ROW_MAJOR, vendor_order, inverse, + vendor_stride, pivot); } else { - return LAPACKE_dgetri( - LAPACK_ROW_MAJOR, - vendor_order, - inverse, vendor_stride, - pivot - ); + return LAPACKE_dgetri(LAPACK_ROW_MAJOR, vendor_order, inverse, + vendor_stride, pivot); } } template -inline auto lapacke_getrs( - const T* RESTRICT lower_upper, - const lapack_int* RESTRICT pivot, - T* RESTRICT solution, - std::size_t order, - std::size_t stride -) -> lapack_int { +inline auto lapacke_getrs(const T *RESTRICT lower_upper, + const lapack_int *RESTRICT pivot, + T *RESTRICT solution, std::size_t order, + std::size_t stride) -> lapack_int { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; if constexpr (std::same_as) { - return LAPACKE_sgetrs( - LAPACK_ROW_MAJOR, 'N', - vendor_order, 1, - lower_upper, vendor_stride, - pivot, - solution, 1 - ); + return LAPACKE_sgetrs(LAPACK_ROW_MAJOR, 'N', vendor_order, 1, lower_upper, + vendor_stride, pivot, solution, 1); } else { - return LAPACKE_dgetrs( - LAPACK_ROW_MAJOR, 'N', - vendor_order, 1, - lower_upper, vendor_stride, - pivot, - solution, 1 - ); + return LAPACKE_dgetrs(LAPACK_ROW_MAJOR, 'N', vendor_order, 1, lower_upper, + vendor_stride, pivot, solution, 1); + } +} + +template +inline auto lapacke_potrs(const T *RESTRICT lower_upper, T *RESTRICT solution, + std::size_t order, std::size_t stride) -> lapack_int { + const auto vendor_order{xpu::detail::checked_cast(order)}; + const auto vendor_stride{xpu::detail::checked_cast(stride)}; + + if constexpr (std::same_as) { + return LAPACKE_spotrs(LAPACK_ROW_MAJOR, 'L', vendor_order, 1, lower_upper, + vendor_stride, solution, 1); + } else { + return LAPACKE_dpotrs(LAPACK_ROW_MAJOR, 'L', vendor_order, 1, lower_upper, + vendor_stride, solution, 1); } } #endif -} // namespace xpu::linalg::detail +} // namespace detail template -inline auto transpose_square( - T* RESTRICT matrix, - std::size_t order, - std::size_t stride -) noexcept -> void { +inline auto transpose_square(T *RESTRICT matrix, std::size_t order, + std::size_t stride) noexcept -> void { const auto empty_matrix{order == 0uz}; if (empty_matrix) { @@ -305,16 +317,9 @@ inline auto transpose_square( static_cast(xpu::detail::checked_bytes(matrix_size)); #if defined(XPU_CUDA) - const auto range = xpu::range<2uz>{ - {0uz, 0uz}, - {order, order}, - {1uz, 1uz} - }; + const auto range = xpu::range<2uz>{{0uz, 0uz}, {order, order}, {1uz, 1uz}}; - const auto transpose = detail::transpose_square{ - matrix, - stride - }; + const auto transpose = detail::transpose_square{matrix, stride}; xpu::parallel_for(range, transpose); #else @@ -323,21 +328,152 @@ inline auto transpose_square( const auto col_idx{row * stride + column}; const auto row_idx{column * stride + row}; - std::swap( - matrix[col_idx], - matrix[row_idx] - ); + std::swap(matrix[col_idx], matrix[row_idx]); } } #endif } -template -class lu_factorization { +template class cholesky_factorization { +private: + std::size_t order_; + std::size_t stride_; + T *factor_{}; + +#if defined(XPU_CUDA) + using dimension_type = int; +#else + using dimension_type = lapack_int; +#endif + +#if defined(XPU_CUDA) + cusolverDnHandle_t handle_; + xpu::buffer workspace_; + xpu::buffer info_; +#endif + + [[nodiscard]] + static auto checked_order(std::size_t order, std::size_t stride) noexcept + -> std::size_t { + static_cast(xpu::detail::checked_cast(order)); + static_cast(xpu::detail::checked_cast(stride)); + + const auto invalid_stride{stride < order}; + + if (invalid_stride) { + detail::linalg_error("matrix stride is smaller than its order"); + } + + const auto matrix_size{xpu::detail::checked_mul(order, stride)}; + + static_cast(xpu::detail::checked_bytes(matrix_size)); + + return order; + } + +public: + cholesky_factorization(std::size_t order, std::size_t stride) + : order_{checked_order(order, stride)}, stride_{stride} +#if defined(XPU_CUDA) + , + handle_{detail::create_cusolver_handle()}, + workspace_{detail::potrf_workspace_size(handle_, order, stride)}, + info_{1uz} +#endif + { + } + + ~cholesky_factorization() { +#if defined(XPU_CUDA) + xpu::cu_check(cusolverDnDestroy(handle_)); +#endif + } + + [[nodiscard]] + constexpr auto order() const noexcept -> std::size_t { + return order_; + } + + [[nodiscard]] + constexpr auto stride() const noexcept -> std::size_t { + return stride_; + } + + [[nodiscard]] + auto factorize(T *RESTRICT matrix) noexcept -> status { + factor_ = nullptr; + +#if defined(XPU_CUDA) + detail::cusolver_potrf(handle_, order_, stride_, matrix, workspace_.data(), + workspace_.size(), info_.data()); + + auto info{0}; + xpu::copy_n(&info, info_.data(), 1uz); + if (info < 0) { + detail::linalg_error("cuSOLVER potrf received an invalid argument"); + } + if (info > 0) { + return status::not_pd; + } +#else + const auto info{detail::lapacke_potrf(matrix, order_, stride_)}; + if (info < 0) { + detail::linalg_error("LAPACKE potrf failed"); + } + if (info > 0) { + return status::not_pd; + } +#endif + + factor_ = matrix; + return status::success; + } + + auto solve(const T *RESTRICT factor, const T *RESTRICT rhs, + T *RESTRICT solution) noexcept -> void { + if (!factor_) { + detail::linalg_error("you must factor before you solve"); + } + if (factor != factor_) { + detail::linalg_error( + "solve requires the most recently factorized matrix"); + } + if (rhs == solution) { + detail::linalg_error("rhs and solution must not alias"); + } + + xpu::copy_n(solution, rhs, order_); + +#if defined(XPU_CUDA) + detail::cusolver_potrs(handle_, order_, 1uz, factor, stride_, solution, + order_, info_.data()); + + auto info{0}; + xpu::copy_n(&info, info_.data(), 1uz); + if (info != 0) { + detail::linalg_error("cuSOLVER potrs received an invalid argument"); + } +#else + const auto info{detail::lapacke_potrs(factor, solution, order_, stride_)}; + if (info != 0) { + detail::linalg_error("LAPACKE potrs failed"); + } +#endif + } + + auto operator=(const cholesky_factorization &) + -> cholesky_factorization & = delete; + cholesky_factorization(const cholesky_factorization &) = delete; + auto operator=(cholesky_factorization &&) + -> cholesky_factorization & = delete; + cholesky_factorization(cholesky_factorization &&) = delete; +}; + +template class lu_factorization { private: std::size_t order_; std::size_t stride_; - T* lower_upper_{}; + T *lower_upper_{}; #if defined(XPU_CUDA) using dimension_type = int; @@ -355,7 +491,8 @@ class lu_factorization { #endif [[nodiscard]] - static auto checked_order(std::size_t order, std::size_t stride) noexcept -> std::size_t { + static auto checked_order(std::size_t order, std::size_t stride) noexcept + -> std::size_t { static_cast(xpu::detail::checked_cast(order)); static_cast(xpu::detail::checked_cast(stride)); @@ -374,17 +511,18 @@ class lu_factorization { public: lu_factorization(std::size_t order, std::size_t stride) - : order_{checked_order(order, stride)} - , stride_{stride} + : order_{checked_order(order, stride)}, stride_{stride} #if defined(XPU_CUDA) - , pivot_{order} - , handle_{detail::create_cusolver_handle()} - , workspace_{detail::getrf_workspace_size(handle_, order, stride)} - , info_{1uz} + , + pivot_{order}, handle_{detail::create_cusolver_handle()}, + workspace_{detail::getrf_workspace_size(handle_, order, stride)}, + info_{1uz} #else - , pivot_{order} + , + pivot_{order} #endif - { } + { + } ~lu_factorization() { #if defined(XPU_CUDA) @@ -403,46 +541,41 @@ class lu_factorization { } [[nodiscard]] - auto factorize(T* RESTRICT matrix) noexcept -> status { + auto factorize(T *RESTRICT matrix) noexcept -> status { lower_upper_ = nullptr; #if defined(XPU_CUDA) - detail::cusolver_getrf( - handle_, order_, stride_, - matrix, workspace_.data(), pivot_.data(), info_.data() - ); + detail::cusolver_getrf(handle_, order_, stride_, matrix, workspace_.data(), + pivot_.data(), info_.data()); auto info{0}; xpu::copy_n(&info, info_.data(), 1uz); if (info < 0) { detail::linalg_error("cuSOLVER getrf received an invalid argument"); } - if (info > 0) { return status::singular; } + if (info > 0) { + return status::singular; + } #else const auto info{ - detail::lapacke_getrf( - matrix, pivot_.data(), order_, stride_ - ) - }; + detail::lapacke_getrf(matrix, pivot_.data(), order_, stride_)}; if (info < 0) { detail::linalg_error("LAPACKE getrf failed"); } - if (info > 0) { return status::singular; } + if (info > 0) { + return status::singular; + } #endif lower_upper_ = matrix; return status::success; } - auto solve( - const T* RESTRICT lower_upper, - const T* RESTRICT rhs, - T* RESTRICT solution - ) noexcept -> void { + auto solve(const T *RESTRICT lower_upper, const T *RESTRICT rhs, + T *RESTRICT solution) noexcept -> void { if (lower_upper != lower_upper_) { detail::linalg_error( - "solve requires the most recently factorized matrix" - ); + "solve requires the most recently factorized matrix"); } if (rhs == solution) { detail::linalg_error("rhs and solution must not alias"); @@ -451,14 +584,9 @@ class lu_factorization { xpu::copy_n(solution, rhs, order_); #if defined(XPU_CUDA) - detail::cusolver_getrs( - handle_, CUBLAS_OP_T, - order_, 1uz, - lower_upper, stride_, - pivot_.data(), - solution, order_, - info_.data() - ); + detail::cusolver_getrs(handle_, CUBLAS_OP_T, order_, 1uz, lower_upper, + stride_, pivot_.data(), solution, order_, + info_.data()); auto info{0}; xpu::copy_n(&info, info_.data(), 1uz); @@ -466,53 +594,35 @@ class lu_factorization { detail::linalg_error("cuSOLVER getrs received an invalid argument"); } #else - const auto info{ - detail::lapacke_getrs( - lower_upper, pivot_.data(), solution, - order_, stride_ - ) - }; + const auto info{detail::lapacke_getrs(lower_upper, pivot_.data(), solution, + order_, stride_)}; if (info != 0) { detail::linalg_error("LAPACKE getrs failed"); } #endif } - auto invert( - const T* RESTRICT lower_upper, - T* RESTRICT inverse - ) noexcept -> void { + auto invert(const T *RESTRICT lower_upper, T *RESTRICT inverse) noexcept + -> void { if (lower_upper != lower_upper_) { detail::linalg_error( - "invert requires the most recently factorized matrix" - ); + "invert requires the most recently factorized matrix"); } if (lower_upper == inverse) { detail::linalg_error("lower_upper and inverse must not alias"); } #if defined(XPU_CUDA) - const auto range = xpu::range<2uz>{ - {0uz, 0uz}, - {order_, order_}, - {1uz, 1uz} - }; + const auto range = + xpu::range<2uz>{{0uz, 0uz}, {order_, order_}, {1uz, 1uz}}; - const auto initialize = detail::build_identity{ - inverse, - stride_ - }; + const auto initialize = detail::build_identity{inverse, stride_}; xpu::parallel_for(range, initialize); - detail::cusolver_getrs( - handle_, CUBLAS_OP_N, - order_, order_, - lower_upper, stride_, - pivot_.data(), - inverse, stride_, - info_.data() - ); + detail::cusolver_getrs(handle_, CUBLAS_OP_N, order_, order_, lower_upper, + stride_, pivot_.data(), inverse, stride_, + info_.data()); auto info{0}; xpu::copy_n(&info, info_.data(), 1uz); @@ -521,30 +631,23 @@ class lu_factorization { } #else for (auto row{0uz}; row < order_; ++row) { - xpu::copy_n( - inverse + row * stride_, - lower_upper + row * stride_, - order_ - ); + xpu::copy_n(inverse + row * stride_, lower_upper + row * stride_, order_); } const auto info{ - detail::lapacke_getri( - inverse, pivot_.data(), order_, stride_ - ) - }; + detail::lapacke_getri(inverse, pivot_.data(), order_, stride_)}; if (info != 0) { detail::linalg_error("LAPACKE getri failed"); } #endif } - auto operator=(const lu_factorization&) -> lu_factorization& = delete; - lu_factorization(const lu_factorization&) = delete; - auto operator=(lu_factorization&&) -> lu_factorization& = delete; - lu_factorization(lu_factorization&&) = delete; + auto operator=(const lu_factorization &) -> lu_factorization & = delete; + lu_factorization(const lu_factorization &) = delete; + auto operator=(lu_factorization &&) -> lu_factorization & = delete; + lu_factorization(lu_factorization &&) = delete; }; -} // namespace xpu::linalg +} // namespace linalg } // namespace xpu From 77831aa5eff1a73f6d514631de3930c045fae1c3 Mon Sep 17 00:00:00 2001 From: alonso130r <131227801+alonso130r@users.noreply.github.com> Date: Wed, 9 Sep 2026 15:35:13 -0400 Subject: [PATCH 2/3] escaping karl's wrath... --- xpu/linear_algebra.hpp | 579 +++++++++++++++++++++++++++-------------- 1 file changed, 385 insertions(+), 194 deletions(-) diff --git a/xpu/linear_algebra.hpp b/xpu/linear_algebra.hpp index 3e4dacf..af246b4 100644 --- a/xpu/linear_algebra.hpp +++ b/xpu/linear_algebra.hpp @@ -7,10 +7,10 @@ #include #if defined(XPU_CUDA) -#include -#include + #include + #include #else -#include + #include #endif #include @@ -22,24 +22,33 @@ namespace xpu { namespace linalg { -enum class status { success, singular, not_pd }; +enum class status { + success, + singular, + not_pd +}; namespace detail { [[noreturn]] -inline auto linalg_error(const char *message) noexcept -> void { - std::fprintf(stderr, "xpu: %s\n", message); +inline auto linalg_error(const char* message) noexcept -> void { + std::fprintf( + stderr, + "xpu: %s\n", + message + ); std::abort(); } #if defined(XPU_CUDA) -template struct build_identity { - T *matrix; +template +struct build_identity { + T* matrix; std::size_t stride; DEVICE_ONLY - auto operator()(const xpu::array &index) const -> void { + auto operator()(const xpu::array& index) const -> void { const auto row{index[0]}; const auto column{index[1]}; @@ -51,12 +60,13 @@ template struct build_identity { } }; -template struct transpose_square { - T *matrix; +template +struct transpose_square { + T* matrix; std::size_t stride; DEVICE_ONLY - auto operator()(const xpu::array &index) const -> void { + auto operator()(const xpu::array& index) const -> void { const auto row{index[0]}; const auto column{index[1]}; const auto skip_pair{row >= column}; @@ -82,8 +92,11 @@ inline auto create_cusolver_handle() -> cusolverDnHandle_t { } template -inline auto getrf_workspace_size(cusolverDnHandle_t handle, std::size_t order, - std::size_t stride) -> std::size_t { +inline auto getrf_workspace_size( + cusolverDnHandle_t handle, + std::size_t order, + std::size_t stride +) -> std::size_t { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; @@ -91,10 +104,18 @@ inline auto getrf_workspace_size(cusolverDnHandle_t handle, std::size_t order, if constexpr (std::same_as) { xpu::cu_check(cusolverDnSgetrf_bufferSize( - handle, vendor_order, vendor_order, nullptr, vendor_stride, &size)); + handle, + vendor_order, vendor_order, + nullptr, vendor_stride, + &size + )); } else { xpu::cu_check(cusolverDnDgetrf_bufferSize( - handle, vendor_order, vendor_order, nullptr, vendor_stride, &size)); + handle, + vendor_order, vendor_order, + nullptr, vendor_stride, + &size + )); } const auto workspace_size{xpu::detail::checked_cast(size)}; @@ -103,21 +124,30 @@ inline auto getrf_workspace_size(cusolverDnHandle_t handle, std::size_t order, } template -inline auto potrf_workspace_size(cusolverDnHandle_t handle, std::size_t order, - std::size_t stride) -> std::size_t { +inline auto potrf_workspace_size( + cusolverDnHandle_t handle, + std::size_t order, + std::size_t stride +) -> std::size_t { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; auto size{0}; if constexpr (std::same_as) { - xpu::cu_check(cusolverDnSpotrf_bufferSize(handle, CUBLAS_FILL_MODE_UPPER, - vendor_order, nullptr, - vendor_stride, &size)); + xpu::cu_check(cusolverDnSpotrf_bufferSize( + handle, CUBLAS_FILL_MODE_UPPER, + vendor_order, + nullptr, vendor_stride, + &size + )); } else { - xpu::cu_check(cusolverDnDpotrf_bufferSize(handle, CUBLAS_FILL_MODE_UPPER, - vendor_order, nullptr, - vendor_stride, &size)); + xpu::cu_check(cusolverDnDpotrf_bufferSize( + handle, CUBLAS_FILL_MODE_UPPER, + vendor_order, + nullptr, vendor_stride, + &size + )); } const auto workspace_size{xpu::detail::checked_cast(size)}; @@ -126,179 +156,289 @@ inline auto potrf_workspace_size(cusolverDnHandle_t handle, std::size_t order, } template -inline auto cusolver_getrf(cusolverDnHandle_t handle, std::size_t order, - std::size_t stride, T *matrix, T *workspace, - int *pivot, int *info) -> void { +inline auto cusolver_getrf( + cusolverDnHandle_t handle, + std::size_t order, + std::size_t stride, + T* matrix, + T* workspace, + int* pivot, + int* info +) -> void { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; if constexpr (std::same_as) { - xpu::cu_check(cusolverDnSgetrf(handle, vendor_order, vendor_order, matrix, - vendor_stride, workspace, pivot, info)); + xpu::cu_check(cusolverDnSgetrf( + handle, + vendor_order, vendor_order, + matrix, vendor_stride, + workspace, pivot, info + )); } else { - xpu::cu_check(cusolverDnDgetrf(handle, vendor_order, vendor_order, matrix, - vendor_stride, workspace, pivot, info)); + xpu::cu_check(cusolverDnDgetrf( + handle, + vendor_order, vendor_order, + matrix, vendor_stride, + workspace, pivot, info + )); } } template -inline auto cusolver_potrf(cusolverDnHandle_t handle, std::size_t order, - std::size_t stride, T *matrix, T *workspace, - std::size_t workspace_size, int *info) -> void { +inline auto cusolver_potrf( + cusolverDnHandle_t handle, + std::size_t order, + std::size_t stride, + T* matrix, + T* workspace, + std::size_t workspace_size, + int* info +) -> void { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; - const auto vendor_workspace_size{ - xpu::detail::checked_cast(workspace_size)}; + const auto vendor_workspace_size{xpu::detail::checked_cast(workspace_size)}; if constexpr (std::same_as) { - xpu::cu_check(cusolverDnSpotrf(handle, CUBLAS_FILL_MODE_UPPER, vendor_order, - matrix, vendor_stride, workspace, - vendor_workspace_size, info)); + xpu::cu_check(cusolverDnSpotrf( + handle, CUBLAS_FILL_MODE_UPPER, + vendor_order, + matrix, vendor_stride, + workspace, vendor_workspace_size, + info + )); } else { - xpu::cu_check(cusolverDnDpotrf(handle, CUBLAS_FILL_MODE_UPPER, vendor_order, - matrix, vendor_stride, workspace, - vendor_workspace_size, info)); + xpu::cu_check(cusolverDnDpotrf( + handle, CUBLAS_FILL_MODE_UPPER, + vendor_order, + matrix, vendor_stride, + workspace, vendor_workspace_size, + info + )); } } template -inline auto cusolver_getrs(cusolverDnHandle_t handle, - cublasOperation_t operation, std::size_t order, - std::size_t right_hand_sides, const T *lower_upper, - std::size_t lower_upper_stride, const int *pivot, - T *solution, std::size_t solution_stride, int *info) - -> void { +inline auto cusolver_getrs( + cusolverDnHandle_t handle, + cublasOperation_t operation, + std::size_t order, + std::size_t right_hand_sides, + const T* lower_upper, + std::size_t lower_upper_stride, + const int* pivot, + T* solution, + std::size_t solution_stride, + int* info +) -> void { const auto vendor_order{xpu::detail::checked_cast(order)}; - const auto vendor_right_hand_sides{ - xpu::detail::checked_cast(right_hand_sides)}; - const auto vendor_lower_upper_stride{ - xpu::detail::checked_cast(lower_upper_stride)}; - const auto vendor_solution_stride{ - xpu::detail::checked_cast(solution_stride)}; + const auto vendor_right_hand_sides{xpu::detail::checked_cast(right_hand_sides)}; + const auto vendor_lower_upper_stride{xpu::detail::checked_cast(lower_upper_stride)}; + const auto vendor_solution_stride{xpu::detail::checked_cast(solution_stride)}; if constexpr (std::same_as) { - xpu::cu_check(cusolverDnSgetrs(handle, operation, vendor_order, - vendor_right_hand_sides, lower_upper, - vendor_lower_upper_stride, pivot, solution, - vendor_solution_stride, info)); + xpu::cu_check(cusolverDnSgetrs( + handle, operation, + vendor_order, vendor_right_hand_sides, + lower_upper, vendor_lower_upper_stride, + pivot, + solution, vendor_solution_stride, + info + )); } else { - xpu::cu_check(cusolverDnDgetrs(handle, operation, vendor_order, - vendor_right_hand_sides, lower_upper, - vendor_lower_upper_stride, pivot, solution, - vendor_solution_stride, info)); + xpu::cu_check(cusolverDnDgetrs( + handle, operation, + vendor_order, vendor_right_hand_sides, + lower_upper, vendor_lower_upper_stride, + pivot, + solution, vendor_solution_stride, + info + )); } } template -inline auto cusolver_potrs(cusolverDnHandle_t handle, std::size_t order, - std::size_t right_hand_sides, const T *factor, - std::size_t factor_stride, T *solution, - std::size_t solution_stride, int *info) -> void { +inline auto cusolver_potrs( + cusolverDnHandle_t handle, + std::size_t order, + std::size_t right_hand_sides, + const T* factor, + std::size_t factor_stride, + T* solution, + std::size_t solution_stride, + int* info +) -> void { const auto vendor_order{xpu::detail::checked_cast(order)}; - const auto vendor_right_hand_sides{ - xpu::detail::checked_cast(right_hand_sides)}; - const auto vendor_factor_stride{ - xpu::detail::checked_cast(factor_stride)}; - const auto vendor_solution_stride{ - xpu::detail::checked_cast(solution_stride)}; + const auto vendor_right_hand_sides{xpu::detail::checked_cast(right_hand_sides)}; + const auto vendor_factor_stride{xpu::detail::checked_cast(factor_stride)}; + const auto vendor_solution_stride{xpu::detail::checked_cast(solution_stride)}; if constexpr (std::same_as) { xpu::cu_check(cusolverDnSpotrs( - handle, CUBLAS_FILL_MODE_UPPER, vendor_order, vendor_right_hand_sides, - factor, vendor_factor_stride, solution, vendor_solution_stride, info)); + handle, CUBLAS_FILL_MODE_UPPER, + vendor_order, vendor_right_hand_sides, + factor, vendor_factor_stride, + solution, vendor_solution_stride, + info + )); } else { xpu::cu_check(cusolverDnDpotrs( - handle, CUBLAS_FILL_MODE_UPPER, vendor_order, vendor_right_hand_sides, - factor, vendor_factor_stride, solution, vendor_solution_stride, info)); + handle, CUBLAS_FILL_MODE_UPPER, + vendor_order, vendor_right_hand_sides, + factor, vendor_factor_stride, + solution, vendor_solution_stride, + info + )); } } #else template -inline auto lapacke_getrf(T *RESTRICT matrix, lapack_int *RESTRICT pivot, - std::size_t order, std::size_t stride) -> lapack_int { +inline auto lapacke_getrf( + T* RESTRICT matrix, + lapack_int* RESTRICT pivot, + std::size_t order, + std::size_t stride +) -> lapack_int { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; if constexpr (std::same_as) { - return LAPACKE_sgetrf(LAPACK_ROW_MAJOR, vendor_order, vendor_order, matrix, - vendor_stride, pivot); + return LAPACKE_sgetrf( + LAPACK_ROW_MAJOR, + vendor_order, vendor_order, + matrix, vendor_stride, + pivot + ); } else { - return LAPACKE_dgetrf(LAPACK_ROW_MAJOR, vendor_order, vendor_order, matrix, - vendor_stride, pivot); + return LAPACKE_dgetrf( + LAPACK_ROW_MAJOR, + vendor_order, vendor_order, + matrix, vendor_stride, + pivot + ); } } template -inline auto lapacke_potrf(T *RESTRICT matrix, std::size_t order, - std::size_t stride) -> lapack_int { +inline auto lapacke_potrf( + T* RESTRICT matrix, + std::size_t order, + std::size_t stride +) -> lapack_int { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; if constexpr (std::same_as) { - return LAPACKE_spotrf(LAPACK_ROW_MAJOR, 'L', vendor_order, matrix, - vendor_stride); // use lower triangular + return LAPACKE_spotrf( + LAPACK_ROW_MAJOR, 'L', + vendor_order, + matrix, vendor_stride + ); // use lower triangular } else { - return LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'L', vendor_order, matrix, - vendor_stride); + return LAPACKE_dpotrf( + LAPACK_ROW_MAJOR, 'L', + vendor_order, + matrix, vendor_stride + ); } } template -inline auto lapacke_getri(T *RESTRICT inverse, const lapack_int *RESTRICT pivot, - std::size_t order, std::size_t stride) -> lapack_int { +inline auto lapacke_getri( + T* RESTRICT inverse, + const lapack_int* RESTRICT pivot, + std::size_t order, + std::size_t stride +) -> lapack_int { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; if constexpr (std::same_as) { - return LAPACKE_sgetri(LAPACK_ROW_MAJOR, vendor_order, inverse, - vendor_stride, pivot); + return LAPACKE_sgetri( + LAPACK_ROW_MAJOR, + vendor_order, + inverse, vendor_stride, + pivot + ); } else { - return LAPACKE_dgetri(LAPACK_ROW_MAJOR, vendor_order, inverse, - vendor_stride, pivot); + return LAPACKE_dgetri( + LAPACK_ROW_MAJOR, + vendor_order, + inverse, vendor_stride, + pivot + ); } } template -inline auto lapacke_getrs(const T *RESTRICT lower_upper, - const lapack_int *RESTRICT pivot, - T *RESTRICT solution, std::size_t order, - std::size_t stride) -> lapack_int { +inline auto lapacke_getrs( + const T* RESTRICT lower_upper, + const lapack_int* RESTRICT pivot, + T* RESTRICT solution, + std::size_t order, + std::size_t stride +) -> lapack_int { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; if constexpr (std::same_as) { - return LAPACKE_sgetrs(LAPACK_ROW_MAJOR, 'N', vendor_order, 1, lower_upper, - vendor_stride, pivot, solution, 1); + return LAPACKE_sgetrs( + LAPACK_ROW_MAJOR, 'N', + vendor_order, 1, + lower_upper, vendor_stride, + pivot, + solution, 1 + ); } else { - return LAPACKE_dgetrs(LAPACK_ROW_MAJOR, 'N', vendor_order, 1, lower_upper, - vendor_stride, pivot, solution, 1); + return LAPACKE_dgetrs( + LAPACK_ROW_MAJOR, 'N', + vendor_order, 1, + lower_upper, vendor_stride, + pivot, + solution, 1 + ); } } template -inline auto lapacke_potrs(const T *RESTRICT lower_upper, T *RESTRICT solution, - std::size_t order, std::size_t stride) -> lapack_int { +inline auto lapacke_potrs( + const T* RESTRICT lower_upper, + T* RESTRICT solution, + std::size_t order, + std::size_t stride +) -> lapack_int { const auto vendor_order{xpu::detail::checked_cast(order)}; const auto vendor_stride{xpu::detail::checked_cast(stride)}; if constexpr (std::same_as) { - return LAPACKE_spotrs(LAPACK_ROW_MAJOR, 'L', vendor_order, 1, lower_upper, - vendor_stride, solution, 1); + return LAPACKE_spotrs( + LAPACK_ROW_MAJOR, 'L', + vendor_order, 1, + lower_upper, vendor_stride, + solution, 1 + ); } else { - return LAPACKE_dpotrs(LAPACK_ROW_MAJOR, 'L', vendor_order, 1, lower_upper, - vendor_stride, solution, 1); + return LAPACKE_dpotrs( + LAPACK_ROW_MAJOR, 'L', + vendor_order, 1, + lower_upper, vendor_stride, + solution, 1 + ); } } #endif -} // namespace detail +} // namespace xpu::linalg::detail template -inline auto transpose_square(T *RESTRICT matrix, std::size_t order, - std::size_t stride) noexcept -> void { +inline auto transpose_square( + T* RESTRICT matrix, + std::size_t order, + std::size_t stride +) noexcept -> void { const auto empty_matrix{order == 0uz}; if (empty_matrix) { @@ -317,9 +457,16 @@ inline auto transpose_square(T *RESTRICT matrix, std::size_t order, static_cast(xpu::detail::checked_bytes(matrix_size)); #if defined(XPU_CUDA) - const auto range = xpu::range<2uz>{{0uz, 0uz}, {order, order}, {1uz, 1uz}}; + const auto range = xpu::range<2uz>{ + {0uz, 0uz}, + {order, order}, + {1uz, 1uz} + }; - const auto transpose = detail::transpose_square{matrix, stride}; + const auto transpose = detail::transpose_square{ + matrix, + stride + }; xpu::parallel_for(range, transpose); #else @@ -328,17 +475,21 @@ inline auto transpose_square(T *RESTRICT matrix, std::size_t order, const auto col_idx{row * stride + column}; const auto row_idx{column * stride + row}; - std::swap(matrix[col_idx], matrix[row_idx]); + std::swap( + matrix[col_idx], + matrix[row_idx] + ); } } #endif } -template class cholesky_factorization { +template +class cholesky_factorization { private: std::size_t order_; std::size_t stride_; - T *factor_{}; + T* factor_{}; #if defined(XPU_CUDA) using dimension_type = int; @@ -353,8 +504,7 @@ template class cholesky_factorization { #endif [[nodiscard]] - static auto checked_order(std::size_t order, std::size_t stride) noexcept - -> std::size_t { + static auto checked_order(std::size_t order, std::size_t stride) noexcept -> std::size_t { static_cast(xpu::detail::checked_cast(order)); static_cast(xpu::detail::checked_cast(stride)); @@ -373,15 +523,14 @@ template class cholesky_factorization { public: cholesky_factorization(std::size_t order, std::size_t stride) - : order_{checked_order(order, stride)}, stride_{stride} + : order_{checked_order(order, stride)} + , stride_{stride} #if defined(XPU_CUDA) - , - handle_{detail::create_cusolver_handle()}, - workspace_{detail::potrf_workspace_size(handle_, order, stride)}, - info_{1uz} + , handle_{detail::create_cusolver_handle()} + , workspace_{detail::potrf_workspace_size(handle_, order, stride)} + , info_{1uz} #endif - { - } + { } ~cholesky_factorization() { #if defined(XPU_CUDA) @@ -400,43 +549,47 @@ template class cholesky_factorization { } [[nodiscard]] - auto factorize(T *RESTRICT matrix) noexcept -> status { + auto factorize(T* RESTRICT matrix) noexcept -> status { factor_ = nullptr; #if defined(XPU_CUDA) - detail::cusolver_potrf(handle_, order_, stride_, matrix, workspace_.data(), - workspace_.size(), info_.data()); + detail::cusolver_potrf( + handle_, order_, stride_, + matrix, workspace_.data(), workspace_.size(), info_.data() + ); auto info{0}; xpu::copy_n(&info, info_.data(), 1uz); if (info < 0) { detail::linalg_error("cuSOLVER potrf received an invalid argument"); } - if (info > 0) { - return status::not_pd; - } + if (info > 0) { return status::not_pd; } #else - const auto info{detail::lapacke_potrf(matrix, order_, stride_)}; + const auto info{ + detail::lapacke_potrf(matrix, order_, stride_) + }; if (info < 0) { detail::linalg_error("LAPACKE potrf failed"); } - if (info > 0) { - return status::not_pd; - } + if (info > 0) { return status::not_pd; } #endif factor_ = matrix; return status::success; } - auto solve(const T *RESTRICT factor, const T *RESTRICT rhs, - T *RESTRICT solution) noexcept -> void { + auto solve( + const T* RESTRICT factor, + const T* RESTRICT rhs, + T* RESTRICT solution + ) noexcept -> void { if (!factor_) { detail::linalg_error("you must factor before you solve"); } if (factor != factor_) { detail::linalg_error( - "solve requires the most recently factorized matrix"); + "solve requires the most recently factorized matrix" + ); } if (rhs == solution) { detail::linalg_error("rhs and solution must not alias"); @@ -445,8 +598,12 @@ template class cholesky_factorization { xpu::copy_n(solution, rhs, order_); #if defined(XPU_CUDA) - detail::cusolver_potrs(handle_, order_, 1uz, factor, stride_, solution, - order_, info_.data()); + detail::cusolver_potrs( + handle_, order_, 1uz, + factor, stride_, + solution, order_, + info_.data() + ); auto info{0}; xpu::copy_n(&info, info_.data(), 1uz); @@ -454,26 +611,27 @@ template class cholesky_factorization { detail::linalg_error("cuSOLVER potrs received an invalid argument"); } #else - const auto info{detail::lapacke_potrs(factor, solution, order_, stride_)}; + const auto info{ + detail::lapacke_potrs(factor, solution, order_, stride_) + }; if (info != 0) { detail::linalg_error("LAPACKE potrs failed"); } #endif } - auto operator=(const cholesky_factorization &) - -> cholesky_factorization & = delete; - cholesky_factorization(const cholesky_factorization &) = delete; - auto operator=(cholesky_factorization &&) - -> cholesky_factorization & = delete; - cholesky_factorization(cholesky_factorization &&) = delete; + auto operator=(const cholesky_factorization&) -> cholesky_factorization& = delete; + cholesky_factorization(const cholesky_factorization&) = delete; + auto operator=(cholesky_factorization&&) -> cholesky_factorization& = delete; + cholesky_factorization(cholesky_factorization&&) = delete; }; -template class lu_factorization { +template +class lu_factorization { private: std::size_t order_; std::size_t stride_; - T *lower_upper_{}; + T* lower_upper_{}; #if defined(XPU_CUDA) using dimension_type = int; @@ -491,8 +649,7 @@ template class lu_factorization { #endif [[nodiscard]] - static auto checked_order(std::size_t order, std::size_t stride) noexcept - -> std::size_t { + static auto checked_order(std::size_t order, std::size_t stride) noexcept -> std::size_t { static_cast(xpu::detail::checked_cast(order)); static_cast(xpu::detail::checked_cast(stride)); @@ -511,18 +668,17 @@ template class lu_factorization { public: lu_factorization(std::size_t order, std::size_t stride) - : order_{checked_order(order, stride)}, stride_{stride} + : order_{checked_order(order, stride)} + , stride_{stride} #if defined(XPU_CUDA) - , - pivot_{order}, handle_{detail::create_cusolver_handle()}, - workspace_{detail::getrf_workspace_size(handle_, order, stride)}, - info_{1uz} + , pivot_{order} + , handle_{detail::create_cusolver_handle()} + , workspace_{detail::getrf_workspace_size(handle_, order, stride)} + , info_{1uz} #else - , - pivot_{order} + , pivot_{order} #endif - { - } + { } ~lu_factorization() { #if defined(XPU_CUDA) @@ -541,41 +697,46 @@ template class lu_factorization { } [[nodiscard]] - auto factorize(T *RESTRICT matrix) noexcept -> status { + auto factorize(T* RESTRICT matrix) noexcept -> status { lower_upper_ = nullptr; #if defined(XPU_CUDA) - detail::cusolver_getrf(handle_, order_, stride_, matrix, workspace_.data(), - pivot_.data(), info_.data()); + detail::cusolver_getrf( + handle_, order_, stride_, + matrix, workspace_.data(), pivot_.data(), info_.data() + ); auto info{0}; xpu::copy_n(&info, info_.data(), 1uz); if (info < 0) { detail::linalg_error("cuSOLVER getrf received an invalid argument"); } - if (info > 0) { - return status::singular; - } + if (info > 0) { return status::singular; } #else const auto info{ - detail::lapacke_getrf(matrix, pivot_.data(), order_, stride_)}; + detail::lapacke_getrf( + matrix, pivot_.data(), order_, stride_ + ) + }; if (info < 0) { detail::linalg_error("LAPACKE getrf failed"); } - if (info > 0) { - return status::singular; - } + if (info > 0) { return status::singular; } #endif lower_upper_ = matrix; return status::success; } - auto solve(const T *RESTRICT lower_upper, const T *RESTRICT rhs, - T *RESTRICT solution) noexcept -> void { + auto solve( + const T* RESTRICT lower_upper, + const T* RESTRICT rhs, + T* RESTRICT solution + ) noexcept -> void { if (lower_upper != lower_upper_) { detail::linalg_error( - "solve requires the most recently factorized matrix"); + "solve requires the most recently factorized matrix" + ); } if (rhs == solution) { detail::linalg_error("rhs and solution must not alias"); @@ -584,9 +745,14 @@ template class lu_factorization { xpu::copy_n(solution, rhs, order_); #if defined(XPU_CUDA) - detail::cusolver_getrs(handle_, CUBLAS_OP_T, order_, 1uz, lower_upper, - stride_, pivot_.data(), solution, order_, - info_.data()); + detail::cusolver_getrs( + handle_, CUBLAS_OP_T, + order_, 1uz, + lower_upper, stride_, + pivot_.data(), + solution, order_, + info_.data() + ); auto info{0}; xpu::copy_n(&info, info_.data(), 1uz); @@ -594,35 +760,53 @@ template class lu_factorization { detail::linalg_error("cuSOLVER getrs received an invalid argument"); } #else - const auto info{detail::lapacke_getrs(lower_upper, pivot_.data(), solution, - order_, stride_)}; + const auto info{ + detail::lapacke_getrs( + lower_upper, pivot_.data(), solution, + order_, stride_ + ) + }; if (info != 0) { detail::linalg_error("LAPACKE getrs failed"); } #endif } - auto invert(const T *RESTRICT lower_upper, T *RESTRICT inverse) noexcept - -> void { + auto invert( + const T* RESTRICT lower_upper, + T* RESTRICT inverse + ) noexcept -> void { if (lower_upper != lower_upper_) { detail::linalg_error( - "invert requires the most recently factorized matrix"); + "invert requires the most recently factorized matrix" + ); } if (lower_upper == inverse) { detail::linalg_error("lower_upper and inverse must not alias"); } #if defined(XPU_CUDA) - const auto range = - xpu::range<2uz>{{0uz, 0uz}, {order_, order_}, {1uz, 1uz}}; + const auto range = xpu::range<2uz>{ + {0uz, 0uz}, + {order_, order_}, + {1uz, 1uz} + }; - const auto initialize = detail::build_identity{inverse, stride_}; + const auto initialize = detail::build_identity{ + inverse, + stride_ + }; xpu::parallel_for(range, initialize); - detail::cusolver_getrs(handle_, CUBLAS_OP_N, order_, order_, lower_upper, - stride_, pivot_.data(), inverse, stride_, - info_.data()); + detail::cusolver_getrs( + handle_, CUBLAS_OP_N, + order_, order_, + lower_upper, stride_, + pivot_.data(), + inverse, stride_, + info_.data() + ); auto info{0}; xpu::copy_n(&info, info_.data(), 1uz); @@ -631,23 +815,30 @@ template class lu_factorization { } #else for (auto row{0uz}; row < order_; ++row) { - xpu::copy_n(inverse + row * stride_, lower_upper + row * stride_, order_); + xpu::copy_n( + inverse + row * stride_, + lower_upper + row * stride_, + order_ + ); } const auto info{ - detail::lapacke_getri(inverse, pivot_.data(), order_, stride_)}; + detail::lapacke_getri( + inverse, pivot_.data(), order_, stride_ + ) + }; if (info != 0) { detail::linalg_error("LAPACKE getri failed"); } #endif } - auto operator=(const lu_factorization &) -> lu_factorization & = delete; - lu_factorization(const lu_factorization &) = delete; - auto operator=(lu_factorization &&) -> lu_factorization & = delete; - lu_factorization(lu_factorization &&) = delete; + auto operator=(const lu_factorization&) -> lu_factorization& = delete; + lu_factorization(const lu_factorization&) = delete; + auto operator=(lu_factorization&&) -> lu_factorization& = delete; + lu_factorization(lu_factorization&&) = delete; }; -} // namespace linalg +} // namespace xpu::linalg } // namespace xpu From f4d40e8f23ae1ac858125366c123c147215734e6 Mon Sep 17 00:00:00 2001 From: alonso130r <131227801+alonso130r@users.noreply.github.com> Date: Wed, 9 Sep 2026 15:43:57 -0400 Subject: [PATCH 3/3] added tests so i don't get in trouble --- tests/linear_algebra/cases.hpp | 124 ++++++++++++++++++++++++++++++++- 1 file changed, 123 insertions(+), 1 deletion(-) diff --git a/tests/linear_algebra/cases.hpp b/tests/linear_algebra/cases.hpp index 70c756e..4934f3b 100644 --- a/tests/linear_algebra/cases.hpp +++ b/tests/linear_algebra/cases.hpp @@ -157,8 +157,124 @@ int run_linear_algebra_type() { return 0; } +template +int run_cholesky_type() { + constexpr auto order{3uz}; + constexpr auto stride{5uz}; + constexpr auto size{order * stride}; + constexpr auto padding{T{19}}; + constexpr auto tolerance{T{100} * std::numeric_limits::epsilon()}; + + constexpr auto original = std::array{ + T{4}, T{2}, T{2}, + T{2}, T{5}, T{1}, + T{2}, T{1}, T{6} + }; + + auto matrix_host = std::array{}; + matrix_host.fill(padding); + for (auto row{0uz}; row < order; ++row) { + for (auto column{0uz}; column < order; ++column) { + matrix_host[row * stride + column] = original[row * order + column]; + } + } + + auto matrix{xpu::buffer{size}}; + xpu::copy_n(matrix.data(), matrix_host.data(), size); + + auto factorization{xpu::linalg::cholesky_factorization{order, stride}}; + if ( + factorization.order() != order || + factorization.stride() != stride + ) { + return test::fail("cholesky_factorization shape accessors returned the wrong values"); + } + + if ( + factorization.factorize(matrix.data()) != + xpu::linalg::status::success + ) { + return test::fail("cholesky_factorization reported a positive-definite matrix as not positive definite"); + } + + xpu::copy_n(matrix_host.data(), matrix.data(), size); + auto diagonal_product{T{1}}; + for (auto i{0uz}; i < order; ++i) { + diagonal_product *= matrix_host[i * stride + i]; + for (auto column{order}; column < stride; ++column) { + if (matrix_host[i * stride + column] != padding) { + return test::fail("Cholesky factorize modified row padding"); + } + } + } + + if (!test::near(diagonal_product * diagonal_product, T{80}, tolerance)) { + return test::fail("Cholesky diagonal has the wrong determinant"); + } + + constexpr auto rhs_host = std::array{T{14}, T{15}, T{22}}; + constexpr auto expected_solution = std::array{T{1}, T{2}, T{3}}; + auto solution_host = std::array{}; + auto rhs{xpu::buffer{order}}; + auto solution{xpu::buffer{order}}; + xpu::copy_n(rhs.data(), rhs_host.data(), order); + + factorization.solve(matrix.data(), rhs.data(), solution.data()); + xpu::copy_n(solution_host.data(), solution.data(), order); + + for (auto i{0uz}; i < order; ++i) { + if (!test::near(solution_host[i], expected_solution[i], tolerance)) { + return test::fail("Cholesky solve produced the wrong solution"); + } + } + + matrix_host.fill(padding); + constexpr auto not_positive_definite = std::array{ + T{1}, T{0}, T{0}, + T{0}, T{1}, T{0}, + T{0}, T{0}, T{-1} + }; + for (auto row{0uz}; row < order; ++row) { + for (auto column{0uz}; column < order; ++column) { + matrix_host[row * stride + column] = + not_positive_definite[row * order + column]; + } + } + xpu::copy_n(matrix.data(), matrix_host.data(), size); + + if ( + factorization.factorize(matrix.data()) != + xpu::linalg::status::not_pd + ) { + return test::fail("cholesky_factorization did not report a non-positive-definite matrix"); + } + + return 0; +} + inline auto run_linear_algebra_checked_cases() -> int { #if defined(__unix__) + if (const auto failure{test::check_abort([] { + auto factorization{xpu::linalg::cholesky_factorization{std::numeric_limits::max(), 1uz}}; + })}; failure) { + + return failure; + } + + if (const auto failure{test::check_abort([] { + auto factorization{xpu::linalg::cholesky_factorization{1uz, std::numeric_limits::max()}}; + })}; failure) { + + return failure; + } + + if (const auto failure{test::check_abort([] { + auto factorization{xpu::linalg::cholesky_factorization{2uz, 1uz}}; + })}; failure) { + + return failure; + } + if (const auto failure{test::check_abort([] { auto factorization{xpu::linalg::lu_factorization{std::numeric_limits::max(), 1uz}}; })}; failure) { @@ -217,5 +333,11 @@ inline int run_linear_algebra_cases() { if (const auto status{run_linear_algebra_type()}; status != 0) { return status; } - return run_linear_algebra_type(); + if (const auto status{run_linear_algebra_type()}; status != 0) { + return status; + } + if (const auto status{run_cholesky_type()}; status != 0) { + return status; + } + return run_cholesky_type(); }