From 7ac6cdabadc2c41dfd148f09dcf1af5f056de7c0 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 1 Oct 2015 17:04:00 +0200
Subject: [PATCH] Do not create a temporary copy of the matrix even if the
ignoreNodes bitset is not empty
This save some memory when the matrix is large. The dune-istl UMFPack solver
(which we use internally here) actually has support for sending only submatrices
to UMFPack proper. Unfortunately, this feature is not well-documented...
Also, it only works if all dofs in a block-vector block are either completely
ignored or completely un-ignored.
---
dune/solvers/solvers/umfpacksolver.hh | 84 ++++++++++++++++-----------
1 file changed, 51 insertions(+), 33 deletions(-)
diff --git a/dune/solvers/solvers/umfpacksolver.hh b/dune/solvers/solvers/umfpacksolver.hh
index d98a7dac..1c1becd6 100644
--- a/dune/solvers/solvers/umfpacksolver.hh
+++ b/dune/solvers/solvers/umfpacksolver.hh
@@ -68,7 +68,7 @@ public:
virtual void solve()
{
// We may use the original rhs, but ISTL modifies it, so we need a non-const type here
- VectorType modifiedRhs = *rhs_;
+ VectorType mutableRhs = *rhs_;
if (not this->ignoreNodes_)
{
@@ -78,49 +78,67 @@ public:
Dune::InverseOperatorResult statistics;
Dune::UMFPack<MatrixType> solver(*matrix_);
solver.setOption(UMFPACK_PRL, 0); // no output at all
- solver.apply(*x_, modifiedRhs, statistics);
+ solver.apply(*x_, mutableRhs, statistics);
}
else
{
- MatrixType modifiedStiffness = *matrix_;
- //Dune::printmatrix(std::cout, modifiedStiffness, "mod stiffness", "--");
-
- ///////////////////////////////////////////
- // Modify Dirichlet rows
- ///////////////////////////////////////////
+ ///////////////////////////////////////////////////////////////////////////////////////////
+ // Extract the set of matrix rows that do not correspond to ignored degrees of freedom.
+ // Unfortunately, not all cases are handled by the ISTL UMFPack solver. Currently,
+ // you can only remove complete block rows. If a block is only partially ignored,
+ // the dune-istl UMFPack solver cannot do it, and the code here will throw an exception.
+ // All this can be fixed, but it needs going into the istl UMFPack code.
+ ///////////////////////////////////////////////////////////////////////////////////////////
+ std::set<std::size_t> nonIgnoreRows;
+ for (size_t i=0; i<matrix_->N(); i++)
+ {
+ auto ignoreCount = (*this->ignoreNodes_)[i].count();
+ if (ignoreCount==0)
+ nonIgnoreRows.insert(i);
+ else if (ignoreCount!=blocksize)
+ DUNE_THROW(Dune::NotImplemented, "Individual blocks must be either ignored completely, or not at all");
+ }
- for (size_t j=0; j<modifiedStiffness.N(); j++)
+ // Construct the solver
+ Dune::InverseOperatorResult statistics;
+ Dune::UMFPack<MatrixType> solver;
+ solver.setOption(UMFPACK_PRL, 0); // no output at all
+ // We eliminate all rows and columns(!) from the matrix that correspond to ignored degrees of freedom.
+ // Here is where the sparse LU decomposition is happenening.
+ solver.setSubMatrix(*matrix_,nonIgnoreRows);
+
+ // We need to modify the rhs vector by static condensation.
+ VectorType shortRhs(nonIgnoreRows.size());
+ int shortRowCount=0;
+ for (auto it = nonIgnoreRows.begin(); it!=nonIgnoreRows.end(); ++shortRowCount, ++it)
+ shortRhs[shortRowCount] = mutableRhs[*it];
+
+ shortRowCount = 0;
+ for (size_t i=0; i<matrix_->N(); i++)
{
- auto cIt = modifiedStiffness[j].begin();
- auto cEndIt = modifiedStiffness[j].end();
+ if ((*this->ignoreNodes_)[i][0])
+ continue;
+
+ auto cIt = (*matrix_)[i].begin();
+ auto cEndIt = (*matrix_)[i].end();
for (; cIt!=cEndIt; ++cIt)
- {
- for (int k=0; k<blocksize; k++)
- if ((*this->ignoreNodes_)[j][k])
- {
- (*cIt)[k] = 0;
- if (cIt.index()==j)
- (*cIt)[k][k] = 1.0;
- }
-
- }
-
- for (int k=0; k<blocksize; k++)
- if ((*this->ignoreNodes_)[j][k])
- modifiedRhs[j][k] = (*x_)[j][k];
+ if ((*this->ignoreNodes_)[cIt.index()][0])
+ cIt->mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
+
+ shortRowCount++;
}
- //Dune::printmatrix(std::cout, modifiedStiffness, "mod stiffness", "--");
+ // Solve the reduced system
+ VectorType shortX(nonIgnoreRows.size());
+ solver.apply(shortX, shortRhs, statistics);
+
+ // Blow up the solution vector
+ shortRowCount=0;
+ for (auto it = nonIgnoreRows.begin(); it!=nonIgnoreRows.end(); ++shortRowCount, ++it)
+ (*x_)[*it] = shortX[shortRowCount];
- /////////////////////////////////////////////////////////////////
- // Solve the system
- /////////////////////////////////////////////////////////////////
- Dune::InverseOperatorResult statistics;
- Dune::UMFPack<MatrixType> solver(modifiedStiffness);
- solver.setOption(UMFPACK_PRL, 0); // no output at all
- solver.apply(*x_, modifiedRhs, statistics);
}
}
--
GitLab