From 7e74f7bec47818ce4873f754f5dbaaf0955d05fc Mon Sep 17 00:00:00 2001
From: Felix Erben <felix.erben@ed.ac.uk>
Date: Thu, 31 Jan 2019 11:35:05 +0000
Subject: [PATCH] tsrc != 0 now works

---
 Hadrons/Modules/MDistil/DistilVectors.hpp | 11 ++++++-----
 Hadrons/Modules/MDistil/PerambLight.hpp   | 12 ++++++------
 2 files changed, 12 insertions(+), 11 deletions(-)

diff --git a/Hadrons/Modules/MDistil/DistilVectors.hpp b/Hadrons/Modules/MDistil/DistilVectors.hpp
index caba41de..492c9755 100644
--- a/Hadrons/Modules/MDistil/DistilVectors.hpp
+++ b/Hadrons/Modules/MDistil/DistilVectors.hpp
@@ -164,25 +164,26 @@ void TDistilVectors<FImpl>::execute(void)
   bool full_tdil=(TI==Nt);
 
   int vecindex;
+  int t_inv;
   for (int inoise = 0; inoise < nnoise; inoise++) {
     for (int dk = 0; dk < LI; dk++) {
       for (int dt = 0; dt < Nt_inv; dt++) {
-        if(full_tdil) dt=tsrc; //TODO: this works for now, as longs as tsrc=0, but will crash otherwise!
         for (int ds = 0; ds < Ns; ds++) {
           vecindex = inoise + nnoise * dk + nnoise * LI * ds + nnoise *LI * Ns*dt;
           rho[vecindex] = zero;
           tmp3d_nospin = zero;
           for (int it = dt; it < Nt; it += TI){
-	    if( it >= Ntfirst && it < Ntfirst + Ntlocal ) {
+            if (full_tdil) t_inv = tsrc; else t_inv = it;
+            if( t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal ) {
               for (int ik = dk; ik < nvec; ik += LI){
                 for (int is = ds; is < Ns; is += Ns){ //at the moment, full spin dilution is enforced
-                  ExtractSliceLocal(evec3d,epack.evec[ik],0,it,3);
-                  tmp3d_nospin = evec3d * noise[inoise + nnoise*(it + Nt*(ik+nvec*is))];
+                  ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv,3);
+                  tmp3d_nospin = evec3d * noise[inoise + nnoise*(t_inv + Nt*(ik+nvec*is))];
                   //tmp3d_nospin = evec3d * noise[inoise][it][ik]()(is)(); //noises do not have to be a spin vector
                   tmp3d=zero;
                   pokeSpin(tmp3d,tmp3d_nospin,is);
                   tmp2=zero;
-                  InsertSliceLocal(tmp3d,tmp2,0,it-Ntfirst,Grid::QCD::Tdir);
+                  InsertSliceLocal(tmp3d,tmp2,0,t_inv-Ntfirst,Grid::QCD::Tdir);
                   rho[vecindex] += tmp2;
                 }
               }
diff --git a/Hadrons/Modules/MDistil/PerambLight.hpp b/Hadrons/Modules/MDistil/PerambLight.hpp
index 532c007f..7ab7372a 100644
--- a/Hadrons/Modules/MDistil/PerambLight.hpp
+++ b/Hadrons/Modules/MDistil/PerambLight.hpp
@@ -288,26 +288,26 @@ void TPerambLight<FImpl>::execute(void)
     ConjugateGradient<LatticeFermion> CG(CGPrecision,MaxIterations);
     SchurRedBlackDiagMooeeSolve<LatticeFermion> SchurSolver(CG);
 
+    int t_inv;
     for (int inoise = 0; inoise < nnoise; inoise++) {
       for (int dk = 0; dk < LI; dk++) {
         for (int dt = 0; dt < Nt_inv; dt++) {
-          if(full_tdil) dt=tsrc; //this works for now, as longs as tsrc=0, but will crash otherwise!
           for (int ds = 0; ds < Ns; ds++) {
             std::cout <<  "LapH source vector from noise " << inoise << " and dilution component (d_k,d_t,d_alpha) : (" << dk << ","<< dt << "," << ds << ")" << std::endl;
             dist_source = zero;
             tmp3d_nospin = zero;
             evec3d = zero;
             for (int it = dt; it < Nt; it += TI){
-              if( it >= Ntfirst && it < Ntfirst + Ntlocal ) {
+              if (full_tdil) t_inv = tsrc; else t_inv = it;
+              if( t_inv >= Ntfirst && t_inv < Ntfirst + Ntlocal ) {
                 for (int ik = dk; ik < nvec; ik += LI){
                   for (int is = ds; is < Ns; is += Ns){ //at the moment, full spin dilution is enforced
-                    ExtractSliceLocal(evec3d,epack.evec[ik],0,it,3);
-                    tmp3d_nospin = evec3d * noise[inoise + nnoise*(it + Nt*(ik+nvec*is))]; 
-                    //tmp3d_nospin = evec3d * noise[inoise][it][ik]()(is)(); //noises do not have to be a spin vector
+                    ExtractSliceLocal(evec3d,epack.evec[ik],0,t_inv,3);
+                    tmp3d_nospin = evec3d * noise[inoise + nnoise*(t_inv + Nt*(ik+nvec*is))]; 
                     tmp3d=zero;
                     pokeSpin(tmp3d,tmp3d_nospin,is);
                     tmp2=zero;
-                    InsertSliceLocal(tmp3d,tmp2,0,it-Ntfirst,Grid::QCD::Tdir);
+                    InsertSliceLocal(tmp3d,tmp2,0,t_inv-Ntfirst,Grid::QCD::Tdir);
                     dist_source += tmp2;
                   }
                 }