complex.hpp

Go to the documentation of this file.

#ifndef PARGEMSLR_COMPLEX_H
#define PARGEMSLR_COMPLEX_H
 
#include <complex>
#ifdef PARGEMSLR_CUDA
#include <cuda_runtime.h>
#endif
//#include <iostream>
 
namespace pargemslr
{
   
   template <typename T>
   class ComplexValueClass
   {
   private:
      
      T  _real;
      
      T  _imag;
      
   public:
      /* constructor */
      ComplexValueClass()
      {
         _real = T(); 
         _imag = T();
      }
      
      ComplexValueClass(T real)
      {
         _real = real;
         _imag = T();
      }
      
      ComplexValueClass(T real, T imag)
      { 
         _real = real;
         _imag = imag;
      }
      
#ifdef PARGEMSLR_CUDA
__host__ __device__
#endif
      T& Real()
      {
         return _real;
      }
      
#ifdef PARGEMSLR_CUDA
__host__ __device__
#endif
      const T& Real() const
      {
         return _real;
      }
      
#ifdef PARGEMSLR_CUDA
__host__ __device__
#endif
      T& Imag()
      {
         return _imag;
      }
      
#ifdef PARGEMSLR_CUDA
__host__ __device__
#endif
      const T& Imag() const
      {
         return _imag;
      }
      
      T Norm() const
      {
         return _real * _real + _imag * _imag;
      }
      
      T Norm2() const
      {
         return std::sqrt(_real * _real + _imag * _imag);
      }
      
      T Abs() const
      {
         return std::sqrt(_real * _real + _imag * _imag);
      }
      
      ComplexValueClass<T> Conj() const
      {
         return ComplexValueClass<T>(_real, -_imag);
      }
      
      /*
      operator T() const
      {
         if(_imag != T())
         {
            std::cout<<"Error: converting complex type to real type with nonzero imag part, imag part lost."<<std::endl;            
            std::cout<<"This message should never be printed from a correct program."<<std::endl;            
            std::cout<<"Possible cause: operating real value with complex value."<<std::endl;            
         }
         return _real;
      }
      */
      
#ifdef PARGEMSLR_CUDA
__host__ __device__
#endif
      ComplexValueClass<T>& operator=(const T& real)
      {
         _real = real; 
         _imag = T(); 
         return *this;
      }
      
      ComplexValueClass<T>& operator+=(const T& real)
      { 
         _real += real; 
         return *this;
      }
      
      ComplexValueClass<T>& operator-=(const T& real)
      { 
         _real -= real; 
         return *this;
      }
      
      ComplexValueClass<T>& operator*=(const T& real)
      { 
         _real *= real; 
         _imag *= real; 
         return *this;
      }
      
      ComplexValueClass<T>& operator/=(const T& real)
      { 
         _real /= real; 
         _imag /= real; 
         return *this;
      }
      
#ifdef PARGEMSLR_CUDA
__host__ __device__
#endif
      ComplexValueClass<T>& operator=(const ComplexValueClass<T>& val)
      { 
         _real = val.Real(); 
         _imag = val.Imag(); 
         return *this;
      }
      
      ComplexValueClass<T>& operator+=(const ComplexValueClass<T>& val)
      { 
         _real += val.Real(); 
         _imag += val.Imag(); 
         return *this;
      }
      
      ComplexValueClass<T>& operator-=(const ComplexValueClass<T>& val)
      { 
         _real -= val.Real(); 
         _imag -= val.Imag(); 
         return *this;
      }
      
      ComplexValueClass<T>& operator*=(const ComplexValueClass<T>& val)
      { 
         /* (a+bi)*(c+di) = ac-bd + (ad+bc)i */
         const T new_real = _real * val.Real() - _imag * val.Imag();
         _imag = _real * val.Imag() + _imag * val.Real();
         _real = new_real; 
         return *this;
      }
      
      ComplexValueClass<T>& operator/=(const ComplexValueClass<T>& val)
      { 
         /* (a+bi)/(c+di) is slightly more complex
          * [(a+bi)*(c-di)]/[(c+di)*(c-di)] =
          * [ ac+bd + (bc-ad)i ]/(c*c+d*d)
          */
         const T new_lower = val.Norm();
         const T new_real = _real * val.Real() + _imag * val.Imag();
         _imag = (val.Real() * _imag - val.Imag() * _real) / new_lower;
         _real = new_real / new_lower; 
         return *this;
      }
      
   };
   
