Skip to content

Use experimental complex extension for all complex arithmetic #2069

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 9 commits into
base: master
Choose a base branch
from
Draft

Use experimental complex extension for all complex arithmetic #2069

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
ae75273
Use sycl complex extension throughout element-wise and utils
ndgrigorian Apr 26, 2025
88e64c1
Update binary functions multiply and subtract to use experimental SYC…
ndgrigorian Apr 29, 2025
01f22d3
Use experimental SYCL complex in dot product
ndgrigorian Apr 29, 2025
34168bd
Use experimental namespace in sequential dot product
ndgrigorian Apr 30, 2025
0fd49ae
Use experimental complex namespace in gemm
ndgrigorian Apr 30, 2025
e3be74e
Use specialized functor for multiplying or adding complex inputs
ndgrigorian Apr 30, 2025
6fba5f2
Use custom plus operator in gemm and dot product tree reduction kernels
ndgrigorian Apr 30, 2025
da15d0e
Refactor operators when dispatching to tree reductions
ndgrigorian Apr 30, 2025
52bb73e
Refactor sycl_complex indirect include
ndgrigorian Apr 30, 2025
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
21 changes: 11 additions & 10 deletions dpctl/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
#include <sycl/sycl.hpp>
#include <type_traits>

#include "sycl_complex.hpp"
#include "utils/sycl_complex.hpp"
#include "vec_size_util.hpp"

