From b60fa0491709540abd54417181d23869859147f6 Mon Sep 17 00:00:00 2001
From: Vivek Kale <11766050+vlkale@users.noreply.github.com>
Date: Mon, 7 Mar 2022 13:36:17 -0500
Subject: [PATCH 1/4] Update Benchmark_su3.cc

---
 benchmarks/Benchmark_su3.cc | 314 +++++++++++++++++++++++++++---------
 1 file changed, 241 insertions(+), 73 deletions(-)

diff --git a/benchmarks/Benchmark_su3.cc b/benchmarks/Benchmark_su3.cc
index c7ac5d3..e410bdb 100644
--- a/benchmarks/Benchmark_su3.cc
+++ b/benchmarks/Benchmark_su3.cc
@@ -26,6 +26,9 @@ Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
     See the full license in the file "LICENSE" in the top level distribution directory
     *************************************************************************************/
     /*  END LEGAL */
+
+ *************************************************************************************/
+    /*  END LEGAL */
 #include <Grid/GridCore.h>
 #include <Grid/GridStd.h>
 using namespace std;
@@ -33,16 +36,21 @@ using namespace Grid;
 #ifdef OMPTARGET_UVM
 #pragma omp requires unified_shared_memory
 #endif
+
+//#define TILE_UNROLL                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               
+#define UNROLL_TILE
+#define TILE_SIZE 64
+#define UNROLL_FACTOR 16
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
 #define LMAX (48)
-#define LMIN (8)
-#define LADD (8)
+#define LMIN (4)
+#define LADD (4)
 int64_t Nwarm=50;
 int64_t Nloop=1000;
- 
+
   Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::cout<<GridLogMessage << "Grid simd_layout" << simd_layout << std::endl;
   Coordinate mpi_layout  = GridDefaultMpi();
@@ -71,20 +79,19 @@ int64_t Nloop=1000;
       LatticeColourMatrix y(&Grid); random(pRNG,y);
 
       for(int64_t i=0;i<Nwarm;i++){
-	x=x*y;
+        x=x*y;
       }
       double start=usecond();
       for(int64_t i=0;i<Nloop;i++){
-	x=x*y;
+        x=x*y;
       }
       double stop=usecond();
       double time = (stop-start)/Nloop*1000.0;
-      std::cout<<"time: "<<time<<std::endl; 
+      std::cout<<"time: "<<time<<std::endl;
       double bytes=3.0*vol*Nc*Nc*sizeof(Complex);
       double footprint=2.0*vol*Nc*Nc*sizeof(Complex);
       double flops=Nc*Nc*(6.0+8.0+8.0)*vol;
       std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<footprint<<"    \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
-
     }
 
 #endif
@@ -110,77 +117,238 @@ int64_t Nloop=1000;
       auto xv=x.View();
       auto yv=y.View();
       auto zv=z.View();
-#ifdef SPOT_CHECK 
-      LatticeColourMatrix zref(&Grid); 
+#ifdef SPOT_CHECK
+      LatticeColourMatrix zref(&Grid);
       auto zref_v=zref.View();
-      
+
       std::cout<<"SPOT_CHECK mode...."<<std::endl;
-     
-      //CPU calculation
+
+      //CPU calculation                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                             
       printf("=====Beginning reference CPU calculations=====\n");
       for(int64_t s=0;s<vol;s++) {
         zref_v[s]=xv[s]*yv[s];
       }
       printf("=====End reference CPU calculations=====\n");
-     
-      
+
+
       printf("=====Beginning GPU calculations=====\n");
-      //expression template calculation if enabled
+      //expression template calculation if enabled                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  
       z=x*y;
       printf("=====End GPU calculations=====\n");
-      
+
      LatticeColourMatrix diff(&Grid);
      diff=z-zref;
      auto dv=diff.View();
