portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-05-10 12:35:57 +01:00
Code Releases Activity

Stenccil improvement

Browse Source
This commit is contained in:
Peter Boyle 2020-07-16 16:57:13 -04:00
parent 488c79d5a1
commit 569f78c2cf
1 changed files with 12 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
Grid/algorithms/CoarsenedMatrix.h
View File

@ -308,6 +308,9 @@ public:
int osites=Grid()->oSites(); int osites=Grid()->oSites();
autoView(st,Stencil,AcceleratorRead);
siteVector *CBp=Stencil.CommBuf();
accelerator_for(sss, Grid()->oSites()*nbasis, Nsimd, { accelerator_for(sss, Grid()->oSites()*nbasis, Nsimd, {
int ss = sss/nbasis; int ss = sss/nbasis;
int b = sss%nbasis; int b = sss%nbasis;
@ -318,12 +321,12 @@ public:
for(int point=0;point<geom.npoint;point++){ for(int point=0;point<geom.npoint;point++){
SE=Stencil.GetEntry(ptype,point,ss); SE=st.GetEntry(ptype,point,ss);
if(SE->_is_local) { if(SE->_is_local) {
nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute); nbr = coalescedReadPermute(in_v[SE->_offset],ptype,SE->_permute);
} else { } else {
nbr = coalescedRead(Stencil.CommBuf()[SE->_offset]); nbr = coalescedRead(CBp[SE->_offset]);
} }
acceleratorSynchronise(); acceleratorSynchronise();
@ -496,8 +499,10 @@ public:
CoarseScalar InnerProd(Grid()); CoarseScalar InnerProd(Grid());
std::cout << GridLogMessage<< "CoarsenMatrix Orthog " << std::endl;
// Orthogonalise the subblocks over the basis // Orthogonalise the subblocks over the basis
blockOrthogonalise(InnerProd,Subspace.subspace); blockOrthogonalise(InnerProd,Subspace.subspace);
std::cout << GridLogMessage<< "CoarsenMatrix Orthog done " << std::endl;
// Compute the matrix elements of linop between this orthonormal // Compute the matrix elements of linop between this orthonormal
// set of vectors. // set of vectors.
@ -532,6 +537,7 @@ public:
} }
evenmask = where(mod(bcb,2)==(Integer)0,one,zero); evenmask = where(mod(bcb,2)==(Integer)0,one,zero);
oddmask = one-evenmask; oddmask = one-evenmask;
std::cout << GridLogMessage<< "CoarsenMatrix Mask " << std::endl;
assert(self_stencil!=-1); assert(self_stencil!=-1);
@ -545,6 +551,7 @@ public:
for(int p=0;p<geom.npoint;p++){ for(int p=0;p<geom.npoint;p++){
// std::cout << GridLogMessage<< "CoarsenMatrix point "<<p << std::endl;
Mphi = Mphi_p[p]; Mphi = Mphi_p[p];
int dir = geom.directions[p]; int dir = geom.directions[p];
@ -573,6 +580,7 @@ public:
} }
} }
// std::cout << GridLogMessage<< "CoarsenMatrix Diag "<<std::endl;
/////////////////////////////////////////// ///////////////////////////////////////////
// Faster alternate self coupling.. use hermiticity to save 2x // Faster alternate self coupling.. use hermiticity to save 2x
/////////////////////////////////////////// ///////////////////////////////////////////
@ -607,11 +615,11 @@ public:
if(hermitian) { if(hermitian) {
std::cout << GridLogMessage << " ForceHermitian, new code "<<std::endl; std::cout << GridLogMessage << " ForceHermitian, new code "<<std::endl;
ForceHermitian(); ForceHermitian();
std::cout << GridLogMessage << " ForceHermitian, done "<<std::endl;
} }
} }
void ForceHermitian(void) { void ForceHermitian(void) {
CoarseMatrix Diff (Grid());
for(int p=0;p<geom.npoint;p++){ for(int p=0;p<geom.npoint;p++){
int dir = geom.directions[p]; int dir = geom.directions[p];
int disp = geom.displacements[p]; int disp = geom.displacements[p];
@ -621,8 +629,7 @@ public:
int dirp = geom.directions[pp]; int dirp = geom.directions[pp];
int dispp = geom.displacements[pp]; int dispp = geom.displacements[pp];
if ( (dirp==dir) && (dispp==1) ){ if ( (dirp==dir) && (dispp==1) ){
// Diff = adj(Cshift(A[p],dir,1)) - A[pp]; std::cout << GridLogMessage<<" Replacing stencil leg "<<pp<<" with leg "<<p<<std::endl;
// std::cout << GridLogMessage<<" Replacing stencil leg "<<pp<<" with leg "<<p<< " diff "<<norm2(Diff) <<std::endl;
A[pp] = adj(Cshift(A[p],dir,1)); A[pp] = adj(Cshift(A[p],dir,1));
} }
} }