#include "kernels/dpctl_tensor_types.hpp"
Expand All @@ -49,7 +49,9 @@ namespace acos
{

using dpctl::tensor::ssize_t;
namespace su_ns = dpctl::tensor::sycl_utils;
namespace td_ns = dpctl::tensor::type_dispatch;
namespace exprm_ns = sycl::ext::oneapi::experimental;

using dpctl::tensor::type_utils::is_complex;

Expand All @@ -72,9 +74,10 @@ template <typename argT, typename resT> struct AcosFunctor
using realT = typename argT::value_type;

constexpr realT q_nan = std::numeric_limits<realT>::quiet_NaN();

const realT x = std::real(in);
const realT y = std::imag(in);
using sycl_complexT = su_ns::sycl_complex_t<realT>;
sycl_complexT z = sycl_complexT(in);
const realT x = exprm_ns::real(z);
const realT y = exprm_ns::imag(z);

if (std::isnan(x)) {
/* acos(NaN + I*+-Inf) = NaN + I*-+Inf */
Expand Down Expand Up @@ -106,20 +109,18 @@ template <typename argT, typename resT> struct AcosFunctor
constexpr realT r_eps =
realT(1) / std::numeric_limits<realT>::epsilon();
if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) {
using sycl_complexT = exprm_ns::complex<realT>;
sycl_complexT log_in =
exprm_ns::log(exprm_ns::complex<realT>(in));
sycl_complexT log_z = exprm_ns::log(z);

const realT wx = log_in.real();
const realT wy = log_in.imag();
const realT wx = log_z.real();
const realT wy = log_z.imag();
const realT rx = sycl::fabs(wy);

realT ry = wx + sycl::log(realT(2));
return resT{rx, (sycl::signbit(y)) ? ry : -ry};
}

/* ordinary cases */
return exprm_ns::acos(exprm_ns::complex<realT>(in)); // acos(in);
return exprm_ns::acos(z); // acos(z);
}
else {
static_assert(std::is_floating_point_v<argT> ||
Expand Down
40 changes: 21 additions & 19 deletions dpctl/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
#include <sycl/sycl.hpp>
#include <type_traits>

#include "sycl_complex.hpp"
#include "utils/sycl_complex.hpp"
#include "vec_size_util.hpp"

#include "kernels/dpctl_tensor_types.hpp"
Expand All @@ -49,7 +49,9 @@ namespace acosh
{

using dpctl::tensor::ssize_t;
namespace su_ns = dpctl::tensor::sycl_utils;
namespace td_ns = dpctl::tensor::type_dispatch;
namespace exprm_ns = sycl::ext::oneapi::experimental;

using dpctl::tensor::type_utils::is_complex;

Expand Down Expand Up @@ -77,33 +79,35 @@ template <typename argT, typename resT> struct AcoshFunctor
* where the sign is chosen so Re(acosh(in)) >= 0.
* So, we first calculate acos(in) and then acosh(in).
*/
const realT x = std::real(in);
const realT y = std::imag(in);
using sycl_complexT = su_ns::sycl_complex_t<realT>;
sycl_complexT z = sycl_complexT(in);
const realT x = exprm_ns::real(z);
const realT y = exprm_ns::imag(z);

resT acos_in;
sycl_complexT acos_z;
if (std::isnan(x)) {
/* acos(NaN + I*+-Inf) = NaN + I*-+Inf */
if (std::isinf(y)) {
acos_in = resT{q_nan, -y};
acos_z = resT{q_nan, -y};
}
else {
acos_in = resT{q_nan, q_nan};
acos_z = resT{q_nan, q_nan};
}
}
else if (std::isnan(y)) {
/* acos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
constexpr realT inf = std::numeric_limits<realT>::infinity();

if (std::isinf(x)) {
acos_in = resT{q_nan, -inf};
acos_z = resT{q_nan, -inf};
}
/* acos(0 + I*NaN) = Pi/2 + I*NaN with inexact */
else if (x == realT(0)) {
const realT pi_half = sycl::atan(realT(1)) * 2;
acos_in = resT{pi_half, q_nan};
acos_z = resT{pi_half, q_nan};
}
else {
acos_in = resT{q_nan, q_nan};
acos_z = resT{q_nan, q_nan};
}
}

Expand All @@ -113,23 +117,21 @@ template <typename argT, typename resT> struct AcoshFunctor
* For large x or y including acos(+-Inf + I*+-Inf)
*/
if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) {
using sycl_complexT = typename exprm_ns::complex<realT>;
const sycl_complexT log_in = exprm_ns::log(sycl_complexT(in));
const realT wx = log_in.real();
const realT wy = log_in.imag();
const sycl_complexT log_z = exprm_ns::log(z);
const realT wx = log_z.real();
const realT wy = log_z.imag();
const realT rx = sycl::fabs(wy);
realT ry = wx + sycl::log(realT(2));
acos_in = resT{rx, (sycl::signbit(y)) ? ry : -ry};
acos_z = resT{rx, (sycl::signbit(y)) ? ry : -ry};
}
else {
/* ordinary cases */
acos_in =
exprm_ns::acos(exprm_ns::complex<realT>(in)); // acos(in);
acos_z = exprm_ns::acos(z); // acos(z);
}

/* Now we calculate acosh(z) */
const realT rx = std::real(acos_in);
const realT ry = std::imag(acos_in);
const realT rx = exprm_ns::real(acos_z);
const realT ry = exprm_ns::imag(acos_z);

/* acosh(NaN + I*NaN) = NaN + I*NaN */
if (std::isnan(rx) && std::isnan(ry)) {
Expand All @@ -145,7 +147,7 @@ template <typename argT, typename resT> struct AcoshFunctor
return resT{ry, ry};
}
/* ordinary cases */
const realT res_im = sycl::copysign(rx, std::imag(in));
const realT res_im = sycl::copysign(rx, exprm_ns::imag(z));
return resT{sycl::fabs(ry), res_im};
}
else {
Expand Down
27 changes: 22 additions & 5 deletions dpctl/tensor/libtensor/include/kernels/elementwise_functions/add.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
#include <sycl/sycl.hpp>
#include <type_traits>

#include "sycl_complex.hpp"
#include "utils/sycl_complex.hpp"
#include "vec_size_util.hpp"

#include "utils/offset_utils.hpp"
Expand All @@ -50,8 +50,10 @@ namespace add
{

using dpctl::tensor::ssize_t;
namespace su_ns = dpctl::tensor::sycl_utils;
namespace td_ns = dpctl::tensor::type_dispatch;
namespace tu_ns = dpctl::tensor::type_utils;
namespace exprm_ns = sycl::ext::oneapi::experimental;

template <typename argT1, typename argT2, typename resT> struct AddFunctor
{
Expand All @@ -69,21 +71,22 @@ template <typename argT1, typename argT2, typename resT> struct AddFunctor
using rT1 = typename argT1::value_type;
using rT2 = typename argT2::value_type;

return exprm_ns::complex<rT1>(in1) + exprm_ns::complex<rT2>(in2);
return su_ns::sycl_complex_t<rT1>(in1) +
su_ns::sycl_complex_t<rT2>(in2);
}
else if constexpr (tu_ns::is_complex<argT1>::value &&
!tu_ns::is_complex<argT2>::value)
{
using rT1 = typename argT1::value_type;

return exprm_ns::complex<rT1>(in1) + in2;
return su_ns::sycl_complex_t<rT1>(in1) + in2;
}
else if constexpr (!tu_ns::is_complex<argT1>::value &&
tu_ns::is_complex<argT2>::value)
{
using rT2 = typename argT2::value_type;

return in1 + exprm_ns::complex<rT2>(in2);
return in1 + su_ns::sycl_complex_t<rT2>(in2);
}
else {
return in1 + in2;
Expand Down Expand Up @@ -460,7 +463,21 @@ template <typename argT, typename resT> struct AddInplaceFunctor
using supports_vec = std::negation<
std::disjunction<tu_ns::is_complex<argT>, tu_ns::is_complex<resT>>>;

void operator()(resT &res, const argT &in) { res += in; }
void operator()(resT &res, const argT &in)
{
if constexpr (tu_ns::is_complex_v<resT> && tu_ns::is_complex_v<argT>) {
using rT1 = typename resT::value_type;
using rT2 = typename argT::value_type;

auto tmp = su_ns::sycl_complex_t<rT1>(res);
tmp += su_ns::sycl_complex_t<rT2>(in);

res = resT(tmp);
}
else {
res += in;
}
}

template <int vec_sz>
void operator()(sycl::vec<resT, vec_sz> &res,
Expand Down
6 changes: 4 additions & 2 deletions dpctl/tensor/libtensor/include/kernels/elementwise_functions/angle.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@
#include <sycl/sycl.hpp>
#include <type_traits>

#include "sycl_complex.hpp"
#include "utils/sycl_complex.hpp"
#include "vec_size_util.hpp"

#include "kernels/dpctl_tensor_types.hpp"
Expand All @@ -50,7 +50,9 @@ namespace angle
{

using dpctl::tensor::ssize_t;
namespace su_ns = dpctl::tensor::sycl_utils;
namespace td_ns = dpctl::tensor::type_dispatch;
namespace exprm_ns = sycl::ext::oneapi::experimental;

using dpctl::tensor::type_utils::is_complex;

Expand All @@ -71,7 +73,7 @@ template <typename argT, typename resT> struct AngleFunctor
{
using rT = typename argT::value_type;

return exprm_ns::arg(exprm_ns::complex<rT>(in)); // arg(in);
return exprm_ns::arg(su_ns::sycl_complex_t<rT>(in)); // arg(in);
}
};

Expand Down
28 changes: 15 additions & 13 deletions dpctl/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
#include <sycl/sycl.hpp>
#include <type_traits>

#include "sycl_complex.hpp"
#include "utils/sycl_complex.hpp"
#include "vec_size_util.hpp"

#include "kernels/dpctl_tensor_types.hpp"
Expand All @@ -49,7 +49,9 @@ namespace asin
{

using dpctl::tensor::ssize_t;
namespace su_ns = dpctl::tensor::sycl_utils;
namespace td_ns = dpctl::tensor::type_dispatch;
namespace exprm_ns = sycl::ext::oneapi::experimental;

using dpctl::tensor::type_utils::is_complex;

Expand Down Expand Up @@ -80,8 +82,10 @@ template <typename argT, typename resT> struct AsinFunctor
* y = imag(I * conj(in)) = real(in)
* and then return {imag(w), real(w)} which is asin(in)
*/
const realT x = std::imag(in);
const realT y = std::real(in);
using sycl_complexT = su_ns::sycl_complex_t<realT>;
sycl_complexT z = sycl_complexT(in);
const realT x = exprm_ns::imag(z);
const realT y = exprm_ns::real(z);

if (std::isnan(x)) {
/* asinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
Expand Down Expand Up @@ -120,26 +124,24 @@ template <typename argT, typename resT> struct AsinFunctor
constexpr realT r_eps =
realT(1) / std::numeric_limits<realT>::epsilon();
if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) {
using sycl_complexT = exprm_ns::complex<realT>;
const sycl_complexT z{x, y};
const sycl_complexT z1{x, y};
realT wx, wy;
if (!sycl::signbit(x)) {
const auto log_z = exprm_ns::log(z);
wx = log_z.real() + sycl::log(realT(2));
wy = log_z.imag();
const auto log_z1 = exprm_ns::log(z1);
wx = log_z1.real() + sycl::log(realT(2));
wy = log_z1.imag();
}
else {
const auto log_mz = exprm_ns::log(-z);
wx = log_mz.real() + sycl::log(realT(2));
wy = log_mz.imag();
const auto log_mz1 = exprm_ns::log(-z1);
wx = log_mz1.real() + sycl::log(realT(2));
wy = log_mz1.imag();
}
const realT asinh_re = sycl::copysign(wx, x);
const realT asinh_im = sycl::copysign(wy, y);
return resT{asinh_im, asinh_re};
}
/* ordinary cases */
return exprm_ns::asin(
exprm_ns::complex<realT>(in)); // sycl::asin(in);
return exprm_ns::asin(z); // sycl::asin(z);
}
else {
static_assert(std::is_floating_point_v<argT> ||
Expand Down
23 changes: 12 additions & 11 deletions dpctl/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@
#include <sycl/sycl.hpp>
#include <type_traits>

#include "sycl_complex.hpp"
#include "utils/sycl_complex.hpp"
#include "vec_size_util.hpp"

#include "kernels/dpctl_tensor_types.hpp"
Expand All @@ -49,7 +49,9 @@ namespace asinh
{

using dpctl::tensor::ssize_t;
namespace su_ns = dpctl::tensor::sycl_utils;
namespace td_ns = dpctl::tensor::type_dispatch;
namespace exprm_ns = sycl::ext::oneapi::experimental;

using dpctl::tensor::type_utils::is_complex;

Expand All @@ -72,9 +74,10 @@ template <typename argT, typename resT> struct AsinhFunctor
using realT = typename argT::value_type;

constexpr realT q_nan = std::numeric_limits<realT>::quiet_NaN();

const realT x = std::real(in);
const realT y = std::imag(in);
using sycl_complexT = su_ns::sycl_complex_t<realT>;
sycl_complexT z = sycl_complexT(in);
const realT x = exprm_ns::real(z);
const realT y = exprm_ns::imag(z);

if (std::isnan(x)) {
/* asinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
Expand Down Expand Up @@ -109,20 +112,18 @@ template <typename argT, typename resT> struct AsinhFunctor
realT(1) / std::numeric_limits<realT>::epsilon();

if (sycl::fabs(x) > r_eps || sycl::fabs(y) > r_eps) {
using sycl_complexT = exprm_ns::complex<realT>;
sycl_complexT log_in = (sycl::signbit(x))
? exprm_ns::log(sycl_complexT(-in))
: exprm_ns::log(sycl_complexT(in));
realT wx = log_in.real() + sycl::log(realT(2));
realT wy = log_in.imag();
sycl_complexT log_in =
(sycl::signbit(x)) ? exprm_ns::log(-z) : exprm_ns::log(z);
realT wx = exprm_ns::real(log_in) + sycl::log(realT(2));
realT wy = exprm_ns::imag(log_in);

const realT res_re = sycl::copysign(wx, x);
const realT res_im = sycl::copysign(wy, y);
return resT{res_re, res_im};
}

/* ordinary cases */
return exprm_ns::asinh(exprm_ns::complex<realT>(in)); // asinh(in);
return exprm_ns::asinh(z); // asinh(z);
}
else {
static_assert(std::is_floating_point_v<argT> ||
Expand Down
Loading
Loading