From 071664b7c7c205cdbd92532b328ba9a0ec718463 Mon Sep 17 00:00:00 2001
From: Boris Martin <boris.martin@uliege.be>
Date: Thu, 10 Aug 2023 11:27:03 +0200
Subject: [PATCH] Block product implmeentation
---
src/problem/Formulation.cpp | 198 +++++++++++++++++++++++++++++++++++-
src/problem/Formulation.h | 10 ++
2 files changed, 206 insertions(+), 2 deletions(-)
diff --git a/src/problem/Formulation.cpp b/src/problem/Formulation.cpp
index ba242a33..44992937 100644
--- a/src/problem/Formulation.cpp
+++ b/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);
@@ -1166,8 +1169,8 @@ namespace gmshddm
_volume[idom]->removeSystem();
_volume[idom]->initSystem();
_volume[idom]->pre();
- _volume[idom]->assemble();
- _deactivateBilinear(idom);
+ _volume[idom]->assemble();
+ _deactivateBilinear(idom);
gmshfem::common::Options::instance()->verbose = outerVerbosity;
}
@@ -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 > >;
diff --git a/src/problem/Formulation.h b/src/problem/Formulation.h
index 9b98c67e..9cd95bb6 100644
--- a/src/problem/Formulation.h
+++ b/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
--
GitLab