Block product implmeentation

src/problem/Formulation.cpp

@@ -1126,6 +1126,9 @@ namespace gmshddm
     template< class T_Scalar >
     gmshfem::algebra::MatrixCCS< T_Scalar > Formulation< T_Scalar >::computeMatrix()
     {
+      gmshfem::common::Timer time;
+      time.tick();
+
       const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
       const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
 
@@ -1242,9 +1245,150 @@ namespace gmshddm
 
 
       gmshfem::algebra::MatrixCCS< T_Scalar > mat(std::move(ai), std::move(aj), std::move(a));
+
+      time.tock();
+      if(mpi::getMPIRank() == 0) {
+        gmshfem::msg::info << "Computed the matrix in " << time.time() << "s." << gmshfem::msg::endl;
+      }
       return mat;
     }
 
+    template< class T_Scalar >
+    gmshfem::algebra::MatrixCCS< T_Scalar > Formulation< T_Scalar >::computeMatrixBlock()
+    {
+      gmshfem::common::Timer time;
+      time.tick();
+
+      const int outerVerbosity = gmshfem::common::Options::instance()->verbose;
+      const int innerVerbosity = (outerVerbosity != 0 ? outerVerbosity - 1 : 0);
+
+      Mat A;
+      Vec Y, Y_seq;
+      VecScatter     ctx;
+
+      // Create a distributed vector that acts as a RHS
+      PetscInt locSize = _rhs.size();
+      VecCreateMPI(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, &Y);
+      VecSetUp(Y);
+      
+      // Create a sequential full vector and a scatter object to synchronize the global vector of outputs,
+      // Needed to build the matrix
+      PetscInt globSize;
+      VecGetSize(Y, &globSize);
+
+      // Create the identity matrix
+      Mat ID;
+      {
+        MatCreate(PETSC_COMM_WORLD, &ID);
+        MatSetSizes(ID, locSize, PETSC_DECIDE, globSize, globSize);
+        MatSetType(ID, MATMPIDENSE);
+        MatSetUp(ID);
+        MatZeroEntries(ID);
+
+        PetscInt start, end;
+        MatGetOwnershipRange(ID, &start, &end);
+
+        for(PetscInt i = start; i < end; ++i) {
+          PetscScalar val = 1.0;
+          MatSetValues(ID, 1, &i, 1, &i, &val, INSERT_VALUES);
+        }
+
+        MatAssemblyBegin(ID, MAT_FINAL_ASSEMBLY);
+        MatAssemblyEnd(ID, MAT_FINAL_ASSEMBLY);
+      }
+
+      VecCreateSeq(MPI_COMM_SELF, globSize, &Y_seq);
+      VecScatterCreateToAll(Y, &ctx, &Y_seq);
+
+      if (mpi::getMPIRank() == 0)
+        gmshfem::msg::info << "Local size:" << locSize << " and global size:" << globSize << "." << gmshfem::msg::endl;
+
+
+      // Setup the A matrix
+      MatCreateShell(MPI_COMM_WORLD, locSize, locSize, globSize, globSize, this, &A);
+      MatShellSetOperation(A, MATOP_MULT, (void (*)(void)) & MatVectProductI_A< T_Scalar >);
+
+      _physicalCommutator = false;
+      _artificialCommutator = true;
+      togglePhysicalAndArtificialSourceTerms();
+
+  
+      // Compute sub matrices then disable bilinear terms (to assemble only RHS)
+      for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
+        _volume[idom]->setAttribute("ddm::physicalCommutator", _physicalCommutator);
+        _volume[idom]->setAttribute("ddm::artificialCommutator", _artificialCommutator);
+        gmshfem::common::Options::instance()->verbose = innerVerbosity;
+        _volume[idom]->removeSystem();
+        _volume[idom]->initSystem();
+        _volume[idom]->pre();
+        _volume[idom]->assemble();
+        _deactivateBilinear(idom);
+        gmshfem::common::Options::instance()->verbose = outerVerbosity;
+      }
+
+
+      Mat Yblock;
+      MatDuplicate(ID, MAT_DO_NOT_COPY_VALUES, &Yblock);
+      MatMatProductI_A< T_Scalar >(A, ID, Yblock);
+
+      std::vector< PetscInt > aiVector(globSize);
+      std::vector< PetscScalar > aVector(globSize);
+
+      std::vector< unsigned long long > aj(globSize + 1, 0);
+      std::vector< unsigned long long > ai;
+      ai.reserve(globSize);
+      std::vector< T_Scalar > a;
+      a.reserve(globSize);
+
+      for(auto i = 0ULL; i < globSize; ++i) {
+        aiVector[i] = i;
+      }
+
+      for(auto i = 0ULL; i < globSize; ++i) {
+
+        
+        if(mpi::getMPIRank() == 0) {
+
+          gmshfem::msg::info << "Column number " << i << " over " << globSize << "." << gmshfem::msg::endl;
+        }
+
+        // Scatter to rank 0
+        MatGetColumnVector(Yblock, Y, i);
+        VecScatterBegin(ctx, Y, Y_seq, INSERT_VALUES, SCATTER_FORWARD);
+        VecScatterEnd(ctx, Y, Y_seq, INSERT_VALUES, SCATTER_FORWARD);
+        VecGetValues(Y_seq, globSize, &aiVector[0], &aVector[0]);
+
+       
+        aj[i + 1] = aj[i];
+        for(auto j = 0ULL; j < aVector.size(); ++j) {
+          if(aVector[j] != 0) {
+            aj[i + 1]++;
+            ai.push_back(j);
+            a.push_back(T_Scalar(aVector[j]));
+          }
+        }
+
+        ai.reserve(aj.size() + globSize);
+        a.reserve(a.size() + globSize);
+
+      }
+      
+
+      for(auto idom = 0ULL; idom < _volume.size(); ++idom) {
+        _activateBilinear(idom);
+      }
+
+      VecScatterDestroy(&ctx);
+      gmshfem::algebra::MatrixCCS< T_Scalar > mat(std::move(ai), std::move(aj), std::move(a));
+
+      time.tock();
+      if(mpi::getMPIRank() == 0) {
+        gmshfem::msg::info << "Computed the matrix in " << time.time() << "s." << gmshfem::msg::endl;
+      }
+      return mat;
+
+    }
+
     template< class T_Scalar >
     void Formulation< T_Scalar >::setSolutionIntoInterfaceFields(const std::vector< T_Scalar > &solution)
     {
@@ -1275,6 +1419,18 @@ namespace gmshddm
       return MatVectProductImpl<T_Scalar>(A, X, Y, true);
     }
 
