cylqcd · sbacchio · Oct 4, 2022 · Oct 4, 2022 · Oct 4, 2022 · Oct 5, 2022
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -5,7 +5,7 @@ set(CMAKE_INCLUDE_DIRECTORIES_PROJECT_BEFORE ON)
 set(CMAKE_COLOR_MAKEFILE ON)
 set(CMAKE_CUDA_SEPARABLE_COMPILATION ON)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
-set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD 17)
 include(CheckIncludeFileCXX)
 include(CheckSymbolExists)
 
@@ -33,6 +33,16 @@ endif()
 FIND_LIBRARY(QUDA_LIB quda ${QUDA_HOME}/lib)
 include_directories(SYSTEM ${QUDA_HOME}/include)
 
+set(QUDA_PRECISION
+  "14"
+  CACHE STRING "which precisions to instantiate in QUDA (4-bit number - double, single, half, quarter)")
+set(QUDA_RECONSTRUCT
+  "7"
+  CACHE STRING "which reconstructs to instantiate in QUDA (3-bit number - 18, 13/12, 9/8)")
+add_definitions(-DQUDA_PRECISION=${QUDA_PRECISION})
+add_definitions(-DQUDA_RECONSTRUCT=${QUDA_RECONSTRUCT})
+
+
 set(QUDA_SRC "" CACHE PATH "path to QUDA src directory")
 if("${QUDA_SRC}" STREQUAL "")
   message( FATAL_ERROR "QUDA_SRC must be defined" )
