From: heister Date: Mon, 25 Nov 2013 16:30:50 +0000 (+0000) Subject: missing changes in r31796, clean up mumps, add asserts X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=30aba11eaa921009b8c98cec9a404af6457e268a;p=dealii-svn.git missing changes in r31796, clean up mumps, add asserts git-svn-id: https://svn.dealii.org/trunk@31797 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index ce205f7e4d..89ec0b712c 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -94,6 +94,11 @@ inconvenience this causes.
    +
  1. Fixed: Missing instantiations of SparseDirectMUMPS have been added. +
    + (Timo Heister, 2013/11/25) +
  2. +
  3. New: introduced "make test" that runs a minimal set of tests. We encourage every user to run this, especially if they run in to problems. The tests are automatically picked depending on the configuration and diff --git a/deal.II/include/deal.II/lac/sparse_direct.h b/deal.II/include/deal.II/lac/sparse_direct.h index 8ef14a7a38..c9df62991f 100644 --- a/deal.II/include/deal.II/lac/sparse_direct.h +++ b/deal.II/include/deal.II/lac/sparse_direct.h @@ -346,7 +346,7 @@ private: #endif // DEAL_II_WITH_MUMPS double *a; - double *rhs; + std::vector rhs; int *irn; int *jcn; types::global_dof_index n; @@ -364,6 +364,11 @@ private: */ void copy_solution (Vector &vector); + /** + * + */ + void copy_rhs_to_mumps(const Vector &rhs); + /** * Flags storing whether the function initialize () has already * been called. @@ -409,7 +414,8 @@ public: /** * A function in which the linear system is solved and the solution - * vector is copied into the given vector. + * vector is copied into the given vector. The right-hand side + * need to be supplied in initialize(matrix, vector); */ void solve (Vector &vector); diff --git a/deal.II/source/lac/sparse_direct.cc b/deal.II/source/lac/sparse_direct.cc index ed801cea31..4333b64fa1 100644 --- a/deal.II/source/lac/sparse_direct.cc +++ b/deal.II/source/lac/sparse_direct.cc @@ -25,6 +25,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -513,6 +514,9 @@ SparseDirectMUMPS::~SparseDirectMUMPS () template void SparseDirectMUMPS::initialize_matrix (const Matrix &matrix) { + Assert(matrix.n() == matrix.m(), ExcMessage("Matrix needs to be square.")); + + // Check we haven't been here before: Assert (initialize_called == false, ExcInitializeAlreadyCalled()); @@ -586,27 +590,35 @@ void SparseDirectMUMPS::initialize (const Matrix &matrix, // Hand over matrix and right-hand side initialize_matrix (matrix); + copy_rhs_to_mumps(vector); +} + +void SparseDirectMUMPS::copy_rhs_to_mumps(const Vector &new_rhs) +{ + Assert(n == new_rhs.size(), ExcMessage("Matrix size and rhs length must be equal.")); + if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0) { - // Object denoting a MUMPS data structure - rhs = new double[n]; - + rhs.resize(n); for (size_type i = 0; i < n; ++i) - rhs[i] = vector (i); + rhs[i] = new_rhs (i); - id.rhs = rhs; + id.rhs = &rhs[0]; } } void SparseDirectMUMPS::copy_solution (Vector &vector) { + Assert(n == vector.size(), ExcMessage("Matrix size and solution vector length must be equal.")); + Assert(n == rhs.size(), ExcMessage("Class not initialized with a rhs vector.")); + // Copy solution into the given vector if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0) { for (size_type i=0; i &vector) { + //TODO: this could be implemented similar to SparseDirectUMFPACK where + // the given vector will be used as the RHS. Sadly, there is no easy + // way to do this without breaking the interface. + // Check that the solver has been initialized by the routine above: Assert (initialize_called == true, ExcNotInitialized()); @@ -629,7 +645,7 @@ void SparseDirectMUMPS::solve (Vector &vector) Assert (nz != 0, ExcNotInitialized()); // Start solver - id.job = 6; + id.job = 6; // 6 = analysis, factorization, and solve dmumps_c (&id); copy_solution (vector); } @@ -643,17 +659,11 @@ void SparseDirectMUMPS::vmult (Vector &dst, // and that the matrix has at least one nonzero element: Assert (nz != 0, ExcNotInitialized()); - // Hand over right-hand side - if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0) - { - // Object denoting a MUMPS data structure: - rhs = new double[n]; - - for (size_type i = 0; i < n; ++i) - rhs[i] = src (i); + Assert(n == dst.size(), ExcMessage("Destination vector has the wrong size.")); + Assert(n == src.size(), ExcMessage("Source vector has the wrong size.")); - id.rhs = rhs; - } + // Hand over right-hand side + copy_rhs_to_mumps(src); // Start solver id.job = 3; @@ -691,7 +701,9 @@ InstantiateUMFPACK(BlockSparseMatrix) #ifdef DEAL_II_WITH_MUMPS #define InstantiateMUMPS(MATRIX) \ template \ - void SparseDirectMUMPS::initialize (const MATRIX &, const Vector &); + void SparseDirectMUMPS::initialize (const MATRIX &, const Vector &);\ + template \ + void SparseDirectMUMPS::initialize (const MATRIX &); InstantiateMUMPS(SparseMatrix) InstantiateMUMPS(SparseMatrix)