sinft.C 4.7 KB
Newer Older
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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15 16 17
// sinft.C
//
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
18 19 20 21

#include "sinft.h"
#include <math.h>
#include <assert.h>
22 23 24 25 26
#include <vector>
#include <complex>
using namespace std;

#if USE_FFTW2
27 28 29
#if USE_DFFTW
#include "dfftw.h"
#else
30
#include "fftw.h"
31
#endif
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
#elif USE_FFTW3
#include "fftw3.h"
#elif USE_ESSL_FFT
extern "C" {
  void dcft_(int *initflag, std::complex<double> *x, int *inc2x, int *inc3x,
             std::complex<double> *y, int *inc2y, int *inc3y,
             int *length, int *ntrans, int *isign,
             double *scale, double *aux1, int *naux1,
             double *aux2, int *naux2);
}
#else
// no FFT library
void cfft ( int idir, complex<double> *z1, complex<double> *z2, int n,
  int *inzee );
void fftstp ( int idir, complex<double> *zin, int after,
              int now, int before, complex<double> *zout );
48
#endif
Francois Gygi committed
49

50
void sinft(int n, double *y)
Francois Gygi committed
51
{
52 53
  vector<complex<double> > zin(2*n), zout(2*n);
#if defined(USE_FFTW2) || defined(USE_FFTW3)
54
  fftw_plan fwplan;
55 56
#endif
#if USE_FFTW2
57
  fwplan = fftw_create_plan(2*n,FFTW_FORWARD,FFTW_ESTIMATE);
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
#elif USE_FFTW3
  fwplan = fftw_plan_dft_1d(2*n,(fftw_complex*)&zin[0],(fftw_complex*)&zout[0],
                            FFTW_FORWARD, FFTW_ESTIMATE);
#elif USE_ESSL_FFT
  int np = 2 * n;
  int naux1 = (int) (30000 + 2.28 * np);
  std::vector<double> aux1(naux1);
  int ntrans = 1;
  int naux2 = (int) (20000 + 2.28 * np + (256 + 2*np)*min(64,ntrans));
  std::vector<double> aux2(naux2);
#else
  // no FFT library
  // no initialization needed
#endif

73 74
  zin[0] = 0.0;
  for ( int i = 1; i < n; i++ )
Francois Gygi committed
75
  {
76 77 78
    const double t = y[i];
    zin[i] = t;
    zin[2*n-i] = -t;
Francois Gygi committed
79
  }
80
#if USE_FFTW2
81
  fftw_one(fwplan,(fftw_complex*)&zin[0],(fftw_complex*)&zout[0]);
82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
#elif USE_FFTW3
  fftw_execute(fwplan);
#elif USE_ESSL_FFT
  // initialize forward transform
  int initflag = 1;
  int inc1 = 1, inc2 = np;
  int isign = 1;
  double scale = 1.0;
  complex<double> *p = 0;
  dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np,&ntrans,
        &isign,&scale,&aux1[0],&naux1,&aux2[0],&naux2);

  // call transform
  initflag = 0;
  dcft_(&initflag,&zin[0],&inc1,&inc2,&zout[0],&inc1,&inc2,&np,&ntrans,
        &isign,&scale,&aux1[0],&naux1,&aux2[0],&naux2);
#else
  // no FFT library
  int idir = 1, inzee = 1, np = 2*n;
  cfft ( idir, &zin[0],&zout[0],np, &inzee );
#endif
103
  for ( int i = 0; i < n; i++ )
Francois Gygi committed
104
  {
105
    y[i] = -0.5 * imag(zout[i]);
Francois Gygi committed
106
  }
107
#if defined(USE_FFTW2) || defined(USE_FFTW3)
108
  fftw_destroy_plan(fwplan);
109
#endif
Francois Gygi committed
110 111
}

112
void cosft1(int n, double *y)
Francois Gygi committed
113
{
114
  /* Note: the array y contains n+1 elements */
115 116
  vector<complex<double> > zin(2*n), zout(2*n);
#if defined(USE_FFTW2) || defined(USE_FFTW3)
117
  fftw_plan fwplan;
118 119
#endif
#if USE_FFTW2
120
  fwplan = fftw_create_plan(2*n,FFTW_FORWARD,FFTW_ESTIMATE);
121 122 123 124 125 126 127 128 129 130 131 132 133 134
#elif USE_FFTW3
  fwplan = fftw_plan_dft_1d(2*n,(fftw_complex*)&zin[0],(fftw_complex*)&zout[0],
                            FFTW_FORWARD, FFTW_ESTIMATE);
#elif USE_ESSL_FFT
  int np = 2 * n;
  int naux1 = (int) (30000 + 2.28 * np);
  std::vector<double> aux1(naux1);
  int ntrans = 1;
  int naux2 = (int) (20000 + 2.28 * np + (256 + 2*np)*min(64,ntrans));
  std::vector<double> aux2(naux2);
#else
  // no FFT library
  // no initialization needed
#endif
Francois Gygi committed
135

136 137
  zin[0] = y[0];
  for ( int i = 1; i < n+1; i++ )
Francois Gygi committed
138
  {
139 140 141
    const double t = y[i];
    zin[i] = t;
    zin[2*n-i] = t;
Francois Gygi committed
142
  }
143
#if USE_FFTW2
144
  fftw_one(fwplan,(fftw_complex*)&zin[0],(fftw_complex*)&zout[0]);
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
#elif USE_FFTW3
  fftw_execute(fwplan);
#elif USE_ESSL_FFT
  // initialize forward transform
  int initflag = 1;
  int inc1 = 1, inc2 = np;
  int isign = 1;
  double scale = 1.0;
  complex<double> *p = 0;
  dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np,&ntrans,
        &isign,&scale,&aux1[0],&naux1,&aux2[0],&naux2);

  // call transform
  initflag = 0;
  dcft_(&initflag,&zin[0],&inc1,&inc2,&zout[0],&inc1,&inc2,&np,&ntrans,
        &isign,&scale,&aux1[0],&naux1,&aux2[0],&naux2);
#else
  // no FFT library
  int idir = 1, inzee = 1, np = 2*n;
  cfft ( idir, &zin[0],&zout[0],np, &inzee );
#endif
166 167
  y[0] = 0.5 * real(zout[0]);
  for ( int i = 1; i < n; i++ )
Francois Gygi committed
168
  {
169
    y[i] = 0.5 * real(zout[i]);
Francois Gygi committed
170
  }
171
#if defined(USE_FFTW2) || defined(USE_FFTW3)
172
  fftw_destroy_plan(fwplan);
173
#endif
Francois Gygi committed
174
}