1 /*-------------------------------------------------------------------
2 Copyright 2011 Ravishankar Sundararaman
3
4 This file is part of JDFTx.
5
6 JDFTx is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 JDFTx is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with JDFTx. If not, see <http://www.gnu.org/licenses/>.
18 -------------------------------------------------------------------*/
19
20 #include <fluid/TranslationOperator.h>
21 #include <fluid/TranslationOperator_internal.h>
22 #include <core/Operators.h>
23
24
TranslationOperator(const GridInfo & gInfo)25 TranslationOperator::TranslationOperator(const GridInfo& gInfo) : gInfo(gInfo)
26 {
27 }
28
TranslationOperatorSpline(const GridInfo & gInfo,SplineType splineType)29 TranslationOperatorSpline::TranslationOperatorSpline(const GridInfo& gInfo, SplineType splineType)
30 : TranslationOperator(gInfo), splineType(splineType)
31 {
32 }
33
constantSplineTaxpy_sub(size_t iStart,size_t iStop,const vector3<int> S,double alpha,const double * x,double * y,const vector3<int> Tint)34 void constantSplineTaxpy_sub(size_t iStart, size_t iStop, const vector3<int> S,
35 double alpha, const double* x, double* y, const vector3<int> Tint)
36 { THREAD_rLoop(constantSplineTaxpy_calc(i, iv, S, alpha, x, y, Tint);)
37 }
linearSplineTaxpy_sub(size_t iStart,size_t iStop,const vector3<int> S,double alpha,const double * x,double * y,const vector3<int> Tint,const vector3<> Tfrac)38 void linearSplineTaxpy_sub(size_t iStart, size_t iStop, const vector3<int> S,
39 double alpha, const double* x, double* y, const vector3<int> Tint, const vector3<> Tfrac)
40 { THREAD_rLoop(linearSplineTaxpy_calc(i, iv, S, alpha, x, y, Tint, Tfrac);)
41 }
42 #ifdef GPU_ENABLED
43 void constantSplineTaxpy_gpu(const vector3<int> S,
44 double alpha, const double* x, double* y, const vector3<int> Tint);
45 void linearSplineTaxpy_gpu(const vector3<int> S,
46 double alpha, const double* x, double* y, const vector3<int> Tint, const vector3<> Tfrac);
47 #endif
taxpy(const vector3<> & t,double alpha,const ScalarField & x,ScalarField & y) const48 void TranslationOperatorSpline::taxpy(const vector3<>& t, double alpha, const ScalarField& x, ScalarField& y) const
49 { //Perform a gather with the inverse translation (hence negate t),
50 //instead of scatter which is less efficient to parallelize
51 vector3<> Tfrac = Diag(gInfo.S) * inv(gInfo.R) * (-t); //now in grid point units
52 vector3<int> Tint;
53 //Prepare output:
54 nullToZero(y, gInfo);
55 switch(splineType)
56 { case Constant:
57 { for(int k=0; k<3; k++)
58 { //round to nearest integer (and ensure symmetric rounding direction for transpose correctness):
59 Tint[k] = int(copysign(floor(fabs(Tfrac[k])+0.5), Tfrac[k]));
60 //reduce to positive first unit cell:
61 Tint[k] = Tint[k] % gInfo.S[k];
62 if(Tint[k]<0) Tint[k] += gInfo.S[k];
63 }
64 //Launch threads/gpu kernels:
65 #ifdef GPU_ENABLED
66 constantSplineTaxpy_gpu(gInfo.S, alpha*x->scale, x->dataGpu(false), y->dataGpu(), Tint);
67 #else
68 threadLaunch(constantSplineTaxpy_sub, gInfo.nr, gInfo.S, alpha*x->scale, x->data(false), y->data(), Tint);
69 #endif
70 break;
71 }
72 case Linear:
73 { for(int k=0; k<3; k++)
74 { //reduce to positive first unit cell:
75 Tfrac[k] = fmod(Tfrac[k], gInfo.S[k]);
76 if(Tfrac[k]<0) Tfrac[k] += gInfo.S[k];
77 //separate integral and fractional parts:
78 Tint[k] = int(floor(Tfrac[k]));
79 Tfrac[k] -= Tint[k];
80 Tint[k] = Tint[k] % gInfo.S[k];
81 }
82 //Launch threads/gpu kernels:
83 #ifdef GPU_ENABLED
84 linearSplineTaxpy_gpu(gInfo.S, alpha*x->scale, x->dataGpu(false), y->dataGpu(), Tint, Tfrac);
85 #else
86 threadLaunch(linearSplineTaxpy_sub, gInfo.nr, gInfo.S, alpha*x->scale, x->data(false), y->data(), Tint, Tfrac);
87 #endif
88 break;
89 }
90 }
91 }
92
TranslationOperatorFourier(const GridInfo & gInfo)93 TranslationOperatorFourier::TranslationOperatorFourier(const GridInfo& gInfo)
94 : TranslationOperator(gInfo)
95 {
96 }
fourierTranslate_sub(size_t iStart,size_t iStop,const vector3<int> S,const vector3<> Gt,complex * xTilde)97 inline void fourierTranslate_sub(size_t iStart, size_t iStop, const vector3<int> S, const vector3<> Gt, complex* xTilde)
98 { THREAD_halfGspaceLoop( fourierTranslate_calc(i, iG, S, Gt, xTilde); )
99 }
100 #ifdef GPU_ENABLED //implemented in TranslationOperator.cu
101 void fourierTranslate_gpu(const vector3<int> S, const vector3<> Gt, complex* xTilde);
102 #endif
taxpy(const vector3<> & t,double alpha,const ScalarField & x,ScalarField & y) const103 void TranslationOperatorFourier::taxpy(const vector3<>& t, double alpha, const ScalarField& x, ScalarField& y) const
104 { ScalarFieldTilde xTilde = J(x);
105 #ifdef GPU_ENABLED
106 fourierTranslate_gpu(gInfo.S, gInfo.G*t, xTilde->dataGpu(false));
107 #else
108 threadLaunch(fourierTranslate_sub, gInfo.nG, gInfo.S, gInfo.G*t, xTilde->data(false));
109 #endif
110 y += alpha*I(xTilde);
111 }
112