Minor changes to residual projection: (#145)

src/libP/solvers/elliptic/ellipticResidualProjection.h

@@ -47,7 +47,7 @@ class ResidualProjection final
   void reOrthogonalize();
   bool checkOrthogonalize();
   void matvec(occa::memory& o_Ax, const dlong Ax_offset, occa::memory& o_x, const dlong x_offset);
-  void gop(dfloat*,dfloat*, const dlong);
+  void gop(dfloat*, const dlong);
   dfloat computeInnerProduct(occa::memory& o_a,
                              const dlong a_offset,
                              occa::memory& o_b,
@@ -79,8 +79,6 @@ class ResidualProjection final
   occa::kernel scaledAddwOffsetKernel;
   occa::kernel multiScaledAddwOffsetKernel;
   occa::kernel subtractedMultiScaledAddwOffsetKernel;
-  occa::kernel placeVectorKernel;
-  occa::kernel extractVectorKernel;
   occa::kernel weightedInnerProduct2Kernel;
   occa::kernel multiWeightedInnerProduct2Kernel;
   occa::kernel innerProductKernel;

src/libP/solvers/elliptic/okl/ellipticResidualProjection.okl

@@ -33,24 +33,6 @@
     if(n < N)
       x[n + offset * N] = alpha * x[n + offset * N];
 }
-@kernel void extractVector(const dlong N,
-                           @restrict const dfloat*  x,
-                           @restrict dfloat* y,
-                           const dlong column)
-{
-  for(dlong n = 0; n < N; ++n; @tile(p_threadBlockSize,@outer,@inner))
-    if(n < N)
-      y[n] = x[n + column * N];
-}
-@kernel void placeVector(const dlong N,
-                         @restrict const dfloat*  x,
-                         @restrict dfloat* y,
-                         const dlong column)
-{
-  for(dlong n = 0; n < N; ++n; @tile(p_threadBlockSize,@outer,@inner))
-    if(n < N)
-      y[n + column * N] = x[n];
-}
 @kernel void scaledAddwOffset(const dlong N,
                               const dfloat alpha,
                               @restrict dfloat*  x,

src/libP/solvers/elliptic/src/ellipticResidualProjection.c

@@ -63,18 +63,18 @@ void ResidualProjection::reOrthogonalize()
           alpha[j] = 0.5 * (computeInnerProduct(o_xx,j,o_bb,k)
                             + computeInnerProduct(o_bb,j,o_xx,k));
       }
-      gop(alpha.data() + k,work.data(),(m - k) + 1);
+      gop(alpha.data() + k,(m - k) + 1);
       o_alpha.copyFrom(alpha.data(),sizeof(dfloat)*((m-k)+1));
       subtractedMultiScaledAddwOffsetKernel(Ntotal, m, o_alpha, o_xx, one, k);
       subtractedMultiScaledAddwOffsetKernel(Ntotal, m, o_alpha, o_bb, one, k);
       dfloat normp = sqrt(alpha[k]);
       dfloat normk = 0.0;
       if(useWeightedFormulation) {
         normk = weightedInnerProduct(elliptic.mesh->ogs->o_invDegree, o_xx,k, o_bb,k);
-        gop(&normk,work.data(),1);
+        gop(&normk,1);
       } else {
         normk = computeInnerProduct(o_xx,k, o_bb,k);
-        gop(&normk,work.data(),1);
+        gop(&normk,1);
       }
       normk = sqrt(normk);
       if(normk > tol * normp) {
@@ -103,10 +103,10 @@ void ResidualProjection::matvec(occa::memory& o_Ax,
                                 const dlong x_offset)
 {
   // o_x_tmp = o_x[x_offset]
-  extractVectorKernel(Ntotal,o_x,elliptic.o_rtmp,x_offset);
+  elliptic.o_rtmp.copyFrom(o_x, Ntotal * sizeof(dfloat), 0, x_offset * Ntotal * sizeof(dfloat));
   ellipticOperator(&elliptic, elliptic.o_rtmp, elliptic.o_Ap, dfloatString);
   // o_Ax[Ax_offset] = o_Ax_tmp
-  placeVectorKernel(Ntotal,elliptic.o_Ap,o_Ax,Ax_offset);
+  o_Ax.copyFrom(elliptic.o_Ap, Ntotal * sizeof(dfloat), Ntotal * Ax_offset * sizeof(dfloat), 0);
 }
 void ResidualProjection::updateProjectionSpace()
 {
@@ -119,7 +119,7 @@ void ResidualProjection::updateProjectionSpace()
     for(int k = 0; k < m; ++k)
       alpha[k] = computeInnerProduct(o_xx,k,o_bb,m - 1);
   }
-  gop(alpha.data(),work.data(),m);
+  gop(alpha.data(),m);
   o_alpha.copyFrom(alpha.data(),sizeof(dfloat)*m);
   const dfloat norm_orig = alpha[m - 1];
   dfloat norm_new = norm_orig;
@@ -157,7 +157,7 @@ void ResidualProjection::computePreProjection(occa::memory& o_r)
     for(int k = 0; k < m; ++k)
       alpha[k] = computeInnerProduct(o_r,0,o_xx,k);
   }
-  gop(alpha.data(),work.data(),m);
+  gop(alpha.data(),m);
 
   o_alpha.copyFrom(alpha.data(), m * sizeof(dfloat));
 
@@ -233,9 +233,6 @@ ResidualProjection::ResidualProjection(elliptic_t& _elliptic,
       scalarMultiplyKernel = elliptic.mesh->device.buildKernel(fileName,
                                                                "scalarMultiply",
                                                                properties);
-      extractVectorKernel =
-        elliptic.mesh->device.buildKernel(fileName, "extractVector", properties);
-
       scaledAddwOffsetKernel = elliptic.mesh->device.buildKernel(fileName,
                                                                  "scaledAddwOffset",
                                                                  properties);
@@ -245,7 +242,6 @@ ResidualProjection::ResidualProjection(elliptic_t& _elliptic,
       subtractedMultiScaledAddwOffsetKernel = elliptic.mesh->device.buildKernel(fileName,
                                                                  "subtractedMultiScaledAddwOffset",
                                                                  properties);
-      placeVectorKernel = elliptic.mesh->device.buildKernel(fileName, "placeVector", properties);
       weightedInnerProduct2Kernel = elliptic.mesh->device.buildKernel(fileName,
                                                                       "weightedInnerProduct2",
                                                                       properties);
@@ -290,10 +286,9 @@ void ResidualProjection::preSolveProjection(occa::memory& o_r)
               << postResidualNorm << ", "
               << ratio << "\n";
 }
-void ResidualProjection::gop(dfloat* a, dfloat* work, const dlong size)
+void ResidualProjection::gop(dfloat* a, const dlong size)
 {
-  MPI_Allreduce(a, work, size, MPI_DFLOAT, MPI_SUM, elliptic.mesh->comm);
-  memcpy(a,work,size * sizeof(dfloat));
+  MPI_Allreduce(a, MPI_IN_PLACE, size, MPI_DFLOAT, MPI_SUM, elliptic.mesh->comm);
 }
 void ResidualProjection::postSolveProjection(occa::memory& o_x)
 {