Casting from cuDoubleComplex to std::complex<double>

[Danesh Daroui]

Hi,

In our code we use std::complex<complex> in all places. Now if we want to use cuDoubleComplex which is a typedef of double2 and used by MAGMA for complex operations with double precision, we need to do modifications in many places. MKL routines allow you to replace their MKL_COMPLEX16 with std::complex<double> with #define compiler directive. Is there similar trick to force MAGMA to use std::complex<double> instead?

Regards,

Dan

[Natalie Beams]

Hi Dan, 

Theoretically, you can change the definition of magmaDoubleComplex in magma_types.h (https://bitbucket.org/icl/magma/src/3bec2af7a08f32d78e413bf4203a64532764460e/include/magma_types.h#lines-97)

However, this will cause several issues with MAGMA's math macros, and generate tons of compiler warnings since MAGMA's interface is C and 

many routines will now include non-C-friendly types.

In the initial stages of the ongoing SYCL/Intel GPU port, we used std::complex because that is what oneMKL uses. However, we recently changed

it to use a C-friendly struct, instead. If you look at the commit that made this change, you can see some of the things we had to change when using

std::complex in the removed code: e.g., the macros at https://bitbucket.org/icl/magma/commits/6e66cc79b4979be8875c20131be1feeafd3c15de#Linclude/magma_types.hF182

and https://bitbucket.org/icl/magma/commits/6e66cc79b4979be8875c20131be1feeafd3c15de#Linclude/magma_types.hF199

Another potential issue is the operators defined in include/magma_operators.h. When using std::complex, the compiler may not know how to resolve
certain overloaded operators (doesn't know whether to use the one from std::complex or from MAGMA). For the SYCL development, this popped up 

with the `*=` operator in zlascl, for example. This occurs when the std::complex function is not a template function; then the compiler doesn't know

whether to prefer MAGMA's non-template function or std::complex.

Finally, a big issue with changing MAGMA's complex type is all the places that it interacts with cuBLAS, which would require casting. In the updated SYCL

version, we define some macros to help with the casting required to go from the C type to std::complex for oneMKL: 
https://bitbucket.org/icl/magma/commits/6e66cc79b4979be8875c20131be1feeafd3c15de#Linclude/magma_types.hT182. But in this case it would 
need to go the other way, because you'd need to cast from std::complex to cuDoubleComplex.

So to sum up: you could make it work, but it will require quite a bit of modification to MAGMA. If you want to go this route, I am happy to help answer

any more questions you run into based on my experience working with the SYCL branch.  But depending on your code, it may still be a lot simpler to
build MAGMA with cuDoubleComplex, then do your own casting to magmaDoubleComplex/cuDoubleComplex when calling MAGMA.
As MAGMA has always had a C interface, we don't currently have plans to officially support std::complex.

-- Natalie

[Danesh Daroui]

Hi Natalie,

Many thanks for your response. I solved the problem by casting the pointer using reinterpret_cast and it works fine:

std::complex<double>* A;

vector<complex<double> > x(S);

// Memory allocation for A is done here and the coefficient matrix is filled.

cuDoubleComplex* dA;
cuDoubleComplex* dx;

dA = reinterpret_cast<cuDoubleComplex*>(A);
dx = reinterpret_cast<cuDoubleComplex*>(&x[0]);

Hopefully I am doing it correctly. The results from MAGMA is exactly same as MKL so, this seems to work. I think oneMKL uses MKL_Complex16 instead of std::complex<double> but that can easily be changed using compiler directives. MKL_Complex16 is similar to double2 that MAGMA uses since this is a C struct with two parameters for real and imaginary values.

Regards,

Dan

[Natalie Beams]

Hi Dan,

That's great that it will work for you to cast to cuDoubleComplex prior to calling MAGMA. I think that will be the easiest option by far. I know CUDA guarantees

that you can cast between cuDoubleComplex and std::complex<double>, so it should be safe as far as I can tell. 

I apologize if I introduced any confusion with my mention of oneMKL. For a long time I used that to refer to the SYCL version of MKL, vs the standard C/Fortran

MKL interface, but I recently learned (and keep forgetting) that I should specify "the DPC++ interface of oneMKL" (https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-dpcpp/2024-2/overview.html) as the "standard" interface is also part of oneMKL. The DPC++ interface is C++ only and uses std::complex for complex types everywhere.

-- Natalie