From 52e4efae9e3ae31c9b2884ce9b2da5d016b69471 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Sun, 11 Jan 2015 15:21:59 +0100 Subject: [PATCH] Remove SparseDirectMUMPS --- include/deal.II/lac/sparse_direct.h | 116 +--------------- source/lac/sparse_direct.cc | 203 +--------------------------- 2 files changed, 2 insertions(+), 317 deletions(-) diff --git a/include/deal.II/lac/sparse_direct.h b/include/deal.II/lac/sparse_direct.h index d485b273fe..30541d8fc2 100644 --- a/include/deal.II/lac/sparse_direct.h +++ b/include/deal.II/lac/sparse_direct.h @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2001 - 2013 by the deal.II authors +// Copyright (C) 2001 - 2015 by the deal.II authors // // This file is part of the deal.II library. // @@ -17,7 +17,6 @@ #define __deal2__sparse_direct_h - #include #include #include @@ -27,10 +26,6 @@ #include #include -#ifdef DEAL_II_WITH_MUMPS -# include -# include -#endif DEAL_II_NAMESPACE_OPEN @@ -311,115 +306,6 @@ private: std::vector control; }; - -/** - * This class provides an interface to the parallel sparse direct solver MUMPS. MUMPS is direct method based on - * a multifrontal approach, which performs a direct LU factorization. The - * matrix coming in may have either symmetric or nonsymmetric sparsity - * pattern. - * - * @note This class is useable if and only if a working installation of MUMPS exists on your system and was - * detected during configuration of deal.II. - * - *

Instantiations

- * - * There are instantiations of this class for SparseMatrix, - * SparseMatrix, BlockSparseMatrix, and - * BlockSparseMatrix. - * - * @author Markus Buerg, 2010 - */ -class SparseDirectMUMPS -{ -private: - -#ifdef DEAL_II_WITH_MUMPS - DMUMPS_STRUC_C id; -#endif // DEAL_II_WITH_MUMPS - - double *a; - std::vector rhs; - int *irn; - int *jcn; - types::global_dof_index n; - types::global_dof_index 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); - - /** - * - */ - void copy_rhs_to_mumps(const Vector &rhs); - - /** - * Flags storing whether the function initialize () has already - * been called. - */ - bool initialize_called; - -public: - /** - * Declare type for container size. - */ - typedef types::global_dof_index size_type; - - /** - * Constructor - */ - SparseDirectMUMPS (); - - /** - * Destructor - */ - ~SparseDirectMUMPS (); - - /** - * Exception - */ - DeclException0 (ExcInitializeAlreadyCalled); - - /** - * This function initializes a MUMPS instance and hands over the system's - * matrix matrix and right-hand side vector to the solver. - */ - template - 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. The right-hand side need to be - * supplied in initialize(matrix, 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 #endif // __deal2__sparse_direct_h diff --git a/source/lac/sparse_direct.cc b/source/lac/sparse_direct.cc index 85c24e9ab9..fa68f5e210 100644 --- a/source/lac/sparse_direct.cc +++ b/source/lac/sparse_direct.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2001 - 2014 by the deal.II authors +// Copyright (C) 2001 - 2015 by the deal.II authors // // This file is part of the deal.II library. // @@ -488,191 +488,6 @@ SparseDirectUMFPACK::Tvmult_add ( - -#ifdef DEAL_II_WITH_MUMPS - -SparseDirectMUMPS::SparseDirectMUMPS () - : - initialize_called (false) -{} - -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_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()); - - // Initialize MUMPS instance: - id.job = -1; - id.par = 1; - id.sym = 0; - - // Use MPI_COMM_WORLD as communicator - id.comm_fortran = -987654; - dmumps_c (&id); - - // 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 = matrix.n (); - - // number of nonzero elements in matrix - nz = matrix.n_actually_nonzero_elements (); - - // representation of the matrix - a = new double[nz]; - - // matrix indices pointing to the row and column dimensions - // respectively of the matrix representation above (a): ie. a[k] is - // the matrix element (irn[k], jcn[k]) - irn = new int[nz]; - jcn = new int[nz]; - - size_type index = 0; - - // loop over the elements of the matrix row by row, as suggested in - // the documentation of the sparse matrix iterator class - for (size_type row = 0; row < matrix.m(); ++row) - { - for (typename Matrix::const_iterator ptr = matrix.begin (row); - ptr != matrix.end (row); ++ptr) - if (std::abs (ptr->value ()) > 0.0) - { - a[index] = ptr->value (); - irn[index] = row + 1; - jcn[index] = ptr->column () + 1; - ++index; - } - } - - id.n = n; - id.nz = nz; - id.irn = irn; - id.jcn = jcn; - id.a = a; - } - - // No outputs - id.icntl[0] = -1; - id.icntl[1] = -1; - id.icntl[2] = -1; - id.icntl[3] = 0; - - // Exit by setting this flag: - initialize_called = true; -} - -template -void SparseDirectMUMPS::initialize (const Matrix &matrix, - const Vector &vector) -{ - // 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) - { - rhs.resize(n); - for (size_type i = 0; i < n; ++i) - rhs[i] = new_rhs (i); - - 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 -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) -{ - //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()); - - // and that the matrix has at least one nonzero element: - Assert (nz != 0, ExcNotInitialized()); - - // Start solver - id.job = 6; // 6 = analysis, factorization, and solve - 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()); - - // and that the matrix has at least one nonzero element: - Assert (nz != 0, ExcNotInitialized()); - - Assert(n == dst.size(), ExcMessage("Destination vector has the wrong size.")); - Assert(n == src.size(), ExcMessage("Source vector has the wrong size.")); - - // Hand over right-hand side - copy_rhs_to_mumps(src); - - // Start solver - id.job = 3; - dmumps_c (&id); - copy_solution (dst); -} - -#endif // DEAL_II_WITH_MUMPS - - // explicit instantiations for SparseMatrixUMFPACK #define InstantiateUMFPACK(MATRIX) \ template \ @@ -696,20 +511,4 @@ InstantiateUMFPACK(SparseMatrixEZ) InstantiateUMFPACK(BlockSparseMatrix) InstantiateUMFPACK(BlockSparseMatrix) -// explicit instantiations for SparseDirectMUMPS -#ifdef DEAL_II_WITH_MUMPS -#define InstantiateMUMPS(MATRIX) \ - template \ - void SparseDirectMUMPS::initialize (const MATRIX &, const Vector &);\ - template \ - void SparseDirectMUMPS::initialize (const MATRIX &); - -InstantiateMUMPS(SparseMatrix) -InstantiateMUMPS(SparseMatrix) -// InstantiateMUMPS(SparseMatrixEZ) -// InstantiateMUMPS(SparseMatrixEZ) -InstantiateMUMPS(BlockSparseMatrix) -InstantiateMUMPS(BlockSparseMatrix) -#endif - DEAL_II_NAMESPACE_CLOSE -- 2.39.5