From 2452e29888d927c18bbdc177aafdf07ef636cdd4 Mon Sep 17 00:00:00 2001 From: Markus Buerg Date: Thu, 11 Oct 2012 09:23:54 +0000 Subject: [PATCH] Make interface similar to the one for UMFPACK. git-svn-id: https://svn.dealii.org/trunk@27061 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/lac/sparse_direct.h | 31 ++++++ deal.II/source/lac/sparse_direct.cc | 101 ++++++++++++++++---- 2 files changed, 112 insertions(+), 20 deletions(-) diff --git a/deal.II/include/deal.II/lac/sparse_direct.h b/deal.II/include/deal.II/lac/sparse_direct.h index df517a69a9..06fef86186 100644 --- a/deal.II/include/deal.II/lac/sparse_direct.h +++ b/deal.II/include/deal.II/lac/sparse_direct.h @@ -1310,6 +1310,20 @@ class SparseDirectMUMPS unsigned int n; unsigned int nz; + /** + * This function initializes a MUMPS instance + * and hands over the system's matrix + * matrix. + */ + template + void initialize_matrix (const Matrix& matrix); + + /** + * Copy the computed solution into the + * solution vector. + */ + void copy_solution (Vector& vector); + /** * Flags storing whether the function * initialize () has already been @@ -1344,12 +1358,29 @@ class SparseDirectMUMPS void initialize (const Matrix& matrix, const Vector & vector); + /** + * This function initializes a MUMPS instance + * and computes the factorization of the + * system's matrix matrix. + */ + template + void initialize (const Matrix& matrix); + /** * A function in which the linear system is * solved and the solution vector is copied * into the given vector. */ void solve (Vector& vector); + + /** + * A function in which the inverse of the + * matrix is applied to the input vector + * src and the solution is + * written into the output vector + * dst. + */ + void vmult (Vector& dst, const Vector& src); }; DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/source/lac/sparse_direct.cc b/deal.II/source/lac/sparse_direct.cc index 157ca2c383..98ec51a194 100644 --- a/deal.II/source/lac/sparse_direct.cc +++ b/deal.II/source/lac/sparse_direct.cc @@ -1957,11 +1957,21 @@ SparseDirectMUMPS::SparseDirectMUMPS () {} SparseDirectMUMPS::~SparseDirectMUMPS () -{} +{ + id.job = -2; + dmumps_c (&id); + + // Do some cleaning + if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0) + { + delete[] a; + delete[] irn; + delete[] jcn; + } +} template -void SparseDirectMUMPS::initialize (const Matrix& matrix, - const Vector & vector) +void SparseDirectMUMPS::initialize_matrix (const Matrix& matrix) { // Check we haven't been here before: Assert (initialize_called == false, ExcInitializeAlreadyCalled()); @@ -1978,11 +1988,10 @@ void SparseDirectMUMPS::initialize (const Matrix& matrix, // Hand over matrix and right-hand side if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0) { - // Objects denoting a MUMPS data structure: // // Set number of unknowns - n = vector.size (); + n = matrix.n (); // number of nonzero elements in matrix nz = matrix.n_actually_nonzero_elements (); @@ -2009,17 +2018,11 @@ void SparseDirectMUMPS::initialize (const Matrix& matrix, ++index; } - rhs = new double[n]; - - for (unsigned int i = 0; i < n; ++i) - rhs[i] = vector (i); - id.n = n; id.nz = nz; id.irn = irn; id.jcn = jcn; id.a = a; - id.rhs = rhs; } // No outputs @@ -2032,6 +2035,47 @@ void SparseDirectMUMPS::initialize (const Matrix& matrix, initialize_called = true; } +template +void SparseDirectMUMPS::initialize (const Matrix& matrix, + const Vector & vector) +{ + // Hand over matrix and right-hand side + initialize_matrix (matrix); + + if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0) + { + // Object denoting a MUMPS data structure + rhs = new double[n]; + + for (unsigned int i = 0; i < n; ++i) + rhs[i] = vector (i); + + id.rhs = rhs; + } +} + +void SparseDirectMUMPS::copy_solution (Vector& vector) +{ + // Copy solution into the given vector + if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0) + { + for (unsigned int i=0; i +void SparseDirectMUMPS::initialize (const Matrix& matrix) +{ + // Initialize MUMPS instance: + initialize_matrix (matrix); + // Start factorization + id.job = 4; + dmumps_c (&id); +} + void SparseDirectMUMPS::solve (Vector& vector) { // Check that the solver has been initialized @@ -2045,21 +2089,38 @@ void SparseDirectMUMPS::solve (Vector& vector) // Start solver id.job = 6; dmumps_c (&id); - id.job = -2; - dmumps_c (&id); + copy_solution (vector); +} + +void SparseDirectMUMPS::vmult (Vector& dst, + const Vector& src) +{ + // Check that the solver has been initialized + // by the routine above: + Assert (initialize_called == true, ExcNotInitialized()); - // Copy solution into the given vector + // 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) { - for (unsigned int i=0; i