prdc/cpu/complex.h

#ifndef PRDC_CPU_COMPLEX_H
#define PRDC_CPU_COMPLEX_H
 
/*
* PSCF Package - Polymer Self-Consistent Field Theory
*
* Copyright 2016 - 2022, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/
 
#include <pscf/math/complex.h>
 
#include <fftw3.h>
#include <complex>
#include <iostream>
 
namespace Pscf {
 
   /*
   * Types Cpu::Complex and Cpu::Real are defined in prdc/cpu/types.h
   * as aliases for fftw_complex and double, respectively.
   */
 
   // Real and imaginary components
 
   template <> inline 
   double real<fftw_complex, double>(fftw_complex const& a)
   {  return a[0]; }
 
   template <> inline 
   double imag<fftw_complex, double>(fftw_complex const& a)
   {  return a[1]; }
 
   // Absolute magnitude
 
   template <> inline 
   double abs<fftw_complex, double>(fftw_complex const& a)
   {  return sqrt(a[0] * a[0] + a[1] * a[1]); }
 
   template <> inline 
   double absSq<fftw_complex, double>(fftw_complex const& a)
   {  return (a[0] * a[0] + a[1] * a[1]); }
 
   // fftw_complex Conjugation
 
   template <> inline
   void conj(fftw_complex& z, fftw_complex const& a)
   {
      z[0] = a[0];
      z[1] = -a[1];
   }
 
   template <> inline
   void conj(fftw_complex& a)
   {
      a[0] = a[0];
      a[1] = -a[1];
   }
 
   // Assignment
 
   template <> inline
   void assign(fftw_complex& z, double const& a, double const& b)
   {
      z[0] = a;
      z[1] = b;
   }
 
   template <> inline
   void assign(fftw_complex& z, double const& a)
   {
      z[0] = a;
      z[1] = 0;
   }
 
   template <> inline
   void assign(fftw_complex& z, fftw_complex const& a)
   {
      z[0] = a[0];
      z[1] = a[1];
   }
 
   template <> inline
   void assign(fftw_complex & z, std::complex<double> const& a)
   {
      z[0] = a.real();
      z[1] = a.imag();
   }
 
   template <> inline
   void assign(std::complex<double> & z, fftw_complex const& a)
   {  z = std::complex<double>(a[0], a[1]); }
 
   // Addition
 
   template <> inline
   void add(fftw_complex& z, fftw_complex const& a, fftw_complex const& b)
   {
      z[0] = a[0] + b[0];
      z[1] = a[1] + b[1];
   }
 
   template <> inline
   void add(fftw_complex& z, fftw_complex const& a, double const& b)
   {
      z[0] = a[0] + b;
      z[1] = a[1];
   }
 
   template <> inline
   void addEq(fftw_complex& a, fftw_complex const& b)
   {
      a[0] += b[0];
      a[1] += b[1];
   }
 
   template <> inline
   void addEq(fftw_complex& a, double const& b)
   {
      a[0] += b;
   }
 
   // Subtraction
 
   template <> inline
   void sub(fftw_complex& z, fftw_complex const& a, fftw_complex const& b)
   {
      z[0] = a[0] - b[0];
      z[1] = a[1] - b[1];
   }
 
   template <> inline
   void sub(fftw_complex& z, fftw_complex const& a, double const& b)
   {
      z[0] = a[0] - b;
      z[1] = a[1];
   }
 
   template <> inline
   void subEq(fftw_complex & a, fftw_complex const& b)
   {
      a[0] -= b[0];
      a[1] -= b[1];
   }
 
   template <> inline
   void subEq(fftw_complex & a, double const& b)
   {
      a[0] -= b;
   }
 
   template <> inline 
   double absSqDiff(fftw_complex const& a, fftw_complex const& b)
   {
      fftw_complex z;
      sub(z, a, b);
      return absSq<fftw_complex, double>(z);
   }
 
   // Multiplication
 
   template <> inline
   void mul(fftw_complex& z, fftw_complex const& a, fftw_complex const& b)
   {
      z[0] = a[0] * b[0] - a[1] * b[1];
      z[1] = a[1] * b[0] + a[0] * b[1];
   }
 
   template <> inline
   void mul(fftw_complex& z, fftw_complex const& a, double const& b)
   {
      z[0] = a[0]*b;
      z[1] = a[1]*b;
   }
 
   template <> inline
   void mulEq(fftw_complex & a, fftw_complex const& b)
   {
      double a0;
      a0   = a[0] * b[0] - a[1] * b[1];
      a[1] = a[1] * b[0] + a[0] * b[1];
      a[0] = a0;
   }
 
   template <> inline
   void mulEq(fftw_complex & a, double const& b)
   {
      a[0] *= b;
      a[1] *= b;
   }
 
   template <> inline
   void square(fftw_complex& z, fftw_complex const& a)
   {
      z[0] = a[0] * a[0] - a[1] * a[1];
      z[1] = 2.0 * a[1] * a[0];
   }
 
   // Division
 
   template <> inline
   void div(fftw_complex& z, fftw_complex const& a, fftw_complex const& b)
   {
      double bSq = b[0] * b[0] + b[1] * b[1];
      z[0] = (a[0] * b[0] + a[1] * b[1])/bSq;
      z[1] = (a[1] * b[0] - a[0] * b[1])/bSq;
   }
 
   template <> inline
   void div(fftw_complex& z, fftw_complex const& a, double const& b)
   {
      z[0] = a[0]/b;
      z[1] = a[1]/b;
   }
 
   template <> inline
   void divEq(fftw_complex & a, fftw_complex const & b)
   {
      double bSq = b[0] * b[0] + b[1] * b[1];
      double a0 = (a[0] * b[0] + a[1] * b[1])/bSq;
      a[1] = (a[1] * b[0] - a[0] * b[1])/bSq;
      a[0] = a0;
   }
 
   template <> inline
   void divEq(fftw_complex & a, double const& b)
   {
      a[0] /= b;
      a[1] /= b;
   }
 
   std::istream& operator >> (std::istream& is, fftw_complex & z);
 
   std::ostream& operator << (std::ostream& os, fftw_complex const & z);
 
   #if 0
   template <typename Archive>
   void serialize(Archive& ar, double z[2], const unsigned int version = 0)
   {
      serialize(ar, z[0], version);
      serialize(ar, z[1], version);
   }
   #endif
 
} // namespace Pscf
#endif