@@ -149,7 +159,6 @@ set(PLEGMA_LIGHT_BARYONS FALSE CACHE BOOL "compile contractions for all light ba
 if(PLEGMA_LIGHT_BARYONS)
   add_definitions(-DPLEGMA_LIGHT_BARYONS)
 endif()
-
 set(PLEGMA_SCATTERING_CONTRACTIONS FALSE CACHE BOOL "compile contractions for pi-N scattering")
 if (PLEGMA_SCATTERING_CONTRACTIONS)
   add_definitions(-DPLEGMA_SCATTERING_CONTRACTIONS)
@@ -264,7 +273,7 @@ find_package(PythonInterp)
 if(${CMAKE_VERSION} VERSION_GREATER 3.7.99)
   find_package(CUDAWrapper)
   set(USING_CUDA_LANG_SUPPORT True)
-  set(CMAKE_CUDA_STANDARD 14)
+  set(CMAKE_CUDA_STANDARD 17)
   set(CMAKE_CUDA_STANDARD_REQUIRED True)
 else()
   set(CUDA_HOST_COMPILER "${CMAKE_CXX_COMPILER}" CACHE FILEPATH "Host side compiler used by NVCC")
@@ -415,7 +424,7 @@ add_definitions(-D__COMPUTE_CAPABILITY__=${COMP_CAP})
 
 # NVCC FLAGS independet off build type
 if(NOT USING_CUDA_LANG_SUPPORT)
-  set(CUDA_NVCC_FLAGS "-std c++11 -arch=${GPU_ARCH} -ftz=true -prec-div=false -prec-sqrt=false")
+  set(CUDA_NVCC_FLAGS "-std c++17 -arch=${GPU_ARCH} -ftz=true -prec-div=false -prec-sqrt=false")
 else()
   set(CUDA_NVCC_FLAGS "-ftz=true -prec-div=false -prec-sqrt=false")
   set(CMAKE_CUDA_FLAGS "-arch=${GPU_ARCH}" CACHE STRING "Flags used by the CUDA compiler" FORCE)
@@ -471,17 +480,17 @@ set(CMAKE_CUDA_FLAGS_DEBUG "${CUDA_NVCC_FLAGS} -g -DHOST_DEBUG -G" CACHE STRING
     FORCE )
 
 #define CXX FLAGS
-set(CMAKE_CXX_FLAGS_DEVEL  "${OpenMP_CXX_FLAGS} -O3 -Wall ${CLANG_FORCE_COLOR}" CACHE STRING
+set(CMAKE_CXX_FLAGS_DEVEL "${OpenMP_CXX_FLAGS} -std=c++17 -O3 -Wall ${CLANG_FORCE_COLOR}" CACHE STRING
 "Flags used by the C++ compiler during regular development builds.")
-set(CMAKE_CXX_FLAGS_STRICT  "${OpenMP_CXX_FLAGS} -O3 -Wall -Werror ${CLANG_NOERROR}" CACHE STRING
+set(CMAKE_CXX_FLAGS_STRICT  "${OpenMP_CXX_FLAGS} -std=c++17 -O3 -Wall -Werror ${CLANG_NOERROR}" CACHE STRING
 "Flags used by the C++ compiler during strict jenkins builds.")
-set(CMAKE_CXX_FLAGS_RELEASE "${OpenMP_CXX_FLAGS} -O3 -w" CACHE STRING
+set(CMAKE_CXX_FLAGS_RELEASE "${OpenMP_CXX_FLAGS} -std=c++17 -O3 -w" CACHE STRING
     "Flags used by the C++ compiler during release builds.")
-set(CMAKE_CXX_FLAGS_HOSTDEBUG "${OpenMP_CXX_FLAGS} -Wall -Wno-unknown-pragmas -g -fno-inline -DHOST_DEBUG ${CLANG_FORCE_COLOR}" CACHE STRING
+set(CMAKE_CXX_FLAGS_HOSTDEBUG "${OpenMP_CXX_FLAGS} -std=c++17 -Wall -Wno-unknown-pragmas -g -fno-inline -DHOST_DEBUG ${CLANG_FORCE_COLOR}" CACHE STRING
     "Flags used by the C++ compiler during host-debug builds.")
-set(CMAKE_CXX_FLAGS_DEVICEDEBUG "${OpenMP_CXX_FLAGS} -Wall -Wno-unknown-pragmas ${CLANG_FORCE_COLOR}" CACHE STRING
+set(CMAKE_CXX_FLAGS_DEVICEDEBUG "${OpenMP_CXX_FLAGS} -std=c++17 -Wall -Wno-unknown-pragmas ${CLANG_FORCE_COLOR}" CACHE STRING
     "Flags used by the C++ compiler during device-debug builds.")
-set(CMAKE_CXX_FLAGS_DEBUG "${OpenMP_CXX_FLAGS} -Wall -Wno-unknown-pragmas -g -fno-inline -DHOST_DEBUG ${CLANG_FORCE_COLOR}" CACHE STRING
+set(CMAKE_CXX_FLAGS_DEBUG "${OpenMP_CXX_FLAGS} -std=c++17 -Wall -Wno-unknown-pragmas -g -fno-inline -DHOST_DEBUG ${CLANG_FORCE_COLOR}" CACHE STRING
     "Flags used by the C++ compiler during full (host+device) debug builds.")
 
 OPTION(USE_OpenMP "Use OpenMP" ON)
@@ -555,4 +564,3 @@ add_subdirectory(tests)
 ### add PLEGMA directory
 add_subdirectory(utils)
 add_subdirectory(plegma)
-
diff --git a/include/PLEGMA_BLAS.h b/include/PLEGMA_BLAS.h
@@ -13,6 +13,8 @@
 
 #include <cublas_v2.h>
 #include <mpi.h>
+#include <PLEGMA_utils.h>
+#include <quda_api.h>
 #pragma once
 enum OPER_MATR_BLAS {NOTRANS, TRANS, DAGGER};
 namespace cBLAS{
@@ -232,8 +234,8 @@ namespace cuBLAS{
     if(trans == NOTRANS) PLEGMA_error("Use gemv without MPI comm");
     if(comm == MPI_COMM_NULL) PLEGMA_error("Communicator is NULL and cannot be used for MPI reduction");
     cuBLAS::gemv_(trans, m, n, alpha, A, x, beta, y);
-    cudaMemcpy(yHost,y,n*2*sizeof(Float),cudaMemcpyDeviceToHost);
-    checkCudaError();
+    qudaMemcpy(yHost,y,n*2*sizeof(Float),qudaMemcpyDeviceToHost);
+    checkQudaError();
     int mpiErr = MPI_Allreduce(MPI_IN_PLACE,yHost,n*2,MPI_Type(yHost),MPI_SUM,comm);
     if(mpiErr != MPI_SUCCESS) PLEGMA_error("MPI_Allreduce failed with error %d\n", mpiErr);
   }
@@ -244,8 +246,8 @@ namespace cuBLAS{
   inline void gemv(OPER_MATR_BLAS trans, int m, int n, Float alpha[2], Float* A, Float* x, Float beta[2], Float* y, MPI_Comm comm){
     Float yHost[n*2];
     cuBLAS::gemv(trans,m, n, alpha, A, x, beta, y, yHost,comm);
-    cudaMemcpy(y,yHost,sizeof(yHost),cudaMemcpyHostToDevice);
-    checkCudaError();
+    qudaMemcpy(y,yHost,sizeof(yHost),qudaMemcpyHostToDevice);
+    checkQudaError();
   }
 }
 //=================================================================//

diff --git a/include/PLEGMA_Correlator.h b/include/PLEGMA_Correlator.h
@@ -204,15 +204,68 @@ namespace plegma {
     void contractMesonsAll(PLEGMA_Propagator<Float> &prop1,
 			   PLEGMA_Propagator<Float> &prop2);
 
-
     void contractBaryons(PLEGMA_Propagator<Float> &prop1,
 			 PLEGMA_Propagator<Float> &prop2);
+
+    void contractBaryonsWall(PLEGMA_Propagator<Float> &prop1,
+			     PLEGMA_Propagator<Float> &prop2,
+			     PLEGMA_Propagator<Float> &prop3,
+			     PLEGMA_Propagator<Float> &prop4);
 
     void contractBaryonsUDSC(PLEGMA_Propagator<Float> &propUP,
 			     PLEGMA_Propagator<Float> &propDN, 
 			     PLEGMA_Propagator<Float> &propST, 
 			     PLEGMA_Propagator<Float> &propCH,
 			     bool only_st=false, bool only_ch=false);
+
+    template<typename Float2>
+    void contractTetraquarks(PLEGMA_Propagator<Float2> &propLT,
+			     PLEGMA_Propagator<Float2> &propST, 
+			     PLEGMA_Propagator<Float2> &propCH,
+                 PLEGMA_Propagator<Float2> &propBT,
+			     bool only_st=false, bool only_ch=false);
+
+    template<typename Float2>
+    void contractTetraquarksBCUD(PLEGMA_Propagator<Float2> &propLT,
+		    PLEGMA_Propagator<Float2> &propST, 
+		    PLEGMA_Propagator<Float2> &propCH, 
+            PLEGMA_Propagator<Float2> &propBT, 
+		    bool only_st=false, bool only_ch=false);
+
+    template<typename Float2>
+    void contractTetraquarkScatteringOpenIndex(
+                PLEGMA_Propagator<Float2> &prop1,
+                PLEGMA_Propagator<Float2> &prop2, 
+                std::vector<GAMMAS> gammas, int s1,
+                std::string Quarks);
+
+    template<typename Float2>
+    void contractTetraquarksStochastic(
+            PLEGMA_Propagator<Float2> &propLT1,
+		    PLEGMA_Propagator<Float2> &propLT2, 
+            PLEGMA_Propagator<Float2> &propST1, 
+            PLEGMA_Propagator<Float2> &propST2,
+		    PLEGMA_Propagator<Float2> &propBT1, 
+            PLEGMA_Propagator<Float2> &propBT2);
+
+    template<typename Float2>
+    void contractTetraquarksStochasticBCUD(PLEGMA_Propagator<Float2> &propLT1,
+		    PLEGMA_Propagator<Float2> &propLT2, 
+            PLEGMA_Propagator<Float2> &propST1, 
+            PLEGMA_Propagator<Float2> &propST2,
+		    PLEGMA_Propagator<Float2> &propCH1,
+            PLEGMA_Propagator<Float2> &propCH2,
+            PLEGMA_Propagator<Float2> &propBT1, 
+            PLEGMA_Propagator<Float2> &propBT2);
+
+
+    template<typename Float2>
+    void contractTetraquarkScatteringOpenIndexStochastic(
+            PLEGMA_Propagator<Float2> &prop1,
+            PLEGMA_Propagator<Float2> &prop2, 
+            std::vector<GAMMAS> gammas, 
+            int randInd1, int randInd2, int s1,
+            std::string Quarks);
 
     void contractNucleonThrp_local(PLEGMA_Propagator<Float> &bwdProp,
 				   PLEGMA_Propagator<Float> &fwdProp,
@@ -270,12 +323,31 @@ namespace plegma {
 					int signProps, std::vector<GAMMAS> gammas,
 					int z, std::string quark);
 
+    void contractNucleonThrp_staple(PLEGMA_Propagator<Float> &bwdProp,
+				    PLEGMA_Propagator<Float> &fwdProp,
+				    PLEGMA_Su3field<Float> &su3,
+				    int signProps, std::vector<GAMMAS> gammas, bool isZfac, int b, int l, int z);
+
     void contractNucleonThrp_qgq(PLEGMA_Propagator<Float> &bwdProp,
 				 PLEGMA_Propagator<Float> &fwdProp,
 				 PLEGMA_Su3field<Float> &su3_1,
 				 PLEGMA_Fmunu<Float> &Fmunu, std::pair<int,int> munu,
 				 PLEGMA_Su3field<Float> &su3_2,
 				 int signProps, std::vector<GAMMAS> gammas, int z2, int z1);
+    void contractMesonsFourp_ultralocal_oneendtrick(PLEGMA_Propagator<Float> &prop1,
+						    PLEGMA_Propagator<Float> &prop2,int b);
+
+    void contractMesonsFourp_ultralocal(PLEGMA_Propagator<Float> &prop1,
+					PLEGMA_Propagator<Float> &prop2,
+					PLEGMA_Propagator<Float> &prop3,
+					PLEGMA_Propagator<Float> &prop4,int b);
+
+    void contractTMDWFMesonsTrick_Zfac(PLEGMA_Propagator<Float> &prop1,
+				       PLEGMA_Su3field<float> &staple, int b, int l);
+
+    void contractTMDWFMesons_Zfac(PLEGMA_Propagator<Float> &prop1,
+				  PLEGMA_Propagator<Float> &prop2,
+                                  PLEGMA_Su3field<float> &staple, int b, int l);
 
     virtual void writeASCII(std::string filename) const;
     virtual void writeHDF5(std::string filename) const;

diff --git a/include/PLEGMA_Field.h b/include/PLEGMA_Field.h
@@ -341,6 +341,9 @@ namespace plegma {
      * @param munu: the mu,nu pair of lorentz indices of the Plq
      */
     void trPmunu(PLEGMA_Gauge<Float> &gauge, std::pair<int,int> munu);
+
+    void SU3Trace(PLEGMA_Su3field<Float> &su3field);
+
     /**
        @brief Absorbs all  elements from a 3D field and puts it at a specific global time of the 4D field
        @param PLEGMA_Field3D<Float> prop, The 3D field

diff --git a/include/PLEGMA_Hprobing.h b/include/PLEGMA_Hprobing.h
@@ -1,4 +1,6 @@
 #pragma once
+//#include <PLEGMA_utils.h>
+#include <tune_quda.h>
 
 namespace plegma {
   inline int getVecToInd(std::vector<int> x, std::vector<int> L){
@@ -100,15 +102,15 @@ namespace plegma {
       createColLattice();
       //  if(check)checkColoring();
       cudaMalloc((void**)&d_arrVc, HGC_localVolume*sizeof(int));
-      checkCudaError();
+//      checkQudaError();
       cudaMemcpy(d_arrVc, h_arrVc, HGC_localVolume*sizeof(int), cudaMemcpyHostToDevice);
-      checkCudaError();    
+//      checkQudaError();    
     }
     ~PLEGMA_Hprobing(){
       delete[] h_arrVc;
       delete[] arrlc;
       cudaFree(d_arrVc);
-      checkCudaError();
+      //checkQudaError();
     }
     int* H_arrVc() const{return h_arrVc;}
     int* D_arrVc() const{return d_arrVc;}

diff --git a/include/PLEGMA_Propagator.h b/include/PLEGMA_Propagator.h
@@ -53,10 +53,30 @@ namespace plegma {
        @return void
      **/    
     void absorb(PLEGMA_Vector3D<Float> &vec, int global_it, int nu, int c2);
+
+    /** 
+       @brief packing the sinktime slice in all other time slices
+       @param the initial propagator to be broadcasted
+       @param PLEGMA_Propagator<Float> in initial propagator to be broadcasted
+       @param int sinktimeslice  the sinktime slice to be broadcasted
+     **/ 
+    void pack_propagator_as_sink(PLEGMA_Propagator<Float> &in, int sinktimeslice, int max_source_sink_separations, bool initial);
+
+    /**
+       @brief packing the sinktime slice in all other time slices
+       @param the initial propagator to be broadcasted
+       @param PLEGMA_Propagator<Float> in initial propagator to be broadcasted
+       @param int sinktimeslice  the sinktime slice to be broadcasted
+     **/
+    void pack_propagator_from_source_to_sink(PLEGMA_Propagator<Float> &in, int sinktimeslice, int source_sink_separations, bool initial);
+
+
+
 
     void applyBoundaries_device(int t0);
     void rotateToPhysicalBase_host(int sign);
     void rotateToPhysicalBase_device(int sign);
+    void PropmulVVdag(PLEGMA_Vector<Float> &vec1,PLEGMA_Vector<Float> &vec2);
   };
 
   ///////////////////////////////

diff --git a/include/PLEGMA_ScattCorrelator.h b/include/PLEGMA_ScattCorrelator.h
@@ -101,7 +101,7 @@ namespace plegma {
      *          in most cases we do set this to false
      *
      **/
-    void V2( PLEGMA_Vector<Float> &Phi, std::vector<GAMMAS_SCATT> &Gammas, PLEGMA_Propagator<Float> &S1,  PLEGMA_Propagator<Float> &S2,bool conj_v=false );
+    void V2( PLEGMA_Vector<Float> &Phi, std::vector<GAMMAS_SCATT> &Gammas, PLEGMA_Propagator<Float> &S1,  PLEGMA_Propagator<Float> &S2,bool conj_v=true );
     /**
      *
      *  @brief performs V3 type reduction produces one spin and one color indices tensor from a fermion vector and a fermion propagator
@@ -138,7 +138,7 @@ namespace plegma {
      *
      *
      **/
-    void V4( PLEGMA_Vector<Float> &Phi, std::vector<GAMMAS_SCATT> &Gammas, PLEGMA_Propagator<Float> &S1,  PLEGMA_Propagator<Float> &S2, bool conj_v=false );
+    void V4( PLEGMA_Vector<Float> &Phi, std::vector<GAMMAS_SCATT> &Gammas, PLEGMA_Propagator<Float> &S1,  PLEGMA_Propagator<Float> &S2, bool conj_v=true );
     /**
      *
      *  @brief performs V5 type reduction produces two spin and one color indices tensor from two fermion vectors
@@ -190,7 +190,7 @@ namespace plegma {
     void T1( std::vector<GAMMAS_SCATT> &Gammas_i, std::vector<GAMMAS_SCATT> &Gammas_f, PLEGMA_Propagator<Float> &S1, PLEGMA_Propagator<Float> &S2, PLEGMA_Propagator<Float> &S3 );
     /**
      *  @brief performs T2 type reduction to compute baryon 2pt functions
-     *  T2_{alpha,beta}=\epsilon_{a,b,c}\epsilon_{l,m,n}S2^{c,l}_{alpha,beta}\Gamma_{i}_{alpha0,alpha1}S2^{b,m}_{beta0,alpha1}\Gamma_{f}_{beta0,beta1}S3^{a,n}_{beta1,alpha0} 
+     *  T2_{alpha,beta}=\epsilon_{a,b,c}\epsilon_{l,m,n}S1^{c,l}_{alpha,beta}\Gamma_{i}_{alpha0,alpha1}S2^{b,m}_{beta0,alpha1}\Gamma_{f}_{beta0,beta1}S3^{a,n}_{beta1,alpha0} 
      *  @param std::vector<GAMMAS_SCATT> &Gammas_i: list of gammas at the source for the contractions
      *  @param std::vector<GAMMAS_SCATT> &Gammas_f: list of gammas at the sink for the contractions
      *  @param PLEGMA_Propagator<Float> &S1: propagator to the sink spin from the source spin index
@@ -202,7 +202,7 @@ namespace plegma {
 
     //manipulation
     //Note that in the case of oet B, W diagram g0 refers not to gamma_i2 but to gamma_f2
-    void V3V2reduction( PLEGMA_ScattCorrelator<Float> &srcV3,PLEGMA_ScattCorrelator<Float> &srcV2, int index_abs, bool transp, int g0, bool transp_i1=false, Float* factor=NULL, bool transp_f1=false, bool oet=false );
+    void V3V2reduction( PLEGMA_ScattCorrelator<Float> &srcV3,PLEGMA_ScattCorrelator<Float> &srcV2, int index_abs, bool transp, int g0, bool transp_i1=false, Float* factor=NULL, bool transp_f1=false, bool oet=false, bool threept=false );
 
     void V5V6reduction(PLEGMA_ScattCorrelator<Float> &srcV6, 
 	                      std::shared_ptr<Float> &Phi0, 
@@ -212,7 +212,7 @@ namespace plegma {
 			      bool transp=false, bool transpgamma_i1=false, bool transgamma_f1=false, 
 			      Float *factor=NULL);
 
-    void V3V2reduction_matrix( PLEGMA_ScattCorrelator<Float> &srcV3, PLEGMA_ScattCorrelator<Float> &srcV2, int index_abs, bool transp,  int g0, bool transp_i1=false, Float* factor=NULL, bool transp_f1=false, bool oet=false ); 
+    void V3V2reduction_matrix( PLEGMA_ScattCorrelator<Float> &srcV3, PLEGMA_ScattCorrelator<Float> &srcV2, int index_abs, bool transp,  int g0, bool transp_i1=false, Float* factor=NULL, bool transp_f1=false, bool oet=false, bool threept=false ); 
 
     void V5V6reduction_matrix( PLEGMA_ScattCorrelator<Float> &srcV6, 
 	                             std::shared_ptr<Float> &Phi0, 
@@ -258,12 +258,16 @@ namespace plegma {
 
 
     //diagrams
-    void B_diagrams( PLEGMA_ScattCorrelator<Float> &srcV3, PLEGMA_ScattCorrelator<Float> &srcV2, int ig_i2, int diagram_index, bool accum=false, bool oet=false );
-    void W_diagrams( PLEGMA_ScattCorrelator<Float> &srcV3, PLEGMA_ScattCorrelator<Float> &srcV2, int ig_i2, int diagramm_index, bool accum=false );
+    void B_diagrams( PLEGMA_ScattCorrelator<Float> &srcV3, PLEGMA_ScattCorrelator<Float> &srcV2, int ig_i2, int diagram_index, bool accum=false, bool oet=false, bool threept=false );
+    void W_diagrams( PLEGMA_ScattCorrelator<Float> &srcV3, PLEGMA_ScattCorrelator<Float> &srcV2, int ig_i2, int diagramm_index, bool accum=false, bool threept=false );
     void W_diagrams_oet(PLEGMA_ScattCorrelator<Float> &srcV6, std::shared_ptr<Float> &Phi0, std::vector<PLEGMA_Vector<Float>*> &Phi_1, int input_mom_f2, int diagram_index, bool accum=false);
 
     void Z_diagrams( std::array<PLEGMA_ScattCorrelator<Float>,4> (&srcV3), std::array<PLEGMA_ScattCorrelator<Float>,4> (&srcV2),int diagramm_index, bool accum=false );
+    void Z_diagrams(std::vector<PLEGMA_ScattCorrelator<Float>*> (&srcV3), std::vector<PLEGMA_ScattCorrelator<Float>*> (&srcV2),int srcV2index, int diagramm_index, bool accum=false );
+
     void Z_diagrams_without_dilution( PLEGMA_ScattCorrelator<Float> &srcV3, PLEGMA_ScattCorrelator<Float> &srcV2, int ig_i2,  int diagramm_index, bool accum=false );
+    void Z_diagrams_without_dilution( PLEGMA_ScattCorrelator<Float> &srcV3, std::vector<PLEGMA_ScattCorrelator<Float>*> (&srcV2), int srcV2index,  int diagramm_index, bool accum=false, bool threept=false);
+
     void Z_diagrams_without_dilution_check(PLEGMA_ScattCorrelator<Float> &srcV3,
                                                 PLEGMA_ScattCorrelator<Float> &srcV2, int i_g_i2,
                                                 std::string proporder, int diagram_number, bool transp_i1, bool transp_f1, bool accum );
@@ -304,6 +308,8 @@ namespace plegma {
 
     void multiply_by_time_slice(std::shared_ptr<Float>&);
     void applyBoundaryConditions( bool antiperiodic, int n_coherent_source=1, int *attract_look_up_table=NULL );
+    void applyBoundaryConditions_3pt( bool antiperiodic, int n_coherent_source, int *attract_look_up_table, int source_sink_separation ); 
+
     void apply_phase( );
 
 

diff --git a/include/PLEGMA_Vector.h b/include/PLEGMA_Vector.h
@@ -95,6 +95,13 @@ namespace plegma {
     void rotate_uk_ch_g5g4();
     void covD(PLEGMA_Vector<Float> &vecIn, PLEGMA_Gauge<Float> &gauge, int dirOr);
     void mulGV(PLEGMA_Vector<Float> &vecIn, PLEGMA_Su3field<Float> &u);
+    void pack_fermion_to_sink(std::vector<PLEGMA_Vector<Float>*> &stochastic_vector, int sinktime);
+
+    void pack_propagator_from_source_to_sink(PLEGMA_Vector<Float> &in, int sinktimeslice, int source_sink_separation, bool initialize);
+
+    void pack_propagator_as_sink(PLEGMA_Vector<Float> &in, int sinktimeslice, int source_sink_separation, bool initialize);
+
+
   };
 
   template<typename Float> void copyToQUDA(quda::ColorSpinorField *cudaVector, Float* delem, bool isEv = false); // delem is a device pointer
@@ -135,6 +142,9 @@ namespace plegma {
      **/
     void absorb(PLEGMA_Propagator<Float> &prop, int global_it, int nu, int c2, bool  broadcast=false);
 
+    void absorb(PLEGMA_Vector<Float> &vec, int global_it, bool  broadcast=false);
+
+
     void gaussianSmearing(PLEGMA_Vector3D<Float> &vecIn, PLEGMA_Gauge3D<Float> &gauge, int nsmearGauss, Float alphaGauss) {
       this->activeTimeSlice = vecIn.activeTimeSlice;
       return ((PLEGMA_Vector<Float>*) this)->gaussianSmearing(vecIn,gauge,nsmearGauss,alphaGauss);

diff --git a/include/PLEGMA_global.h b/include/PLEGMA_global.h
@@ -43,6 +43,9 @@ using namespace std::chrono_literals;
 
 //======== Preprocessor macros =========//
 #include <global/PLEGMA_macros.hpp>
+#include <tune_quda.h>
+#include <comm_quda.h>
+using namespace quda;
 
 namespace plegma {
   //======== PLEGMA_printf, PLEGMA_error, PLEGMA_warning =========//
@@ -64,5 +67,11 @@ namespace plegma {
 #include <global/PLEGMA_global_constants.h>
 
 
+static inline void PLEGMA_memset(void *ptr, int value, size_t count){
+   return qudaMemset(ptr,value,count);
+}
+static inline void PLEGMA_memcpy(void *dst, const void *src, size_t count,qudaMemcpyKind kind){
+   return qudaMemcpy(dst,src,count,kind);
+}
 }
 using namespace plegma; // TODO: This one shouldn't be here.. But helps avoiding missing namespace.