mytypes.hpp

Go to the documentation of this file.

#pragma once
 
#ifdef _OPENMP
    #define EIGEN_DONT_PARALLELIZE
#endif
 
#if __has_include(<mkl_types.h>)
    #ifndef MKL
        #define MKL
    #endif
    #ifndef EIGEN_USE_MKL_ALL
        #define EIGEN_USE_MKL_ALL
    #endif
#else
    #if __has_include(<Accelerate/Accelerate.h>)
        #ifndef ACCELERATE
            #define ACCELERATE
        #endif
    #endif
#endif
 
#define EIGEN_DEFAULT_IO_FORMAT \
    Eigen::IOFormat(Eigen::StreamPrecision, 0, ", ", ";\n", " ", "", "[", ";\n]")
#ifdef MKL
    #define EIGEN_USE_MKL_ALL
#endif
 
// Define macros for the architecture information
#if defined(__NVCC__) && defined(MAGMA) && defined(MKL) && !defined(CPU)
    #if !defined(CPU) && !defined(GPU)
        #define GPU
    #endif
    #if defined(CPU) && !defined(DOUBLE)
        #define DOUBLE
    #endif
    #ifdef DOUBLE
typedef double Real_t;
    #else
        #ifndef FLOAT
            #define FLOAT
        #endif
typedef float Real_t;
    #endif
 
#else
    #ifndef CPU
        #define CPU
    #endif
    #undef GPU
    #undef FLOAT
    #ifndef DOUBLE
        #define DOUBLE
    #endif
typedef double Real_t;
    #define __host__
    #define __device__
    #define __global__
    #define __constant__
#endif  // (END) defined(__NVCC__) && defined(MAGMA) && defined(MKL)
 
// Define complex number types
#ifdef GPU
    #include <cuda/std/complex>
    #include <magma_types.h>
class libFloatComplex : public magmaFloatComplex {
    public:
        using magmaFloatComplex::magmaFloatComplex;
        template<typename RealScalar>
        __host__ __device__ libFloatComplex(RealScalar real, RealScalar imag) {
            this->x = static_cast<float>(real);
            this->y = static_cast<float>(imag);
        }
        __host__ __device__ float& real() { return this->x; }
        __host__ __device__ float  real() const { return this->x; }
        __host__ __device__ float& imag() { return this->y; }
        __host__ __device__ float  imag() const { return this->y; }
};
class libDoubleComplex : public magmaDoubleComplex {
    public:
        using magmaDoubleComplex::magmaDoubleComplex;
        template<typename RealScalar>
        __host__ __device__ libDoubleComplex(RealScalar real, RealScalar imag) {
            this->x = static_cast<double>(real);
            this->y = static_cast<double>(imag);
        }
        __host__ __device__ double& real() { return this->x; }
        __host__ __device__ double  real() const { return this->x; }
        __host__ __device__ double& imag() { return this->y; }
        __host__ __device__ double  imag() const { return this->y; }
};
template<typename T>
using myComplex       = cuda::std::complex<T>;
using myFloatComplex  = cuda::std::complex<float>;
using myDoubleComplex = cuda::std::complex<double>;
 
#elif !defined(MKL) && defined(ACCELERATE)
    #include <complex>
    #include <Accelerate/Accelerate.h>
class libFloatComplex : public __CLPK_complex {
    public:
        using __CLPK_complex::__CLPK_complex;
        __host__ __device__ float& real() { return this->r; }
        __host__ __device__ float  real() const { return this->r; }
        __host__ __device__ float& imag() { return this->i; }
        __host__ __device__ float  imag() const { return this->i; }
};
class libDoubleComplex : public __CLPK_doublecomplex {
    public:
        using __CLPK_doublecomplex::__CLPK_doublecomplex;
        __host__ __device__ double& real() { return this->r; }
        __host__ __device__ double  real() const { return this->r; }
        __host__ __device__ double& imag() { return this->i; }
        __host__ __device__ double  imag() const { return this->i; }
};
template<typename T>
using myComplex       = std::complex<T>;
using myFloatComplex  = std::complex<float>;
using myDoubleComplex = std::complex<double>;
 
#else
    #include <complex>
using libFloatComplex  = std::complex<float>;
using libDoubleComplex = std::complex<double>;
template<typename T>
using myComplex       = std::complex<T>;
using myFloatComplex  = libFloatComplex;
using myDoubleComplex = libDoubleComplex;
#endif  // (END) GPU
 
