Generalize the classes to use complex or real numbers

[PabloIbannez]

No description provided.

[RaulPPelaez]

Is leaving this as real a conscious decision?

[PabloPalaciosAlonso]

Yes, we believe that since H is associated with a size, it will always be real. Furthermore, in the complex_bvp branch (which is the branch we're trying to merge into the main one with these changes), H is treated as a real number, while all other entities are complex.

[RaulPPelaez]

At first I though this PR was making some things compatible with either float or double, thus some of my comments -.-
Calling the template parameter U seems to signal that this class can be instantiated with any type (like float, float4, complex...). I guess this is why concepts were introduced in C++20, so you can restrict the type to be either thrust::complex or real and nothing else.
I now remember I did not merge the complex branch because I found no way to state this clearly via code.
The code for the complex and real versions is really similar but not quite so. I considered making everything complex, but many usecases are real and that would entail twice the work.
Maybe putting everything into a hidden namespace and just exposing two instantiated types (real and thrust::complex) for the public API is the best way to go. Then if you want to use another complex type or whatever you know you are on your own.

[RaulPPelaez]

I am really scared of the changes to this file. cufftComplex types are just an alias to float2 or double2. thrust::complex actually has the correct complex algebra.
By defining complex to thrust::complex you changed the behavior of any line that does something like:

complex a,b = ...
complex c = a*b;

[RaulPPelaez]

Could you comment on what you did here?
Did you ran the tests for these components locally?

[RaulPPelaez]

I did something really similar when I needed the BVP solver to run with complex:
v2.x...complex_bvp
But in that case I wrote it to be compatible just with complex and using thrust::complex.

My point is this function here:

    //Given A (2x2 matrix) and b as two numbers of an arbitrary type (could be complex, real...)
    template<class T, class T2>
    __device__ thrust::pair<T,T> solve2x2System(T2 A[4], thrust::pair<T,T> b){
      const auto det  = A[0]*A[3] - A[1]*A[2];
      const T c0 = (b.first*A[3] - A[1]*b.second)/det;
      const T d0 = (b.second*A[0] - b.first*A[2])/det;
      return thrust::make_pair(c0, d0);
    }

Is incorrect if the type passed to it is cufftComplex instead of thrust::complex.

[RaulPPelaez]

For instance, these tests might fail, since they define complex with cufftComplex:
https://github.com/RaulPPelaez/UAMMD/blob/v2.x/test/misc/bvp/bvp.cu

[PabloPalaciosAlonso]

I did something really similar when I needed the BVP solver to run with complex: v2.x...complex_bvp But in that case I wrote it to be compatible just with complex and using thrust::complex.

My point is this function here:

    //Given A (2x2 matrix) and b as two numbers of an arbitrary type (could be complex, real...)
    template<class T, class T2>
    __device__ thrust::pair<T,T> solve2x2System(T2 A[4], thrust::pair<T,T> b){
      const auto det  = A[0]*A[3] - A[1]*A[2];
      const T c0 = (b.first*A[3] - A[1]*b.second)/det;
      const T d0 = (b.second*A[0] - b.first*A[2])/det;
      return thrust::make_pair(c0, d0);
    }

Is incorrect if the type passed to it is cufftComplex instead of thrust::complex.

I see your point, probably this function is a copy-paste of what you did in BVP_complex. One solution could be to specialize this function for the most common complex types?

[RaulPPelaez]

I am not sure that is worth it. There are probably other places in which the code is expecting either cufftComplex or thrust::complex, and the difference there are super tricky and subtle.
For instance this code here will be broken in a very subtle and mostly silent way if you do not use thrust::complex and the waveVector is complex:

UAMMD/src/misc/BoundaryValueProblem/BVPSolver.cuh

Line 183 in b1da83f

const real kH2 = k*k*H*H;

Same here:

UAMMD/src/misc/BoundaryValueProblem/BVPSolver.cuh

Lines 70 to 71 in 73335f7

	auto B00 = -k*k*H*H;
	auto B11 = -k*k*H*H;

I would say these classes should work as is with thrust::complex and something should be done to ensure they cannot be instantiated with cufftComplex.

For starters, the tests in here related to the BVP solver should be included as part of this PR, those tests include complex and are super thorough. We cannot trust these changes unless that passes. Sadly they have to be ran locally since github CI machines do not have GPUs.

Most of these classes are internal implementation details and should typically not be instantiated by users of the library. One solution here could be to make these classes a template:

UAMMD/src/misc/BoundaryValueProblem/BVPSolver.cuh

Line 151 in b1da83f

class BoundaryValueProblemSolver{

UAMMD/src/misc/BoundaryValueProblem/BVPSolver.cuh

Line 228 in b1da83f

class BatchedBVPHandler{

Call them something like name_impl and then generate two aliases like:

using BVPSolverComplex = BVPSolver_impl<thrust::complex<real>>;
using BVPSolverReal = BVPSolver_impl<real>;

This at least signals that you are not supposed to instantiate _impl with just about any type.
Any option is ok as long as the documentation reflects it. You can modify the documentation for these components here:
https://github.com/RaulPPelaez/UAMMD/blob/v2.x/docs/BoundaryValueProblemSolver.rst