FourierTransform.h 4.66 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
15 16 17 18 19 20 21 22 23 24
// FourierTransform.h
//
////////////////////////////////////////////////////////////////////////////////

#ifndef FOURIERTRANSFORM_H
#define FOURIERTRANSFORM_H

#include <complex>
#include <vector>

25 26 27 28 29 30 31 32 33
#if !( defined(USE_FFTW2) || defined(USE_FFTW3) || defined(USE_ESSL_FFT) || defined(FFT_NOLIB) )
#error "Must define USE_FFTW2, USE_FFTW3, USE_ESSL_FFT or FFT_NOLIB"
#endif

#if defined(USE_FFTW2) && defined(USE_FFTW3)
#error "Cannot define USE_FFTW2 and USE_FFTW3"
#endif

#if USE_FFTW2
34 35 36
#if USE_DFFTW
#include "dfftw.h"
#else
Francois Gygi committed
37 38
#include "fftw.h"
#endif
39
#endif
Francois Gygi committed
40

41 42 43 44 45 46 47
#if USE_FFTW3
#include "fftw3.h"
#if USE_FFTW3MKL
#include "fftw3_mkl.h"
#endif
#endif

48 49 50 51
#if USE_MPI
#include <mpi.h>
#endif

52 53
#include "Timer.h"

Francois Gygi committed
54 55 56 57 58 59
class Basis;

class FourierTransform
{
  private:

60
  MPI_Comm comm_;
61
  const Basis& basis_;
Francois Gygi committed
62 63 64
  int nprocs_, myproc_;

  int np0_,np1_,np2_;
65
  int ntrans0_,ntrans1_,ntrans2_;
66

Francois Gygi committed
67
  int nvec_;
68

69 70 71
  std::vector<int> np2_loc_; // np2_loc_[iproc], iproc=0, nprocs_-1
  std::vector<int> np2_first_; // np2_first_[iproc], iproc=0, nprocs_-1
  std::vector<std::complex<double> > zvec_;
72

73 74
  std::vector<int> scounts, sdispl, rcounts, rdispl;
  std::vector<std::complex<double> > sbuf, rbuf;
Francois Gygi committed
75

76 77
  std::vector<int> ifftp_, ifftm_;
  std::vector<int> ipack_, iunpack_;
78

Francois Gygi committed
79
  void init_lib(void);
80

81
#if USE_ESSL_FFT
82
#if USE_ESSL_2DFFT
83 84 85 86
  std::vector<double> aux1xyf,aux1zf;
  std::vector<double> aux1xyb,aux1zb;
  std::vector<double> aux2;
  int naux1xy,naux1z,naux2;
Francois Gygi committed
87
#else
88 89 90
  std::vector<double> aux1xf, aux1yf, aux1zf;
  std::vector<double> aux1xb, aux1yb, aux1zb;
  std::vector<double> aux2;
Francois Gygi committed
91 92
  int naux1x,naux1y,naux1z,naux2;
#endif
93
#elif USE_FFTW2
Francois Gygi committed
94
  fftw_plan fwplan0,fwplan1,fwplan2,bwplan0,bwplan1,bwplan2;
95 96 97 98 99
#elif USE_FFTW3
  //plans for np2_
  fftw_plan fwplan, bwplan;
#if defined(USE_FFTW3_2D) || defined(USE_FFTW3_THREADS)
  fftw_plan fwplan2d, bwplan2d;
Francois Gygi committed
100
#else
101 102 103
  fftw_plan fwplanx, fwplany, bwplanx, bwplany;
#endif
#elif defined(FFT_NOLIB)
Francois Gygi committed
104
  // no library
105 106
#else
#error "Must define USE_FFTW2, USE_FFTW3, USE_ESSL_FFT or FFT_NOLIB"
Francois Gygi committed
107 108
#endif

109 110 111 112 113 114 115
  void vector_to_zvec(const std::complex<double>* c);
  void zvec_to_vector(std::complex<double>* c);
  void doublevector_to_zvec(const std::complex<double>* c1,
       const std::complex<double> *c2);
  void zvec_to_doublevector(std::complex<double>* c1, std::complex<double>* c2);
  void fwd(std::complex<double>* val);
  void bwd(std::complex<double>* val);
116

Francois Gygi committed
117 118 119 120
  public:

  FourierTransform (const Basis &basis, int np0, int np1, int np2);
  ~FourierTransform ();
121
  MPI_Comm comm(void) const { return comm_; }
122

Francois Gygi committed
123 124 125 126
  // backward: Fourier synthesis, compute real-space function
  // forward:  Fourier analysis, compute Fourier coefficients
  // forward transform includes scaling by 1/np012
  // single transforms: c -> f, f -> c
127
  void backward (const std::complex<double>* c, std::complex<double>* f);
Francois Gygi committed
128
  // Note: forward transforms overwrite the array f
129
  void forward(std::complex<double>* f, std::complex<double>* c);
130

Francois Gygi committed
131
  // double transforms: c1 + i*c2 -> f, f -> c1 + i*c2
132
  void backward (const std::complex<double>* c1,
133
                 const std::complex<double>* c2, std::complex<double>* f);
Francois Gygi committed
134
  // Note: forward transforms overwrite the array f
135
  void forward(std::complex<double>* f,
136
               std::complex<double>* c1, std::complex<double>* c2);
137

Francois Gygi committed
138 139 140 141 142 143 144 145 146 147 148
  int np0() const { return np0_; }
  int np1() const { return np1_; }
  int np2() const { return np2_; }
  int np2_loc() const { return np2_loc_[myproc_]; }
  int np2_loc(int iproc) const { return np2_loc_[iproc]; }
  int np2_first() const { return np2_first_[myproc_]; }
  int np2_first(int iproc) const { return np2_first_[iproc]; }
  int np012() const { return np0_ * np1_ * np2_; }
  int np012loc(int iproc) const { return np0_ * np1_ * np2_loc_[iproc]; }
  int np012loc() const { return np0_ * np1_ * np2_loc_[myproc_]; }
  int index(int i, int j, int k) const
149
  { return i + np0_ * ( j +  np1_ * k ); }
Francois Gygi committed
150

151 152
  void reset_timers(void);
  Timer tm_f_map, tm_f_fft, tm_f_pack, tm_f_mpi, tm_f_zero, tm_f_unpack,
153 154 155 156 157
        tm_b_map, tm_b_fft, tm_b_pack, tm_b_mpi, tm_b_zero, tm_b_unpack,
        tm_f_xy, tm_f_z, tm_f_x, tm_f_y,
        tm_b_xy, tm_b_z, tm_b_x, tm_b_y,
        tm_init, tm_b_com, tm_f_com;

Francois Gygi committed
158 159
};
#endif