template<class T = Real_t>
using complex = std::conditional_t<
    std::is_same_v<T, float>, myFloatComplex,
    std::conditional_t<std::is_same_v<T, double>, myDoubleComplex, myComplex<T>> >;
template<class T>
class is_complex {
        template<class U, void* dummy = (&U::imag, nullptr)>
        static std::true_type  test(U*);
        static std::false_type test(...);
 
    public:
        static constexpr bool value = decltype(test((T*)nullptr))::value;
};
template<class T>
inline constexpr bool is_complex_v = is_complex<T>::value;
 
template<class T = Real_t>
class Complex_t : public complex<T> {
    public:
        using complex<T>::complex;
 
        __host__ __device__ Complex_t(complex<T> const& z) : complex<T>(z) {}
 
        template<class U>
        __host__ __device__ Complex_t(Complex_t<U> const& z) : complex<T>(z.real(), z.imag()) {}
 
#if defined(GPU) && defined(MAGMA)
        __host__ __device__ Complex_t(magmaFloatComplex const& z) : complex<T>(z.x, z.y) {}
        __host__ __device__ Complex_t(magmaDoubleComplex const& z) : complex<T>(z.x, z.y) {}
#endif
#if defined(GPU) || defined(ACCELERATE)
        __host__ __device__ Complex_t(libFloatComplex const& z) : complex<T>(z.real(), z.imag()) {}
        __host__ __device__ Complex_t(libDoubleComplex const& z) : complex<T>(z.real(), z.imag()) {}
        __host__ __device__ operator libFloatComplex() const {
            return libFloatComplex(this->real(), this->imag());
        }
        __host__ __device__ operator libDoubleComplex() const {
            return libDoubleComplex(this->real(), this->imag());
        }
#endif  // (END) defined(GPU) || defined(ACCELERATE)
 
};  // (END) class Complex_t
#undef MKL_Complex8
#undef MKL_Complex16
#define MKL_Complex8  Complex_t<float>
#define MKL_Complex16 Complex_t<double>
 
#ifdef GPU
using cuda::std::conj;
using cuda::std::exp;
using cuda::std::fabs;
using cuda::std::imag;
using cuda::std::real;
#else
using std::conj;
using std::exp;
using std::fabs;
using std::imag;
using std::real;
#endif  // (END) #ifdef GPU
 
#if defined(GPU) || defined(ACCELERATE)
__host__ __device__ static inline Complex_t<float> exp(libFloatComplex const& z) {
    return exp(complex<float>(z.real(), z.imag()));
}
__host__ __device__ static inline Complex_t<double> exp(libDoubleComplex const& z) {
    return exp(complex<double>(z.real(), z.imag()));
}
#endif  // (END) #if defined(GPU) || defined(ACCELERATE)
 
template<class T = Real_t>
__constant__ constexpr Complex_t<T> ComplexZero = {0, 0};
template<class T = Real_t>
__constant__ constexpr Complex_t<T> ComplexOne = {1, 0};
template<class T = Real_t>
__constant__ constexpr Complex_t<T> ComplexI = {0, 1};
 
#include <iostream>
#include <sstream>
#include <iomanip>
#include <Eigen/Core>
namespace Eigen {
    template<typename T>
    struct NumTraits<Complex_t<T>> : GenericNumTraits<Complex_t<T>> {
            typedef T                              Real;
            typedef typename NumTraits<T>::Literal Literal;
            enum {
                IsComplex             = 1,
                RequireInitialization = NumTraits<Real>::RequireInitialization,
                ReadCost              = 2 * NumTraits<Real>::ReadCost,
                AddCost               = 2 * NumTraits<Real>::AddCost,
                MulCost               = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost
            };
 
            EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR static inline Real epsilon() {
                return NumTraits<Real>::epsilon();
            }
            EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR static inline Real dummy_precision() {
                return NumTraits<Real>::dummy_precision();
            }
            EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR static inline int digits10() {
                return NumTraits<Real>::digits10();
            }
    };
 