-     int s=SPOT_CHECK; 
+     int s=SPOT_CHECK;
      if (s >= vol) {
-	     std::cout<<"Spot check failed; index out of bound"<<std::endl;
-     	     return -1;
+             std::cout<<"Spot check failed; index out of bound"<<std::endl;
+             return -1;
      }
      std::cout<<"s="<<s<<": CPU-GPU = "<<dv[s]<<std::endl;
      std::cout<<"s="<<s<<": GOT:" <<zv[s]<<" \n EXPECTED: "<<zref_v[s]<<std::endl;
-     
+
 
 #endif
 
-#if defined (OMPTARGET)  && defined (OMPTARGET_MAP) 
+#if defined (OMPTARGET)  && defined (OMPTARGET_MAP)
      #pragma omp target enter data map(alloc:zv._odata[ :zv.size()]) \
                                     map(to:xv._odata[ :xv.size()]) \
                                     map(to:yv._odata[ :yv.size()])
 
       for(int64_t i=0;i<Nwarm;i++){
-     #pragma omp target teams distribute parallel for thread_limit(gpu_threads)
-      for(int64_t s=0;s<vol;s++) {
-        zv[s]=xv[s]*yv[s];
-      }
-      }
 
+#ifdef UNROLL_TILE
+
+      #pragma omp target teams distribute parallel for thread_limit(gpu_threads) unroll partial(UNROLL_FACTOR) tile sizes(TILE_SIZE)
+        for(int64_t ss=0;ss<vol;ss+=TILE_SIZE*UNROLL_FACTOR) {
+
+          for(int64_t s = ss; s < ss+TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ TILE_SIZE; s < ss+ 2*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss + 2*TILE_SIZE; s < ss+3*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 3*TILE_SIZE; s < ss+4*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 4*TILE_SIZE; s < ss+5*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 5*TILE_SIZE; s < ss+6*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 6*TILE_SIZE; s < ss+7*TILE_SIZE; s++)
+		  
+          for(int64_t s = ss+ 7*TILE_SIZE; s < ss+8*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 8*TILE_SIZE; s < ss+9*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 9*TILE_SIZE; s < ss+10*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 10*TILE_SIZE; s < ss+11*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 11*TILE_SIZE; s < ss+12*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 12*TILE_SIZE; s < ss+13*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 13*TILE_SIZE; s < ss+14*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+        }
+
+#else
+
+#pragma omp target teams distribute parallel for thread_limit(gpu_threads) tile sizes(TILE_SIZE) unroll partial(UNROLL_FACTOR)
+        for(int64_t ss=0; ss<vol; ss += TILE_SIZE) {
+          for(int64_t s = ss; s < ss + TILE_SIZE; s += UNROLL_FACTOR)
+            {
+
+
+            zv[s]=xv[s]*yv[s];
+            zv[s+1] = xv[s+1]*yv[s+1];
+            zv[s+2] = xv[s+2]*yv[s+2];
+            zv[s+3] =  xv[s+3]*yv[s+3];
+            zv[s+4] = xv[s+4]*yv[s+4];
+            zv[s+5] = xv[s+5]*yv[s+5];
+            zv[s+6] = xv[s+6]*yv[s+6];
+            zv[s+7] = xv[s+7]*yv[s+7];
+            zv[s+8] = xv[s+8]*yv[s+8];
+            zv[s+9] = xv[s+9]*yv[s+9];
+            zv[s+10] =  xv[s+10]*yv[s+10];
+            zv[s+11] = xv[s+11]*yv[s+11];
+            zv[s+12] = xv[s+12]*yv[s+12];
+            zv[s+13] = xv[s+13]*yv[s+13];
+            zv[s+14] =  xv[s+14]*yv[s+14];
+            zv[s+15] =  xv[s+15]*yv[s+15];
+
+            }
+        }
+      #endif
+      }
 
       double start=usecond();
       for(int64_t i=0;i<Nloop;i++){
-      #pragma omp target teams distribute parallel for thread_limit(gpu_threads)
-      for(int64_t s=0;s<vol;s++) {
-        zv[s]=xv[s]*yv[s];
+
+#ifdef UNROLL_TILE
+
+      #pragma omp target teams distribute parallel for thread_limit(gpu_threads) unroll partial(UNROLL_FACTOR) tile sizes(TILE_SIZE)
+        for(int64_t ss=0;ss<vol;ss += TILE_SIZE*UNROLL_FACTOR) {
+
+          for(int64_t s = ss; s < ss + TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ TILE_SIZE; s < ss + 2*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+	  for(int64_t s = ss + 2*TILE_SIZE; s < ss + 3*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 3*TILE_SIZE; s < ss+4*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+		
+		
+		
+          for(int64_t s = ss+ 4*TILE_SIZE; s < ss+5*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 5*TILE_SIZE; s < ss+6*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 6*TILE_SIZE; s < ss+7*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 7*TILE_SIZE; s < ss+8*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 8*TILE_SIZE; s < ss+9*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 9*TILE_SIZE; s < ss+10*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 10*TILE_SIZE; s < ss+11*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 11*TILE_SIZE; s < ss+12*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 12*TILE_SIZE; s < ss+13*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 13*TILE_SIZE; s < ss+14*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 14*TILE_SIZE; s < ss+15*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+
+          for(int64_t s = ss+ 15*TILE_SIZE; s < ss+16*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+      }
+
+#else
+
+#pragma omp target teams distribute parallel for thread_limit(gpu_threads) tile sizes(TILE_SIZE) unroll partial(UNROLL_FACTOR)
+      for(int64_t ss=0;ss<vol;ss+=TILE_SIZE) {
+        for(int64_t s = ss; s < ss+TILE_SIZE; s+=UNROLL_FACTOR)
+          {
+            zv[s]=xv[s]*yv[s];
+            zv[s+1] = xv[s+1]*yv[s+1];
+            zv[s+2] = xv[s+2]*yv[s+2];
+            zv[s+3] =  xv[s+3]*yv[s+3];
+
+            zv[s+4]=xv[s+4]*yv[s+4];
+            zv[s+5] = xv[s+5]*yv[s+5];
+            zv[s+6] = xv[s+6]*yv[s+6];
+
+            zv[s+7]=xv[s+7]*yv[s+7];
+            zv[s+8] = xv[s+8]*yv[s+8];
+            zv[s+9] = xv[s+9]*yv[s+9];
+            zv[s+10] =  xv[s+10]*yv[s+10];
+            zv[s+11]=xv[s+11]*yv[s+11];
+            zv[s+6] = xv[s+6]*yv[s+6];
+
+            zv[s+7]=xv[s+7]*yv[s+7];
+            zv[s+8] = xv[s+8]*yv[s+8];
+            zv[s+9] = xv[s+9]*yv[s+9];
+            zv[s+10] =  xv[s+10]*yv[s+10];
+            zv[s+11]=xv[s+11]*yv[s+11];
+            zv[s+12] = xv[s+12]*yv[s+12];
+            zv[s+13] = xv[s+13]*yv[s+13];
+            zv[s+14] =  xv[s+14]*yv[s+14];
+            zv[s+15] =  xv[s+15]*yv[s+15];
+          }
+
       }
+      #endif
       }
+
       double stop=usecond();
        #pragma omp target exit data map (from:zv._odata[ :zv.size()])
        #pragma omp target exit data map (delete:yv._odata[ :yv.size()])
        #pragma omp target exit data map (delete:xv._odata[ :xv.size()])
 #else
       for(int64_t i=0;i<Nwarm;i++){
-	z=x*y;
+        z=x*y;
       }
 
 
       double start=usecond();
       for(int64_t i=0;i<Nloop;i++){
-	z=x*y;
+        z=x*y;
       }
       double stop=usecond();
 #endif
       double time = (stop-start)/Nloop*1000.0;
-      
+
       double bytes=3*vol*Nc*Nc*sizeof(Complex);
       double flops=Nc*Nc*(6+8+8)*vol;
       std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"    \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
@@ -208,16 +376,16 @@ int64_t Nloop=1000;
       LatticeColourMatrix y(&Grid); random(pRNG,y);
 
       for(int64_t i=0;i<Nwarm;i++){
-	mult(z,x,y);
+        mult(z,x,y);
       }
-      double start=usecond();
+  double start=usecond();
       for(int64_t i=0;i<Nloop;i++){
-	mult(z,x,y);
-	//	mac(z,x,y);
+        mult(z,x,y);
+        //      mac(z,x,y);                                                                                                                                                                                                                                                                                                                                                                                                                                                                  
       }
       double stop=usecond();
       double time = (stop-start)/Nloop*1000.0;
-      
+
       double bytes=3*vol*Nc*Nc*sizeof(Complex);
       double flops=Nc*Nc*(6+8+8)*vol;
       std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"    \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
@@ -243,22 +411,22 @@ int64_t Nloop=1000;
       LatticeColourMatrix y(&Grid); random(pRNG,y);
 
       for(int64_t i=0;i<Nwarm;i++){
-	mac(z,x,y);
+        mac(z,x,y);
       }
       double start=usecond();
       for(int64_t i=0;i<Nloop;i++){
-	mac(z,x,y);
+        mac(z,x,y);
       }
       double stop=usecond();
       double time = (stop-start)/Nloop*1000.0;
-      
+
       double bytes=4*vol*Nc*Nc*sizeof(Complex);
       double flops=Nc*Nc*(8+8+8)*vol;
       std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
 
     }
 
-#endif 
+#endif
 
 #if 0
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
@@ -267,6 +435,7 @@ int64_t Nloop=1000;
   std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
   std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
 
+  for(int lat=LMIN;lat<=LMAX;lat+=LADD){
   for(int lat=LMIN;lat<=LMAX;lat+=LADD){
 
       std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
@@ -280,17 +449,16 @@ int64_t Nloop=1000;
       LatticeColourMatrix y(&Grid); random(pRNG,y);
 
       for(int mu=0;mu<4;mu++){
-	      double start=usecond();
-	      for(int64_t i=0;i<Nloop;i++){
-	        z = PeriodicBC::CovShiftForward(x,mu,y);
-	    }
-	    double stop=usecond();
-	    double time = (stop-start)/Nloop*1000.0;
-	
-	
-	    double bytes=3*vol*Nc*Nc*sizeof(Complex);
-	    double flops=Nc*Nc*(6+8+8)*vol;
-	    std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
+              double start=usecond();
+              for(int64_t i=0;i<Nloop;i++){
+                z = PeriodicBC::CovShiftForward(x,mu,y);
+            }
+            double stop=usecond();
+            double time = (stop-start)/Nloop*1000.0;
+
+            double bytes=3*vol*Nc*Nc*sizeof(Complex);
+            double flops=Nc*Nc*(6+8+8)*vol;
+            std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
       }
   }
 
@@ -315,28 +483,28 @@ int64_t Nloop=1000;
       LatticeColourMatrix tmp(&Grid);
 
       for(int mu=0;mu<4;mu++){
-	double tshift=0;
-	double tmult =0;
-
-	double start=usecond();
-	for(int64_t i=0;i<Nloop;i++){
-	  tshift-=usecond();
-	  tmp = Cshift(y,mu,-1);
-	  tshift+=usecond();
-	  tmult-=usecond();
-	  z   = x*tmp;
-	  tmult+=usecond();
+        double tshift=0;
+        double tmult =0;
+
+        double start=usecond();
+        for(int64_t i=0;i<Nloop;i++){
+          tshift-=usecond();
+          tmp = Cshift(y,mu,-1);
+          tshift+=usecond();
+          tmult-=usecond();
+          z   = x*tmp;
+          tmult+=usecond();
 	}
-	double stop=usecond();
-	double time = (stop-start)/Nloop;
-	tshift = tshift/Nloop;
-	tmult  = tmult /Nloop;
-	
-	double bytes=3*vol*Nc*Nc*sizeof(Complex);
-	double flops=Nc*Nc*(6+8+8)*vol;
-	std::cout<<GridLogMessage<<std::setprecision(3) << "total us "<<time<<" shift "<<tshift <<" mult "<<tmult<<std::endl;
-	time = time * 1000; // convert to NS for GB/s
-	std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
+        double stop=usecond();
+        double time = (stop-start)/Nloop;
+        tshift = tshift/Nloop;
+        tmult  = tmult /Nloop;
+
+        double bytes=3*vol*Nc*Nc*sizeof(Complex);
+        double flops=Nc*Nc*(6+8+8)*vol;
+        std::cout<<GridLogMessage<<std::setprecision(3) << "total us "<<time<<" shift "<<tshift <<" mult "<<tmult<<std::endl;
+        time = time * 1000; // convert to NS for GB/s                                                                                                                                                                                                                                                                                                                                                                                                                                        
+        std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
       }
     }
 #endif

From 5ac5358565774a2255580b6691cca34dc386ba14 Mon Sep 17 00:00:00 2001
From: Vivek Kale <11766050+vlkale@users.noreply.github.com>
Date: Mon, 7 Mar 2022 13:38:44 -0500
Subject: [PATCH 2/4] Update Benchmark_su3.cc

---
 benchmarks/Benchmark_su3.cc | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/benchmarks/Benchmark_su3.cc b/benchmarks/Benchmark_su3.cc
index e410bdb..f334b9e 100644
--- a/benchmarks/Benchmark_su3.cc
+++ b/benchmarks/Benchmark_su3.cc
@@ -202,7 +202,12 @@ int64_t Nloop=1000;
 
           for(int64_t s = ss+ 13*TILE_SIZE; s < ss+14*TILE_SIZE; s++)
             zv[s]=xv[s]*yv[s];
-
+		
+        for(int64_t s = ss+ 14*TILE_SIZE; s < ss+15*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+         for(int64_t s = ss+ 15*TILE_SIZE; s < ss+16*TILE_SIZE; s++)
+            zv[s]=xv[s]*yv[s];
+		
         }
 
 #else
@@ -211,8 +216,6 @@ int64_t Nloop=1000;
         for(int64_t ss=0; ss<vol; ss += TILE_SIZE) {
           for(int64_t s = ss; s < ss + TILE_SIZE; s += UNROLL_FACTOR)
             {
-
-
             zv[s]=xv[s]*yv[s];
             zv[s+1] = xv[s+1]*yv[s+1];
             zv[s+2] = xv[s+2]*yv[s+2];

From 6f2ba873b57e408beaad471228b977b7fb58adc2 Mon Sep 17 00:00:00 2001
From: Vivek Kale <11766050+vlkale@users.noreply.github.com>
Date: Fri, 18 Mar 2022 16:05:24 -0400
Subject: [PATCH 3/4] Update Benchmark_su3.cc

---
 benchmarks/Benchmark_su3.cc | 156 ++----------------------------------
 1 file changed, 7 insertions(+), 149 deletions(-)

diff --git a/benchmarks/Benchmark_su3.cc b/benchmarks/Benchmark_su3.cc
index f334b9e..6731c39 100644
--- a/benchmarks/Benchmark_su3.cc
+++ b/benchmarks/Benchmark_su3.cc
@@ -160,80 +160,15 @@ int64_t Nloop=1000;
 #ifdef UNROLL_TILE
 
       #pragma omp target teams distribute parallel for thread_limit(gpu_threads) unroll partial(UNROLL_FACTOR) tile sizes(TILE_SIZE)
-        for(int64_t ss=0;ss<vol;ss+=TILE_SIZE*UNROLL_FACTOR) {
-
-          for(int64_t s = ss; s < ss+TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ TILE_SIZE; s < ss+ 2*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss + 2*TILE_SIZE; s < ss+3*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 3*TILE_SIZE; s < ss+4*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 4*TILE_SIZE; s < ss+5*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 5*TILE_SIZE; s < ss+6*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 6*TILE_SIZE; s < ss+7*TILE_SIZE; s++)
-		  
-          for(int64_t s = ss+ 7*TILE_SIZE; s < ss+8*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 8*TILE_SIZE; s < ss+9*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 9*TILE_SIZE; s < ss+10*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 10*TILE_SIZE; s < ss+11*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 11*TILE_SIZE; s < ss+12*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 12*TILE_SIZE; s < ss+13*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 13*TILE_SIZE; s < ss+14*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-		
-        for(int64_t s = ss+ 14*TILE_SIZE; s < ss+15*TILE_SIZE; s++)
+        for(int64_t s=0;s<vol;s++) { 
             zv[s]=xv[s]*yv[s];
-         for(int64_t s = ss+ 15*TILE_SIZE; s < ss+16*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-		
         }
 
 #else
 
 #pragma omp target teams distribute parallel for thread_limit(gpu_threads) tile sizes(TILE_SIZE) unroll partial(UNROLL_FACTOR)
-        for(int64_t ss=0; ss<vol; ss += TILE_SIZE) {
-          for(int64_t s = ss; s < ss + TILE_SIZE; s += UNROLL_FACTOR)
-            {
-            zv[s]=xv[s]*yv[s];
-            zv[s+1] = xv[s+1]*yv[s+1];
-            zv[s+2] = xv[s+2]*yv[s+2];
-            zv[s+3] =  xv[s+3]*yv[s+3];
-            zv[s+4] = xv[s+4]*yv[s+4];
-            zv[s+5] = xv[s+5]*yv[s+5];
-            zv[s+6] = xv[s+6]*yv[s+6];
-            zv[s+7] = xv[s+7]*yv[s+7];
-            zv[s+8] = xv[s+8]*yv[s+8];
-            zv[s+9] = xv[s+9]*yv[s+9];
-            zv[s+10] =  xv[s+10]*yv[s+10];
-            zv[s+11] = xv[s+11]*yv[s+11];
-            zv[s+12] = xv[s+12]*yv[s+12];
-            zv[s+13] = xv[s+13]*yv[s+13];
-            zv[s+14] =  xv[s+14]*yv[s+14];
-            zv[s+15] =  xv[s+15]*yv[s+15];
-
-            }
+        for(int64_t s=0; s<vol; s++) {
+		zv[s]=xv[s]*yv[s]; 
         }
       #endif
       }
@@ -244,95 +179,18 @@ int64_t Nloop=1000;
 #ifdef UNROLL_TILE
 
       #pragma omp target teams distribute parallel for thread_limit(gpu_threads) unroll partial(UNROLL_FACTOR) tile sizes(TILE_SIZE)
-        for(int64_t ss=0;ss<vol;ss += TILE_SIZE*UNROLL_FACTOR) {
-
-          for(int64_t s = ss; s < ss + TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ TILE_SIZE; s < ss + 2*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-	  for(int64_t s = ss + 2*TILE_SIZE; s < ss + 3*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 3*TILE_SIZE; s < ss+4*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-		
-		
-		
-          for(int64_t s = ss+ 4*TILE_SIZE; s < ss+5*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 5*TILE_SIZE; s < ss+6*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 6*TILE_SIZE; s < ss+7*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 7*TILE_SIZE; s < ss+8*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 8*TILE_SIZE; s < ss+9*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 9*TILE_SIZE; s < ss+10*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 10*TILE_SIZE; s < ss+11*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 11*TILE_SIZE; s < ss+12*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 12*TILE_SIZE; s < ss+13*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 13*TILE_SIZE; s < ss+14*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 14*TILE_SIZE; s < ss+15*TILE_SIZE; s++)
-            zv[s]=xv[s]*yv[s];
-
-          for(int64_t s = ss+ 15*TILE_SIZE; s < ss+16*TILE_SIZE; s++)
+        for(int64_t s=0;s<vol;s ++) {
             zv[s]=xv[s]*yv[s];
       }
 
 #else
 
 #pragma omp target teams distribute parallel for thread_limit(gpu_threads) tile sizes(TILE_SIZE) unroll partial(UNROLL_FACTOR)
-      for(int64_t ss=0;ss<vol;ss+=TILE_SIZE) {
-        for(int64_t s = ss; s < ss+TILE_SIZE; s+=UNROLL_FACTOR)
-          {
-            zv[s]=xv[s]*yv[s];
-            zv[s+1] = xv[s+1]*yv[s+1];
-            zv[s+2] = xv[s+2]*yv[s+2];
-            zv[s+3] =  xv[s+3]*yv[s+3];
-
-            zv[s+4]=xv[s+4]*yv[s+4];
-            zv[s+5] = xv[s+5]*yv[s+5];
-            zv[s+6] = xv[s+6]*yv[s+6];
-
-            zv[s+7]=xv[s+7]*yv[s+7];
-            zv[s+8] = xv[s+8]*yv[s+8];
-            zv[s+9] = xv[s+9]*yv[s+9];
-            zv[s+10] =  xv[s+10]*yv[s+10];
-            zv[s+11]=xv[s+11]*yv[s+11];
-            zv[s+6] = xv[s+6]*yv[s+6];
-
-            zv[s+7]=xv[s+7]*yv[s+7];
-            zv[s+8] = xv[s+8]*yv[s+8];
-            zv[s+9] = xv[s+9]*yv[s+9];
-            zv[s+10] =  xv[s+10]*yv[s+10];
-            zv[s+11]=xv[s+11]*yv[s+11];
-            zv[s+12] = xv[s+12]*yv[s+12];
-            zv[s+13] = xv[s+13]*yv[s+13];
-            zv[s+14] =  xv[s+14]*yv[s+14];
-            zv[s+15] =  xv[s+15]*yv[s+15];
-          }
-
+      for(int64_t s=0;s<vol;s++) {
+      	zv[s]=xv[s]*yv[s];
       }
       #endif
-      }
+      }  
 
       double stop=usecond();
        #pragma omp target exit data map (from:zv._odata[ :zv.size()])

From e48def8a7affc1c57938896322274dcd36ffd434 Mon Sep 17 00:00:00 2001
From: Vivek Kale <11766050+vlkale@users.noreply.github.com>
Date: Tue, 26 Apr 2022 23:43:21 -0400
Subject: [PATCH 4/4] Adding OpenMP multi-GPU parallelization for OpenMP
 offload

Adding OpenMP multi-GPU parallelization for OpenMP offload in the version without function overloading of the star operator.
---
 benchmarks/Benchmark_su3.cc | 52 +++++++++++++++++++++++++------------
 1 file changed, 35 insertions(+), 17 deletions(-)

diff --git a/benchmarks/Benchmark_su3.cc b/benchmarks/Benchmark_su3.cc
index 6731c39..18aa1e3 100644
--- a/benchmarks/Benchmark_su3.cc
+++ b/benchmarks/Benchmark_su3.cc
@@ -155,42 +155,60 @@ int64_t Nloop=1000;
                                     map(to:xv._odata[ :xv.size()]) \
                                     map(to:yv._odata[ :yv.size()])
 
-      for(int64_t i=0;i<Nwarm;i++){
-
+  for(int64_t i=0;i<Nwarm;i++) {  
+     int chunkSz = 4;  // TODO: this should be global 
+     int taskNum = 0;  // TODO: this should be a global 
+     int chosenDev = 0; 
+     int numDevs = omp_get_num_devices(); 
+     int numTasks =  (ceil) (vol/chunkSz); 
+	      
+    #pragma omp parallel for schedule(dynamic) private(chosenDev)
+	     for (taskNum = 0; taskNum < numTasks; taskNum++) 
+	      {
+		      chosenDev = omp_get_thread_num()%numDevs; 
+		     
 #ifdef UNROLL_TILE
-
-      #pragma omp target teams distribute parallel for thread_limit(gpu_threads) unroll partial(UNROLL_FACTOR) tile sizes(TILE_SIZE)
-        for(int64_t s=0;s<vol;s++) { 
+      #pragma omp target teams distribute parallel for simd thread_limit(gpu_threads) unroll partial(UNROLL_FACTOR) tile sizes(TILE_SIZE) is_device_ptr(taskNum, numTasks) device(chosenDev)
+	for(int64_t s=(taskNum)*(vol/numTasks);s< (taskNum+1)*(vol/numTasks) ;s++) {
             zv[s]=xv[s]*yv[s];
-        }
-
+        } 
 #else
 
-#pragma omp target teams distribute parallel for thread_limit(gpu_threads) tile sizes(TILE_SIZE) unroll partial(UNROLL_FACTOR)
-        for(int64_t s=0; s<vol; s++) {
+#pragma omp target teams distribute parallel for simd thread_limit(gpu_threads) tile sizes(TILE_SIZE) unroll partial(UNROLL_FACTOR) is_device_ptr(taskNum, numTasks) device(chosenDev)
+	for(int64_t s=(taskNum)*(vol/numTasks); s< (taskNum+1)*(vol/numTasks) ;s++) { 
 		zv[s]=xv[s]*yv[s]; 
         }
       #endif
       }
 
+     int chunkSz = 4;  // TODO: this should be global  
+     int taskNum = 0; 
+     int chosenDev = 0; 
+     int numDevs = omp_get_num_devices(); 
+     int numTasks =  (ceil) (vol/chunkSz); 
+	  
       double start=usecond();
       for(int64_t i=0;i<Nloop;i++){
-
+	      
+	      #pragma omp parallel for schedule(dynamic) private(chosenDev)
+	      for (taskNum = 0; taskNum < numTasks; taskNum++) 
+	      {
+		      chosenDev = omp_get_thread_num()%numDevs; 
 #ifdef UNROLL_TILE
-
-      #pragma omp target teams distribute parallel for thread_limit(gpu_threads) unroll partial(UNROLL_FACTOR) tile sizes(TILE_SIZE)
-        for(int64_t s=0;s<vol;s ++) {
+      #pragma omp target teams distribute parallel for simd thread_limit(gpu_threads)  unroll partial(UNROLL_FACTOR) tile sizes(TILE_SIZE) is_device_ptr(taskNum, numTasks) device(chosenDev)
+      for(int64_t s=(taskNum)*(vol/numTasks);s< (taskNum+1)*(vol/numTasks) ;s++) {
             zv[s]=xv[s]*yv[s];
       }
 
 #else
-
-#pragma omp target teams distribute parallel for thread_limit(gpu_threads) tile sizes(TILE_SIZE) unroll partial(UNROLL_FACTOR)
-      for(int64_t s=0;s<vol;s++) {
+	      
+#pragma omp target teams distribute parallel for simd thread_limit(gpu_threads) tile sizes(TILE_SIZE) unroll partial(UNROLL_FACTOR) is_device_ptr(taskNum, numTasks) device(chosenDev)
+      for(int64_t s=(taskNum)*(vol/numTasks);s< (taskNum+1)*(vol/numTasks) ;s++) {
       	zv[s]=xv[s]*yv[s];
       }
       #endif
-      }  
+     }
+   }  
 
       double stop=usecond();
        #pragma omp target exit data map (from:zv._odata[ :zv.size()])