From 1e0bccbfe8364a69e97668bdb10ef2859ecee0d3 Mon Sep 17 00:00:00 2001 From: Oliver Sander <sander@igpm.rwth-aachen.de> Date: Thu, 1 Oct 2015 12:39:16 +0200 Subject: [PATCH] Make UMFPack available as a dune-solvers solver This class simply forwards all calls to the ISTL UMFPack solver. As an additional feature, this new implementation honours the ignoreNodes_ field. Currently, the price for this is an additional copy of the matrix. I'll fix that in a later release. --- dune/solvers/solvers/umfpacksolver.hh | 149 ++++++++++++++++++++++++++ 1 file changed, 149 insertions(+) create mode 100644 dune/solvers/solvers/umfpacksolver.hh diff --git a/dune/solvers/solvers/umfpacksolver.hh b/dune/solvers/solvers/umfpacksolver.hh new file mode 100644 index 00000000..90c43c3f --- /dev/null +++ b/dune/solvers/solvers/umfpacksolver.hh @@ -0,0 +1,149 @@ +// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=8 sw=4 sts=4: +#ifndef DUNE_SOLVERS_SOLVERS_UMFPACKSOLVER_HH +#define DUNE_SOLVERS_SOLVERS_UMFPACKSOLVER_HH + +/** \file + \brief A wrapper for the UMFPack sparse direct solver +*/ +#include <algorithm> +#include <vector> + +#include <dune/common/exceptions.hh> +#include <dune/common/bitsetvector.hh> +#include <dune/common/fmatrix.hh> + +#include <dune/istl/solver.hh> +#include <dune/istl/umfpack.hh> +#include <dune/istl/io.hh> + +#include <dune/solvers/common/canignore.hh> +#include <dune/solvers/solvers/solver.hh> + +namespace Dune +{ + +namespace Solvers +{ + +/** \brief Wraps the UMFPack sparse direct solver */ +template <class MatrixType, class VectorType> +class UMFPackSolver + : public Solver, public CanIgnore<Dune::BitSetVector<VectorType::block_type::dimension> > +{ + enum {blocksize = VectorType::block_type::dimension}; + + typedef typename VectorType::field_type field_type; + + typedef Dune::FieldMatrix<field_type, blocksize, blocksize> MatrixBlock; + typedef Dune::FieldVector<field_type, blocksize> VectorBlock; + +public: + + /** \brief Default constructor */ + UMFPackSolver () + : Solver(NumProc::FULL) + {} + + /** \brief Constructor for a linear problem */ + UMFPackSolver (const MatrixType& matrix, + VectorType& x, + const VectorType& rhs, + NumProc::VerbosityMode verbosity=NumProc::FULL) + : Solver(verbosity), + matrix_(&matrix), + x_(&x), + rhs_(&rhs) + {} + + void setProblem(const MatrixType& matrix, + VectorType& x, + const VectorType& rhs) + { + matrix_ = &matrix; + x_ = &x; + rhs_ = &rhs; + } + + virtual void solve() + { + // We may use the original rhs, but ISTL modifies it, so we need a non-const type here + VectorType modifiedRhs = *rhs_; + + if (not this->ignoreNodes_) + { + ///////////////////////////////////////////////////////////////// + // Solve the system + ///////////////////////////////////////////////////////////////// + Dune::InverseOperatorResult statistics; + Dune::UMFPack<MatrixType> solver(*matrix_); + solver.setOption(UMFPACK_PRL, 0); // no output at all + solver.apply(*x_, modifiedRhs, statistics); + + } + else + { + MatrixType modifiedStiffness = *matrix_; + //Dune::printmatrix(std::cout, modifiedStiffness, "mod stiffness", "--"); + + /////////////////////////////////////////// + // Modify Dirichlet rows + /////////////////////////////////////////// + + for (size_t j=0; j<modifiedStiffness.N(); j++) + for (int k=0; k<blocksize; k++) + modifiedStiffness[j][j][k][k] += 0.0; + + for (size_t j=0; j<modifiedStiffness.N(); j++) + { + auto cIt = modifiedStiffness[j].begin(); + auto cEndIt = modifiedStiffness[j].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]; + } + + //Dune::printmatrix(std::cout, modifiedStiffness, "mod stiffness", "--"); + + ///////////////////////////////////////////////////////////////// + // 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); + } + } + + /////////////////////////////////////////////////////// + // Data + /////////////////////////////////////////////////////// + + //! The quadratic term in the quadratic energy + const MatrixType* matrix_; + + //! Vector to store the solution + VectorType* x_; + + //! The linear term in the quadratic energy + const VectorType* rhs_; +}; + +} // namespace Solvers + +} // namespace Dune + +#endif -- GitLab