Commit bacb4b94 by Francois Gygi

Test for ntrans=0 case to avoid ESSL DCFT error when performing zero transforms.


git-svn-id: http://qboxcode.org/svn/qb/trunk@1898 cba15fb0-1239-40c8-b417-11db7ca47a34
parent 519b2295
......@@ -589,8 +589,9 @@ void FourierTransform::bwd(complex<double>* val)
int inc1 = 1, inc2 = np2_, ntrans = nvec_, isign = -1, initflag = 0;
double scale = 1.0;
dcft_(&initflag,&zvec_[0],&inc1,&inc2,&zvec_[0],&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zb[0],&naux1z,&aux2[0],&naux2);
if ( ntrans > 0 )
dcft_(&initflag,&zvec_[0],&inc1,&inc2,&zvec_[0],&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zb[0],&naux1z,&aux2[0],&naux2);
#elif USE_FFTW2
/*
* void fftw(fftw_plan plan, int howmany,
......@@ -834,29 +835,35 @@ void FourierTransform::bwd(complex<double>* val)
// transform only non-zero vectors along x
// First block: positive y indices: [0,ntrans0_]
ntrans = ntrans0_;
inc1 = 1;
inc2 = np0_;
istart = k * np0_ * np1_;
length = np0_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1xb[0],&naux1x,&aux2[0],&naux2);
// Second block: negative y indices: [np1-ntrans0_,np1-1]
inc1 = 1;
inc2 = np0_;
istart = np0_ * ( (np1_-ntrans) + k * np1_ );
length = np0_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1xb[0],&naux1x,&aux2[0],&naux2);
if ( ntrans > 0 )
{
inc1 = 1;
inc2 = np0_;
istart = k * np0_ * np1_;
length = np0_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1xb[0],&naux1x,&aux2[0],&naux2);
// Second block: negative y indices: [np1-ntrans0_,np1-1]
inc1 = 1;
inc2 = np0_;
istart = np0_ * ( (np1_-ntrans) + k * np1_ );
length = np0_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1xb[0],&naux1x,&aux2[0],&naux2);
}
// transform along y for all values of x
ntrans = np0_;
inc1 = np0_;
inc2 = 1;
istart = k * np0_ * np1_;
length = np1_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1yb[0],&naux1y,&aux2[0],&naux2);
if ( ntrans > 0 )
{
inc1 = np0_;
inc2 = 1;
istart = k * np0_ * np1_;
length = np1_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1yb[0],&naux1y,&aux2[0],&naux2);
}
#endif // USE_ESSL_2DFFT
} // k
......@@ -1062,28 +1069,34 @@ void FourierTransform::fwd(complex<double>* val)
double scale = 1.0;
// transform along y for all values of x
ntrans = np0_;
inc1 = np0_;
inc2 = 1;
istart = k * np0_ * np1_;
length = np1_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1yf[0],&naux1y,&aux2[0],&naux2);
if ( ntrans > 0 )
{
inc1 = np0_;
inc2 = 1;
istart = k * np0_ * np1_;
length = np1_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1yf[0],&naux1y,&aux2[0],&naux2);
}
// transform only non-zero vectors along x
ntrans = ntrans0_;
inc1 = 1;
inc2 = np0_;
istart = k * np0_ * np1_;
length = np0_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1xf[0],&naux1x,&aux2[0],&naux2);
inc1 = 1;
inc2 = np0_;
istart = np0_ * ( (np1_-ntrans) + k * np1_ );
length = np0_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1xf[0],&naux1x,&aux2[0],&naux2);
if ( ntrans > 0 )
{
inc1 = 1;
inc2 = np0_;
istart = k * np0_ * np1_;
length = np0_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1xf[0],&naux1x,&aux2[0],&naux2);
inc1 = 1;
inc2 = np0_;
istart = np0_ * ( (np1_-ntrans) + k * np1_ );
length = np0_;
dcft_(&initflag,&val[istart],&inc1,&inc2,&val[istart],&inc1,&inc2,
&length,&ntrans,&isign,&scale,&aux1xf[0],&naux1x,&aux2[0],&naux2);
}
#endif // USE_ESSL_2DFFT
} // k
#elif USE_FFTW2
......@@ -1275,8 +1288,9 @@ void FourierTransform::fwd(complex<double>* val)
int inc1 = 1, inc2 = np2_, ntrans = nvec_, isign = 1, initflag = 0;
double scale = 1.0 / np012();
dcft_(&initflag,&zvec_[0],&inc1,&inc2,&zvec_[0],&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zf[0],&naux1z,&aux2[0],&naux2);
if ( ntrans > 0 )
dcft_(&initflag,&zvec_[0],&inc1,&inc2,&zvec_[0],&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zf[0],&naux1z,&aux2[0],&naux2);
#elif USE_FFTW2
#if _OPENMP
......@@ -1386,13 +1400,16 @@ void FourierTransform::init_lib(void)
// initialize z transforms
int ntrans = nvec_;
inc1 = 1; inc2 = np2_;
isign = -1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zb[0],&naux1z,&aux2[0],&naux2);
isign = 1; scale = 1.0 / np012();
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zf[0],&naux1z,&aux2[0],&naux2);
if ( ntrans > 0 )
{
inc1 = 1; inc2 = np2_;
isign = -1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zb[0],&naux1z,&aux2[0],&naux2);
isign = 1; scale = 1.0 / np012();
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zf[0],&naux1z,&aux2[0],&naux2);
}
#else // USE_ESSL_2DFFT
naux1x = (int) (20000 + 2.28 * np0_);
......@@ -1420,30 +1437,39 @@ void FourierTransform::init_lib(void)
// x transforms
inc1 = 1; inc2 = np0_; ntrans = ntrans0_;
isign = -1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np0_,&ntrans,
&isign,&scale,&aux1xb[0],&naux1x,&aux2[0],&naux2);
isign = 1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np0_,&ntrans,
&isign,&scale,&aux1xf[0],&naux1x,&aux2[0],&naux2);
if ( ntrans > 0 )
{
isign = -1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np0_,&ntrans,
&isign,&scale,&aux1xb[0],&naux1x,&aux2[0],&naux2);
isign = 1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np0_,&ntrans,
&isign,&scale,&aux1xf[0],&naux1x,&aux2[0],&naux2);
}
// y transforms
inc1 = np0_; inc2 = 1; ntrans = ntrans1_;
isign = -1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np1_,&ntrans,
&isign,&scale,&aux1yb[0],&naux1y,&aux2[0],&naux2);
isign = 1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np1_,&ntrans,
&isign,&scale,&aux1yf[0],&naux1y,&aux2[0],&naux2);
if ( ntrans > 0 )
{
isign = -1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np1_,&ntrans,
&isign,&scale,&aux1yb[0],&naux1y,&aux2[0],&naux2);
isign = 1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np1_,&ntrans,
&isign,&scale,&aux1yf[0],&naux1y,&aux2[0],&naux2);
}
// z transforms
inc1 = 1; inc2 = np2_; ntrans = ntrans2_;
isign = -1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zb[0],&naux1z,&aux2[0],&naux2);
isign = 1; scale = 1.0 / np012();
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zf[0],&naux1z,&aux2[0],&naux2);
if ( ntrans > 0 )
{
isign = -1;
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zb[0],&naux1z,&aux2[0],&naux2);
isign = 1; scale = 1.0 / np012();
dcft_(&initflag,p,&inc1,&inc2,p,&inc1,&inc2,&np2_,&ntrans,
&isign,&scale,&aux1zf[0],&naux1z,&aux2[0],&naux2);
}
#endif // USE_ESSL_2DFFT
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment