FourierTransform.h 4.87 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 48 49 50 51
#if USE_FFTW3
#include "fftw3.h"
#if USE_FFTW3MKL
#include "fftw3_mkl.h"
#endif
#endif

#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
  bool basis_fits_in_grid_;
69

70 71 72
  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_;
73

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

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

Francois Gygi committed
80
  void init_lib(void);
81

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

110 111 112 113 114 115 116
  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);
117

Francois Gygi committed
118 119 120 121
  public:

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

Francois Gygi committed
124 125 126 127
  // 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
128
  void backward (const std::complex<double>* c, std::complex<double>* f);
Francois Gygi committed
129
  // Note: forward transforms overwrite the array f
130
  void forward(std::complex<double>* f, std::complex<double>* c);
131

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

Francois Gygi committed
139 140 141 142 143 144 145
  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]; }
Francois Gygi committed
146
  long int np012() const { return ((long int)np0_) * np1_ * np2_; }
Francois Gygi committed
147 148 149
  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
150
  { return i + np0_ * ( j +  np1_ * k ); }
151 152 153
  int i(int ind) const { return ind % np0_; }
  int j(int ind) const { return (ind / np0_) % np1_; }
  int k(int ind) const { return (ind / np0_) / np1_ + np2_first(); }
Francois Gygi committed
154

155 156
  void reset_timers(void);
  Timer tm_f_map, tm_f_fft, tm_f_pack, tm_f_mpi, tm_f_zero, tm_f_unpack,
157 158 159 160 161
        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
162 163
};
#endif