Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 34 additions & 34 deletions xpu/linear_algebra.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ struct build_identity {
std::size_t stride;

DEVICE_ONLY
auto operator()(const xpu::array<std::size_t, 2uz>& index) const -> void {
auto operator()(const array<std::size_t, 2uz>& index) const -> void {
const auto row{index[0]};
const auto column{index[1]};

Expand All @@ -66,7 +66,7 @@ struct transpose_square {
std::size_t stride;

DEVICE_ONLY
auto operator()(const xpu::array<std::size_t, 2uz>& index) const -> void {
auto operator()(const array<std::size_t, 2uz>& index) const -> void {
const auto row{index[0]};
const auto column{index[1]};
const auto skip_pair{row >= column};
Expand All @@ -87,7 +87,7 @@ struct transpose_square {

inline auto create_cusolver_handle() -> cusolverDnHandle_t {
auto handle{cusolverDnHandle_t{}};
xpu::cu_check(cusolverDnCreate(&handle));
cu_check(cusolverDnCreate(&handle));
return handle;
}

Expand All @@ -103,14 +103,14 @@ inline auto getrf_workspace_size(
auto size{0};

if constexpr (std::same_as<T, float>) {
xpu::cu_check(cusolverDnSgetrf_bufferSize(
cu_check(cusolverDnSgetrf_bufferSize(
handle,
vendor_order, vendor_order,
nullptr, vendor_stride,
&size
));
} else {
xpu::cu_check(cusolverDnDgetrf_bufferSize(
cu_check(cusolverDnDgetrf_bufferSize(
handle,
vendor_order, vendor_order,
nullptr, vendor_stride,
Expand All @@ -135,14 +135,14 @@ inline auto potrf_workspace_size(
auto size{0};

if constexpr (std::same_as<T, float>) {
xpu::cu_check(cusolverDnSpotrf_bufferSize(
cu_check(cusolverDnSpotrf_bufferSize(
handle, CUBLAS_FILL_MODE_UPPER,
vendor_order,
nullptr, vendor_stride,
&size
));
} else {
xpu::cu_check(cusolverDnDpotrf_bufferSize(
cu_check(cusolverDnDpotrf_bufferSize(
handle, CUBLAS_FILL_MODE_UPPER,
vendor_order,
nullptr, vendor_stride,
Expand All @@ -169,14 +169,14 @@ inline auto cusolver_getrf(
const auto vendor_stride{xpu::detail::checked_cast<int>(stride)};

if constexpr (std::same_as<T, float>) {
xpu::cu_check(cusolverDnSgetrf(
cu_check(cusolverDnSgetrf(
handle,
vendor_order, vendor_order,
matrix, vendor_stride,
workspace, pivot, info
));
} else {
xpu::cu_check(cusolverDnDgetrf(
cu_check(cusolverDnDgetrf(
handle,
vendor_order, vendor_order,
matrix, vendor_stride,
Expand All @@ -200,15 +200,15 @@ inline auto cusolver_potrf(
const auto vendor_workspace_size{xpu::detail::checked_cast<int>(workspace_size)};

if constexpr (std::same_as<T, float>) {
xpu::cu_check(cusolverDnSpotrf(
cu_check(cusolverDnSpotrf(
handle, CUBLAS_FILL_MODE_UPPER,
vendor_order,
matrix, vendor_stride,
workspace, vendor_workspace_size,
info
));
} else {
xpu::cu_check(cusolverDnDpotrf(
cu_check(cusolverDnDpotrf(
handle, CUBLAS_FILL_MODE_UPPER,
vendor_order,
matrix, vendor_stride,
Expand Down Expand Up @@ -237,7 +237,7 @@ inline auto cusolver_getrs(
const auto vendor_solution_stride{xpu::detail::checked_cast<int>(solution_stride)};

if constexpr (std::same_as<T, float>) {
xpu::cu_check(cusolverDnSgetrs(
cu_check(cusolverDnSgetrs(
handle, operation,
vendor_order, vendor_right_hand_sides,
lower_upper, vendor_lower_upper_stride,
Expand All @@ -246,7 +246,7 @@ inline auto cusolver_getrs(
info
));
} else {
xpu::cu_check(cusolverDnDgetrs(
cu_check(cusolverDnDgetrs(
handle, operation,
vendor_order, vendor_right_hand_sides,
lower_upper, vendor_lower_upper_stride,
Expand Down Expand Up @@ -274,15 +274,15 @@ inline auto cusolver_potrs(
const auto vendor_solution_stride{xpu::detail::checked_cast<int>(solution_stride)};

if constexpr (std::same_as<T, float>) {
xpu::cu_check(cusolverDnSpotrs(
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(
cu_check(cusolverDnDpotrs(
handle, CUBLAS_FILL_MODE_UPPER,
vendor_order, vendor_right_hand_sides,
factor, vendor_factor_stride,
Expand Down Expand Up @@ -468,7 +468,7 @@ inline auto transpose_square(
stride
};

xpu::parallel_for(range, transpose);
parallel_for(range, transpose);
#else
for (auto row{0uz}; row < order; ++row) {
for (auto column{row + 1uz}; column < order; ++column) {
Expand Down Expand Up @@ -499,8 +499,8 @@ class cholesky_factorization {

#if defined(XPU_CUDA)
cusolverDnHandle_t handle_;
xpu::buffer<T> workspace_;
xpu::buffer<int> info_;
buffer<T> workspace_;
buffer<int> info_;
#endif

[[nodiscard]]
Expand Down Expand Up @@ -534,7 +534,7 @@ class cholesky_factorization {

~cholesky_factorization() {
#if defined(XPU_CUDA)
xpu::cu_check(cusolverDnDestroy(handle_));
cu_check(cusolverDnDestroy(handle_));
#endif
}

Expand All @@ -559,7 +559,7 @@ class cholesky_factorization {
);

auto info{0};
xpu::copy_n(&info, info_.data(), 1uz);
copy_n(&info, info_.data(), 1uz);
if (info < 0) {
detail::linalg_error("cuSOLVER potrf received an invalid argument");
}
Expand Down Expand Up @@ -595,7 +595,7 @@ class cholesky_factorization {
detail::linalg_error("rhs and solution must not alias");
}

xpu::copy_n(solution, rhs, order_);
copy_n(solution, rhs, order_);

#if defined(XPU_CUDA)
detail::cusolver_potrs(
Expand All @@ -606,7 +606,7 @@ class cholesky_factorization {
);

auto info{0};
xpu::copy_n(&info, info_.data(), 1uz);
copy_n(&info, info_.data(), 1uz);
if (info != 0) {
detail::linalg_error("cuSOLVER potrs received an invalid argument");
}
Expand Down Expand Up @@ -640,12 +640,12 @@ class lu_factorization {
#endif

#if defined(XPU_CUDA)
xpu::buffer<int> pivot_;
buffer<int> pivot_;
cusolverDnHandle_t handle_;
xpu::buffer<T> workspace_;
xpu::buffer<int> info_;
buffer<T> workspace_;
buffer<int> info_;
#else
xpu::buffer<lapack_int> pivot_;
buffer<lapack_int> pivot_;
#endif

[[nodiscard]]
Expand All @@ -661,7 +661,7 @@ class lu_factorization {

const auto matrix_size{xpu::detail::checked_mul(order, stride)};

static_cast<void>(xpu::detail::checked_bytes<T>(matrix_size));
static_cast<void>(bytes<T>(matrix_size));

return order;
}
Expand All @@ -682,7 +682,7 @@ class lu_factorization {

~lu_factorization() {
#if defined(XPU_CUDA)
xpu::cu_check(cusolverDnDestroy(handle_));
cu_check(cusolverDnDestroy(handle_));
#endif
}

Expand All @@ -707,7 +707,7 @@ class lu_factorization {
);

auto info{0};
xpu::copy_n(&info, info_.data(), 1uz);
copy_n(&info, info_.data(), 1uz);
if (info < 0) {
detail::linalg_error("cuSOLVER getrf received an invalid argument");
}
Expand Down Expand Up @@ -742,7 +742,7 @@ class lu_factorization {
detail::linalg_error("rhs and solution must not alias");
}

xpu::copy_n(solution, rhs, order_);
copy_n(solution, rhs, order_);

#if defined(XPU_CUDA)
detail::cusolver_getrs(
Expand All @@ -755,7 +755,7 @@ class lu_factorization {
);

auto info{0};
xpu::copy_n(&info, info_.data(), 1uz);
copy_n(&info, info_.data(), 1uz);
if (info != 0) {
detail::linalg_error("cuSOLVER getrs received an invalid argument");
}
Expand Down Expand Up @@ -797,7 +797,7 @@ class lu_factorization {
stride_
};

xpu::parallel_for(range, initialize);
parallel_for(range, initialize);

detail::cusolver_getrs(
handle_, CUBLAS_OP_N,
Expand All @@ -809,13 +809,13 @@ class lu_factorization {
);

auto info{0};
xpu::copy_n(&info, info_.data(), 1uz);
copy_n(&info, info_.data(), 1uz);
if (info != 0) {
detail::linalg_error("cuSOLVER getrs received an invalid argument");
}
#else
for (auto row{0uz}; row < order_; ++row) {
xpu::copy_n(
copy_n(
inverse + row * stride_,
lower_upper + row * stride_,
order_
Expand Down
13 changes: 7 additions & 6 deletions xpu/math.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -67,25 +67,26 @@ using xstd::ldexp; using xstd::frexp; using xstd::modf;

// Overflow safe version of (num + den - 1) / den
template <std::unsigned_integral T> [[nodiscard]] CUDA_CALLABLE
inline constexpr auto ceiling_div(T num, T den) noexcept -> T {
constexpr auto ceiling_div(T num, T den) noexcept -> T {
const auto zero_denominator{den == 0};

if (zero_denominator) {
xpu::detail::checked_error("zero ceiling division denominator");
}

const auto quotient{num / den + (num % den != 0)};
const auto quotient{ ( num / den ) + (num % den != 0)};

return quotient;
}


template <std::floating_point T> CUDA_CALLABLE
inline auto sincos(T arg, T* RESTRICT s, T* RESTRICT c) noexcept -> void {
inline auto sincos(T arg, T* RESTRICT sin, T* RESTRICT cos) noexcept -> void {
#if defined(__CUDA_ARCH__)
if constexpr (std::is_same_v<T, float>) { ::sincosf(arg, s, c); }
else { ::sincos(arg, s, c); }
if constexpr (std::is_same_v<T, float>) { ::sincosf(arg, sin, cos); }
else { ::sincos(arg, sin, cos); }
#else
*s = xpu::sin(arg); *c = xpu::cos(arg);
*sin = xpu::sin(arg); *cos = xpu::cos(arg);
#endif
}

Expand Down
Loading
Loading