    template<typename charT, typename traits, typename RealType>
    std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& strm,
                                                  const Complex_t<RealType>&         c) {
        constexpr size_t  precision = 4;
        constexpr size_t  realWidth = precision + 7;
        constexpr size_t  imagWidth = realWidth;
        constexpr size_t  width     = realWidth + imagWidth + 2;
        constexpr double  epsilon   = 1.0e-8;
        std::stringstream buff("");
        buff << std::showpos << std::scientific << std::setprecision(precision);
        if(std::abs(c.real()) < epsilon && std::abs(c.imag()) < 1.0e-8)
            buff << 0;
        else if(std::abs(c.real()) > epsilon && std::abs(c.imag()) < epsilon) { buff << c.real(); }
        else if(std::abs(c.real()) < epsilon && std::abs(c.imag()) > epsilon) {
            buff << c.imag() << "*i";
        }
        else { buff << c.real() << c.imag() << "*i"; }
        strm << std::setw(width) << buff.str();
        return strm;
    }
}  // namespace Eigen
// (END) Define complex number types
 
// Class for GPU configs and traits
template<class T>
class on_GPU {
        template<class U, void* dummy = (&U::onGPU, nullptr)>
        static std::true_type  test(U*);
        static std::false_type test(...);
 
    public:
        static constexpr bool value = decltype(test((T*)nullptr))::value;
};
template<class T>
inline constexpr bool on_GPU_v = on_GPU<T>::value;
 
#ifdef GPU
class GPUconfig {
    private:
        dim3          m_dimGrid;
        dim3          m_dimBlock;
        size_t        m_shared;
        magma_queue_t m_queue;
 
    public:
        GPUconfig(){};
        GPUconfig(dim3 A, dim3 B, size_t C = 0, magma_queue_t D = NULL)
            : m_dimGrid(A), m_dimBlock(B), m_shared(C), m_queue(D) {}
        GPUconfig(dim3 A, dim3 B, magma_queue_t D)
            : m_dimGrid(A), m_dimBlock(B), m_shared(), m_queue(D) {}
 
        void dimGrid(int x, int y, int z) {
            m_dimGrid.x = x;
            m_dimGrid.y = y;
            m_dimGrid.z = z;
        }
        void dimBlock(int x, int y, int z) {
            m_dimBlock.x = x;
            m_dimBlock.y = y;
            m_dimBlock.z = z;
        }
        void shared(size_t x) { m_shared = x; }
        void queue(magma_queue_t x) { m_queue = x; }
 
        dim3          dimGrid() const { return m_dimGrid; }
        dim3          dimBlock() const { return m_dimBlock; }
        size_t        shared() const { return m_shared; }
        magma_queue_t queue() const { return m_queue; }
        cudaStream_t  stream() const { return magma_queue_get_cuda_stream(m_queue); }
};
 
    #ifndef CHECK
        #include <iostream>
        #define CHECK(err)                                                              \
            do {                                                                        \
                magma_int_t err_ = (err);                                               \
                if(err_ != 0) {                                                         \
                    std::cerr << "Error: " << #err << "\nfailed at " << __FILE__ << ":" \
                              << __LINE__ << ": error " << (long long)err_ << ": "      \
                              << magma_strerror(err_) << std::endl;                     \
                    std::exit(EXIT_FAILURE);                                            \
                }                                                                       \
            } while(0);
    #endif  // ifndef CHECK
    #ifndef cuCHECK
        #define cuCHECK(call)                                                          \
            {                                                                          \
                const cudaError_t error = call;                                        \
                if(error != cudaSuccess) {                                             \
                    printf("Error: %s:%d,  ", __FILE__, __LINE__);                     \
                    printf("code:%d, reason: %s\n", error, cudaGetErrorString(error)); \
                    cudaDeviceSynchronize();                                           \
                    assert(false);                                                     \
                }                                                                      \
            };
    #endif  // ifndef cuCHECK
#endif      // (END) #ifdef GPU, Class for GPU configs and traits
 
// #ifdef GPU
//  #include <thrust/host_vector.h>
//  #include <thrust/device_vector.h>
// template<class T>
// using host_vector = thrust::host_vector<T>;
// template<class T>
// using device_vector = thrust::device_vector<T>;
// #else
// template<class T>
// using host_vector = std::vector<T>;
// #endif
 
#ifdef MKL
    #include "mkl.h"
    #ifndef zgemm_
        #define zgemm_ zgemm
    #endif
    #ifndef zhemv_
        #define zhemv_ zhemv
    #endif
    #ifndef zdotc_
        #define zdotc_ zdotc
    #endif
#else
    #define MKL_INT int
#endif  // (END) #ifdef MKL
typedef MKL_INT Integer_t;