   template <typename T>
   inline ComplexValueClass<T> operator+(const ComplexValueClass<T> &val, const T& real)
   {
      ComplexValueClass<T> val1 = val;
      val1 += real; 
      return val1;
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator+( const T& real, const ComplexValueClass<T> &val)
   {
      ComplexValueClass<T> val1 = val;
      val1 += real; 
      return val1;
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator+( const ComplexValueClass<T> &val1, const ComplexValueClass<T> &val2)
   {
      ComplexValueClass<T> val = val1;
      val += val2; 
      return val;
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator+( const ComplexValueClass<T> &val)
   {
      return ComplexValueClass<T>( val.Real(), val.Imag());
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator-(const ComplexValueClass<T> &val, const T& real)
   {
      ComplexValueClass<T> val1 = val;
      val1 -= real; 
      return val1;
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator-( const T& real, const ComplexValueClass<T> &val)
   {
      ComplexValueClass<T> val1 = ComplexValueClass<T>(real, T());
      val1 -= val; 
      return val1;
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator-( const ComplexValueClass<T> &val1, const ComplexValueClass<T> &val2)
   {
      ComplexValueClass<T> val = val1;
      val -= val2; 
      return val;
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator-( const ComplexValueClass<T> &val)
   {
      return ComplexValueClass<T>( -val.Real(), -val.Imag());
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator*(const ComplexValueClass<T> &val, const T& real)
   {
      ComplexValueClass<T> val1 = val;
      val1 *= real; 
      return val1;
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator*( const T& real, const ComplexValueClass<T> &val)
   {
      ComplexValueClass<T> val1 = val;
      val1 *= real; 
      return val1;
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator*( const ComplexValueClass<T> &val1, const ComplexValueClass<T> &val2)
   {
      ComplexValueClass<T> val = val1;
      val *= val2; 
      return val;
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator/(const ComplexValueClass<T> &val, const T& real)
   {
      ComplexValueClass<T> val1 = val;
      val1 /= real; 
      return val1;
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator/( const T& real, const ComplexValueClass<T> &val)
   {
      ComplexValueClass<T> val1 = ComplexValueClass<T>(real, T());
      val1 /= val; 
      return val1;
   }
   
   template <typename T>
   inline ComplexValueClass<T> operator/( const ComplexValueClass<T> &val1, const ComplexValueClass<T> &val2)
   {
      ComplexValueClass<T> val = val1;
      val /= val2; 
      return val;
   }
   
   template <typename T>
   inline bool operator==( const ComplexValueClass<T> &val1, const T &val2)
   {
      return (val1.Real() == val2) && (val1.Imag() == T());
   }
   
   template <typename T>
   inline bool operator==( const T &val1, const ComplexValueClass<T> &val2)
   {
      return (val1 == val2.Real()) && (T() == val2.Imag());
   }
   
   template <typename T>
   inline bool operator==( const ComplexValueClass<T> &val1, const ComplexValueClass<T> &val2)
   {
      return (val1.Real() == val2.Real()) && (val1.Imag() == val2.Imag());
   }
   
   template <typename T>
   inline bool operator!=( const ComplexValueClass<T> &val1, const T &val2)
   {
      return (val1.Real() != val2) || (val1.Imag() != T());
   }
   
   template <typename T>
   inline bool operator!=( const T &val1, const ComplexValueClass<T> &val2)
   {
      return (val1 != val2.Real()) || (T() != val2.Imag());
   }
   
   template <typename T>
   inline bool operator!=( const ComplexValueClass<T> &val1, const ComplexValueClass<T> &val2)
   {
      return (val1.Real() != val2.Real()) || (val1.Imag() != val2.Imag());
   }
   
   template<typename T, typename CharT, class Traits>
   std::basic_ostream<CharT, Traits>&
   operator<<(std::basic_ostream<CharT, Traits>& basic_os, const ComplexValueClass<T>& val)
   {
      std::basic_ostringstream<CharT, Traits> basic_ostr;
      basic_ostr.flags(basic_os.flags());
      basic_ostr.imbue(basic_os.getloc());
      basic_ostr.precision(basic_os.precision());
      
      T val_r = val.Real();
      T val_i = val.Imag();
      
      if(val_r < 0)
      {
         basic_ostr<<"-";
      }
      else
      {
         basic_ostr<<"+";
      }
      basic_ostr<<std::abs(val_r);
      if(val_i < 0)
      {
         basic_ostr<<"-";
      }
      else
      {
         basic_ostr<<"+";
      }
      basic_ostr<<std::abs(val_i)<<"i";
      return basic_os << basic_ostr.str();
   }
   
   typedef ComplexValueClass<float>    complexs;
   typedef ComplexValueClass<double>   complexd;
 
#ifdef PARGEMSLR_OPENMP
 
#pragma omp declare reduction(+:complexs:omp_out += omp_in) initializer (omp_priv=complexs())
#pragma omp declare reduction(+:complexd:omp_out += omp_in) initializer (omp_priv=complexd())
 
#endif
 
   typedef struct CComplexSingleStruct
   {
      float real, imag;
   }ccomplexs;
   
   typedef struct CComplexDoubleStruct
   {
      double real, imag;
   }ccomplexd;
   
   template <class T> struct PargemslrIsComplex : public std::false_type {};
   template <class T> struct PargemslrIsComplex<const T> : public PargemslrIsComplex<T> {};
   template <class T> struct PargemslrIsComplex<volatile const T> : public PargemslrIsComplex<T>{};
   template <class T> struct PargemslrIsComplex<volatile T> : public PargemslrIsComplex<T>{};
   template<> struct PargemslrIsComplex<complexs> : public std::true_type {};
   template<> struct PargemslrIsComplex<complexd> : public std::true_type {};
   
   template <class T> struct PargemslrIsReal : public std::false_type {};
   template <class T> struct PargemslrIsReal<const T> : public PargemslrIsReal<T> {};
   template <class T> struct PargemslrIsReal<volatile const T> : public PargemslrIsReal<T>{};
   template <class T> struct PargemslrIsReal<volatile T> : public PargemslrIsReal<T>{};
   template<> struct PargemslrIsReal<float> : public std::true_type {};
   template<> struct PargemslrIsReal<double> : public std::true_type {};
   
}
 
#endif