Commit ec447c3c authored by Florian Goth's avatar Florian Goth
Browse files

current things

parent 86183979
Pipeline #9979 failed with stage
in 5 minutes and 56 seconds
......@@ -157,14 +157,20 @@ class SSH_multisite
int L;
double* gdat1D;
std::complex<double>* dexpk;
double* ftexp;
int ntau;
SSH_multisite(double g, double b, double d, int myL) : k(-0.5*g*g), beta(b), dtau(d), L(myL), g(g)
{
ntau = std::round(beta/dtau); // ntau should be nothing else than Ltime
gdat1D = new double[ntau*L];
dexpk = new std::complex<double>[L];
ftexp = new double[L*L];
for(int d = 0; d < L; ++d)
{
dexpk[d] = std::exp(std::complex<double>(0.0, -2.0*M_PI/L*d ));
for(int k = 0; k < L; ++k)
ftexp[d*L + k] = std::real(std::exp(std::complex<double>(0.0, -2.0*M_PI*k/L*d )));
}
for(int j = 0; j < L; ++j)
{
......@@ -293,7 +299,6 @@ class SSH_multisite
template <typename VertexType>
double green(const VertexType& c1, const VertexType& c2)
{
std::complex<double> retval = 0;
std::complex<double> jj(0,1);
// space
......@@ -322,20 +327,23 @@ class SSH_multisite
// compute Greens function in Fourier space
std::complex<double> dexpkx = dexpk[dist];
std::complex<double> expk = 1.0;
const double *const __restrict__ dat = ftexp + dist*L;
const double *const __restrict__ gdat = gdat1D + dti*L;
double retval = 0;
// #pragma omp simd
for (int j = 0; j < L; ++j)
{
retval = std::fma(dat[j], gdat[j], retval);
// perform sum
retval += expk.real()*gdat1D[dti*L + j];
// increment Fourier step
expk = std::complex<double>(
expk.real() * dexpkx.real() - expk.imag() * dexpkx.imag(),
expk.real() * dexpkx.imag() + expk.imag() * dexpkx.real()
);
// expk *= dexpkx; //loop-carried dependency. breaks vectorization.
// retval[0] = std::fma(dat[j], gdat[j], retval[0]);//currently fastest and simplest on my CPU...
// retval[1] = std::fma(dat[j+1], gdat[j+1], retval[1]);//currently fastest and simplest on my CPU...
// retval[2] = std::fma(dat[j+2], gdat[j+2], retval[2]);//currently fastest and simplest on my CPU...
// retval[3] = std::fma(dat[j+3], gdat[j+3], retval[3]);//currently fastest and simplest on my CPU...
// retval += dat[j] * gdat[j];
}
const double norml = 1.0/L;
return norml*signum*retval.real();
return norml*signum*(retval /*+ retval[1] + retval[2] + retval[3]*/);
}
#endif
......
Supports Markdown
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