+    template< class T_Scalar >
+    int MatMatProductI_A(Mat A, Mat X, Mat Y)
+    {
+      return MatMatProductImpl<T_Scalar>(A, X, Y, true);
+    }
+
+    template< class T_Scalar >
+    int MatMatProductA(Mat A, Mat X, Mat Y)
+    {
+      return MatMatProductImpl<T_Scalar>(A, X, Y, false);
+    }
+
     template< class T_Scalar >
     void Formulation<T_Scalar>::_fromPetscToInterfaceFields(Vec G)
     {
@@ -1518,6 +1674,44 @@ namespace gmshddm
       PetscFunctionReturn(0);
     }
 
+    template< class T_Scalar >
+    int MatMatProductImpl(Mat A, Mat X, Mat Y, bool IA)
+    {
+
+      // EXTREMLY UNOPTIMIZED
+      PetscFunctionBegin;
+
+      PetscInt n, p, n_local, p_local;
+      MatGetSize(X, &n, &p);
+      MatGetLocalSize(X, &n_local, &p_local);
+      Vec Xcol, Ycol;
+      VecCreateMPI(MPI_COMM_WORLD, n_local, n, &Xcol);
+      VecDuplicate(Xcol, &Ycol);
+
+    for (auto j = 0; j < p; ++j) {
+      MatGetColumnVector(X, Xcol, j);
+      MatVectProductImpl<T_Scalar>(A, Xcol, Ycol, IA);
+
+      // Copy to matrix
+      PetscInt low, high, i;
+      PetscScalar val;
+
+      // Get the ownership range for the current rank
+      VecGetOwnershipRange(Ycol, &low, &high);
+
+      // Loop over local rows and set the value for the specified column
+      for (i = low; i < high; ++i) {
+          VecGetValues(Ycol, 1, &i, &val);
+          MatSetValue(Y, i, j, val, INSERT_VALUES);
+      }
+
+    }
+
+      PetscFunctionReturn(0);
+    }
+
+    
+
 
     template class Formulation< std::complex< double > >;
     template class Formulation< std::complex< float > >;

src/problem/Formulation.h

@@ -32,6 +32,12 @@ namespace gmshddm
     int MatVectProductI_A(Mat A, Vec X, Vec Y);
     template< class T_Scalar >
     int MatVectProductImpl(Mat A, Vec X, Vec Y, bool IA);
+    template< class T_Scalar >
+    int MatMatProductImpl(Mat A, Mat X, Mat Y, bool IA);
+    template< class T_Scalar >
+    int MatMatProductA(Mat A, Mat X, Mat Y);
+    template< class T_Scalar >
+    int MatMatProductI_A(Mat A, Mat X, Mat Y);
 
     template< class T_Scalar >
     class AbstractIterativeSolver;
@@ -121,6 +127,7 @@ namespace gmshddm
       gmshfem::common::Timer solve(const std::string &solver, const double tolerance = 1e-6, const int iterMax = 1000, const bool sameMatrixWithArtificialAndPhysicalSources = false, const bool skipFinalSolutionComputation = false);
       gmshfem::common::Timer solve(AbstractIterativeSolver<T_Scalar> &solver, const double tolerance = 1e-6, const int iterMax = 1000, const bool sameMatrixWithArtificialAndPhysicalSources = false, const bool skipFinalSolutionComputation = false);
       gmshfem::algebra::MatrixCCS< T_Scalar > computeMatrix();
+      gmshfem::algebra::MatrixCCS< T_Scalar > computeMatrixBlock();
 
       void setSolutionIntoInterfaceFields(const std::vector< T_Scalar > &solution);
       const gmshfem::algebra::Vector< T_Scalar > &getRHS() const;
@@ -130,7 +137,10 @@ namespace gmshddm
 
       friend int MatVectProductA< T_Scalar >(Mat A, Vec X, Vec Y);
       friend int MatVectProductI_A< T_Scalar >(Mat A, Vec X, Vec Y);
+      friend int MatMatProductI_A< T_Scalar >(Mat A, Mat X, Mat Y);
+      friend int MatMatProductA< T_Scalar >(Mat A, Mat X, Mat Y);
       friend int MatVectProductImpl< T_Scalar >(Mat A, Vec X, Vec Y, bool IA);
+      friend int MatMatProductImpl< T_Scalar >(Mat A, Mat X, Mat Y, bool IA);
 
 
       // T_Scalar, PetscScalar, T_Interger