#include <base/config.h>
+#include <lac/exceptions.h>
#include <lac/trilinos_vector.h>
#include <lac/trilinos_sparse_matrix.h>
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
#include <boost/shared_ptr.hpp>
#ifdef DEAL_II_USE_TRILINOS
#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
# include <Epetra_MpiComm.h>
-# include <mpi.h>
#else
# include <Epetra_SerialComm.h>
#endif
* fact that some of the entries in the preconditioner matrix are zero
* and hence can be neglected.
*
+ * The implementation is able to distinguish between matrices from
+ * elliptic problems and convection dominated problems. We use the standard
+ * options provided by Trilinos ML for elliptic problems, except that we use a
+ * Chebyshev smoother instead of a symmetric Gauss-Seidel smoother.
+ * For most elliptic problems, Chebyshev is better than Gauss-Seidel (SSOR).
+ *
* @author Martin Kronbichler, 2008
*/
class PreconditionAMG : public Subscriptor
* Let Trilinos compute a
* multilevel hierarchy for the
* solution of a linear system
- * with the given matrix.
+ * with the given matrix. The
+ * function uses the matrix
+ * format specified in
+ * TrilinosWrappers::SparseMatrix.
+ */
+ void initialize (const SparseMatrix &matrix,
+ const std::vector<double> &null_space,
+ const unsigned int null_space_dimension,
+ const bool higher_order_elements,
+ const bool elliptic,
+ const bool output_details);
+
+ /**
+ * Let Trilinos compute a
+ * multilevel hierarchy for the
+ * solution of a linear system
+ * with the given matrix. This
+ * function takes a deal.ii matrix
+ * and copies the content into a
+ * Trilinos matrix, so the function
+ * can be considered rather
+ * inefficient.
*/
void initialize (const dealii::SparseMatrix<double> &matrix,
const std::vector<double> &null_space,
const unsigned int null_space_dimension,
const bool higher_order_elements,
const bool elliptic,
- const bool output_details,
- const double drop_tolerance = 1e-13);
+ const bool output_details);
+
+ /**
+ * This function can be used
+ * for a faster recalculation of
+ * the preconditioner construction
+ * when the matrix entries
+ * underlying the preconditioner
+ * have changed,
+ * but the matrix sparsity pattern
+ * has remained the same. What this
+ * function does is to take the
+ * already generated coarsening
+ * structure, compute the AMG
+ * prolongation and restriction
+ * according to a smoothed aggregation
+ * strategy and then builds the whole
+ * multilevel hiearchy. This function
+ * can be considerably faster in that
+ * case, since the coarsening pattern
+ * is usually the most difficult thing
+ * to do when setting up the
+ * AMG ML preconditioner.
+ */
+ void reinit ();
/**
* Multiply the source vector
* and return it in the
* destination vector.
*/
+ void vmult (Vector &dst,
+ const Vector &src) const;
+
+ /**
+ * Do the same as before, but
+ * now use deal.II vectors instead
+ * of the ones provided in the
+ * Trilinos wrapper class.
+ */
void vmult (dealii::Vector<double> &dst,
const dealii::Vector<double> &src) const;
* A copy of the deal.II matrix
* into Trilinos format.
*/
- boost::shared_ptr<Epetra_CrsMatrix> Matrix;
+ boost::shared_ptr<SparseMatrix> Matrix;
};
}
#include <base/config.h>
#include <base/subscriptor.h>
-#include <lac/compressed_sparsity_pattern.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/sparse_matrix.h>
#include <lac/exceptions.h>
#include <lac/trilinos_vector.h>
* not directly available, use one of the
* other reinit functions.
*/
- void reinit (const CompressedSparsityPattern &sparsity_pattern,
- const unsigned int n_max_entries_per_row);
+ void reinit (const SparsityPattern &sparsity_pattern,
+ const unsigned int n_max_entries_per_row);
/**
* This function initializes the
* of nonzeros from the sparsity
* pattern internally.
*/
- void reinit (const CompressedSparsityPattern &sparsity_pattern);
+ void reinit (const SparsityPattern &sparsity_pattern);
+
+ /**
+ * This function initializes the
+ * Trilinos matrix using the deal.II
+ * sparse matrix and the entries stored
+ * therein. It uses a threshold
+ * to copy only elements whose
+ * modulus is larger than the
+ * threshold (so zeros in the
+ * deal.II matrix can be filtered
+ * away).
+ */
+ void reinit (const Epetra_Map &input_map,
+ const ::dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
+ const double drop_tolerance=1e-15);
/**
* This function is similar to the
* when the matrix structure changes,
* e.g. when the grid is refined.
*/
- void reinit (const Epetra_Map &input_map,
- const CompressedSparsityPattern &sparsity_pattern);
+ void reinit (const Epetra_Map &input_map,
+ const SparsityPattern &sparsity_pattern);
/**
* This function is similar to the
* To be used e.g. for rectangular
* matrices after remeshing.
*/
- void reinit (const Epetra_Map &input_row_map,
- const Epetra_Map &input_col_map,
- const CompressedSparsityPattern &sparsity_pattern);
+ void reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
+ const SparsityPattern &sparsity_pattern);
/**
* Release all memory and return
#include <lac/trilinos_precondition_amg.h>
-#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
-
#ifdef DEAL_II_USE_TRILINOS
#include <Epetra_Map.h>
-#include <Epetra_CrsMatrix.h>
#include <Epetra_Vector.h>
#include <Teuchos_ParameterList.hpp>
#include <ml_include.h>
PreconditionAMG::~PreconditionAMG ()
{}
-
-
+
+
void
PreconditionAMG::
- initialize (const dealii::SparseMatrix<double> &matrix,
- const std::vector<double> &null_space,
- const unsigned int null_space_dimension,
- const bool elliptic,
- const bool higher_order_elements,
- const bool output_details,
- const double drop_tolerance)
+ initialize (const SparseMatrix &matrix,
+ const std::vector<double> &null_space,
+ const unsigned int null_space_dimension,
+ const bool elliptic,
+ const bool higher_order_elements,
+ const bool output_details)
{
- Assert (drop_tolerance >= 0,
- ExcMessage ("Drop tolerance must be a non-negative number."));
-
const unsigned int n_rows = matrix.m();
-
- // Init Epetra Matrix, avoid
- // storing the nonzero elements.
- {
- Map.reset (new Epetra_Map(n_rows, 0, communicator));
- std::vector<int> row_lengths (n_rows);
- for (dealii::SparseMatrix<double>::const_iterator p = matrix.begin();
- p != matrix.end(); ++p)
- if (std::abs(p->value()) > drop_tolerance)
- ++row_lengths[p->row()];
-
- Matrix.reset (new Epetra_CrsMatrix(Copy, *Map, &row_lengths[0], true));
-
- const unsigned int max_nonzero_entries
- = *std::max_element (row_lengths.begin(), row_lengths.end());
-
- std::vector<double> values(max_nonzero_entries, 0);
- std::vector<int> row_indices(max_nonzero_entries);
-
- for (unsigned int row=0; row<n_rows; ++row)
- {
- unsigned int index = 0;
- for (dealii::SparseMatrix<double>::const_iterator p = matrix.begin(row);
- p != matrix.end(row); ++p)
- if (std::abs(p->value()) > drop_tolerance)
- {
- row_indices[index] = p->column();
- values[index] = p->value();
- ++index;
- }
-
- Assert (index == static_cast<unsigned int>(row_lengths[row]),
- ExcMessage("Filtering out zeros could not "
- "be successfully finished!"));
-
- Matrix->InsertGlobalValues(row, row_lengths[row],
- &values[0], &row_indices[0]);
- }
-
- Matrix->FillComplete();
- }
-
+ if (!Matrix)
+ {
+ Matrix = boost::shared_ptr<SparseMatrix>
+ (const_cast<SparseMatrix*>(&matrix));
+ Map = boost::shared_ptr<Epetra_Map>
+ (const_cast<Epetra_Map*>(&matrix.row_map));
+ }
+
// Build the AMG preconditioner.
Teuchos::ParameterList parameter_list;
- // The implementation is able
- // to distinguish between
- // matrices from elliptic problems
- // and convection dominated
- // problems. We use the standard
- // options for elliptic problems,
- // except that we use a
- // Chebyshev smoother instead
- // of a symmetric Gauss-Seidel
- // smoother. For most elliptic
- // problems, Chebyshev is better
- // than Gauss-Seidel (SSOR).
if (elliptic)
{
ML_Epetra::SetDefaults("SA",parameter_list);
parameter_list.set("null space: dimension", int(null_space_dimension));
parameter_list.set("null space: vectors", (double *)&null_space[0]);
}
-
+
multigrid_operator = boost::shared_ptr<ML_Epetra::MultiLevelPreconditioner>
- (new ML_Epetra::MultiLevelPreconditioner(*Matrix, parameter_list, true));
+ (new ML_Epetra::MultiLevelPreconditioner(
+ *matrix.matrix, parameter_list, true));
if (output_details)
multigrid_operator->PrintUnused(0);
}
+
+
+ void
+ PreconditionAMG::
+ initialize (const dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
+ const std::vector<double> &null_space,
+ const unsigned int null_space_dimension,
+ const bool elliptic,
+ const bool higher_order_elements,
+ const bool output_details)
+ {
+ const unsigned int n_rows = deal_ii_sparse_matrix.m();
+
+ // Init Epetra Matrix, avoid
+ // storing the nonzero elements.
+
+ Map.reset (new Epetra_Map(n_rows, 0, communicator));
+
+ Matrix.reset();
+ Matrix = boost::shared_ptr<SparseMatrix> (new SparseMatrix());
+
+ Matrix->reinit (*Map, deal_ii_sparse_matrix);
+ Matrix->compress();
+
+ initialize (*Matrix, null_space, null_space_dimension, elliptic,
+ higher_order_elements, output_details);
+ }
+
+
+ void PreconditionAMG::
+ reinit ()
+ {
+ multigrid_operator->ReComputePreconditioner();
+ }
+
+
+
+ void PreconditionAMG::vmult (Vector &dst,
+ const Vector &src) const
+ {
+ const int ierr = multigrid_operator->ApplyInverse (*src.vector, *dst.vector);
+
+ Assert (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
// For the implementation of the
// <code>vmult</code> function we
// note that invoking a call of
// the Trilinos preconditioner
// requires us to use Epetra vectors
- // as well. Luckily, it is sufficient
+ // as well. It is faster
// to provide a view, i.e., feed
// Trilinos with a pointer to the
// data, so we avoid copying the
void PreconditionAMG::vmult (dealii::Vector<double> &dst,
const dealii::Vector<double> &src) const
{
+ Assert (Map->SameAs(Matrix->row_map),
+ ExcMessage("The sparse matrix given to the preconditioner uses "
+ "a map that is not compatible. Check ML preconditioner "
+ "setup."));
+ Assert (Map->SameAs(Matrix->col_map),
+ ExcMessage("The sparse matrix given to the preconditioner uses "
+ "a map that is not compatible. Check ML preconditioner "
+ "setup."));
+
Epetra_Vector LHS (View, *Map, dst.begin());
Epetra_Vector RHS (View, *Map, const_cast<double*>(src.begin()));
void
- SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern,
- const unsigned int n_max_entries_per_row)
+ SparseMatrix::reinit (const SparsityPattern &sparsity_pattern,
+ const unsigned int n_max_entries_per_row)
{
unsigned int n_rows = sparsity_pattern.n_rows();
void
- SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern)
+ SparseMatrix::reinit (const SparsityPattern &sparsity_pattern)
{
unsigned int n_rows = sparsity_pattern.n_rows();
-
+
Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
ExcDimensionMismatch (matrix->NumGlobalRows(),
sparsity_pattern.n_rows()));
-
+
// Trilinos seems to have a bug for
// rectangular matrices at this point,
// so do not check for consistent
// sparsity_pattern.n_cols()));
std::vector<int> n_entries_per_row(n_rows);
-
+
for (unsigned int row=0; row<n_rows; ++row)
n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
-
+
const unsigned int n_max_entries_per_row = *std::max_element (
&n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
void
- SparseMatrix::reinit (const Epetra_Map &input_map,
- const CompressedSparsityPattern &sparsity_pattern)
+ SparseMatrix::reinit (const Epetra_Map &input_map,
+ const SparsityPattern &sparsity_pattern)
{
+ matrix.reset();
unsigned int n_rows = sparsity_pattern.n_rows();
- matrix.reset();
- row_map = input_map;
- col_map = row_map;
-
+
Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
ExcDimensionMismatch (input_map.NumGlobalElements(),
sparsity_pattern.n_rows()));
Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
ExcDimensionMismatch (input_map.NumGlobalElements(),
sparsity_pattern.n_cols()));
-
+
+ row_map = input_map;
+ col_map = row_map;
+
std::vector<int> n_entries_per_row(n_rows);
-
+
for (unsigned int row=0; row<n_rows; ++row)
n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
-
+
matrix = std::auto_ptr<Epetra_FECrsMatrix> (new Epetra_FECrsMatrix (
Copy, row_map, &n_entries_per_row[0], true));
-
+
const unsigned int n_max_entries_per_row = *std::max_element (
&n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
void
- SparseMatrix::reinit (const Epetra_Map &input_row_map,
- const Epetra_Map &input_col_map,
- const CompressedSparsityPattern &sparsity_pattern)
+ SparseMatrix::reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
+ const SparsityPattern &sparsity_pattern)
{
+ matrix.reset();
unsigned int n_rows = sparsity_pattern.n_rows();
- matrix.reset();
- row_map = input_row_map;
- col_map = input_col_map;
-
+
Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
ExcDimensionMismatch (input_row_map.NumGlobalElements(),
sparsity_pattern.n_rows()));
Assert (input_col_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
ExcDimensionMismatch (input_col_map.NumGlobalElements(),
sparsity_pattern.n_cols()));
-
+
+ row_map = input_row_map;
+ col_map = input_col_map;
+
std::vector<int> n_entries_per_row(n_rows);
-
+
for (unsigned int row=0; row<n_rows; ++row)
n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
-
+
matrix = std::auto_ptr<Epetra_FECrsMatrix>
(new Epetra_FECrsMatrix(Copy, row_map, col_map,
&n_entries_per_row[0], true));
-
+
const unsigned int n_max_entries_per_row = *std::max_element (
&n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
+ void
+ SparseMatrix::reinit (const Epetra_Map &input_map,
+ const ::dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
+ const double drop_tolerance)
+ {
+ matrix.reset();
+
+ unsigned int n_rows = deal_ii_sparse_matrix.m();
+
+ Assert (input_map.NumGlobalElements() == (int)n_rows,
+ ExcDimensionMismatch (input_map.NumGlobalElements(),
+ n_rows));
+ Assert (input_map.NumGlobalElements() == (int)deal_ii_sparse_matrix.n(),
+ ExcDimensionMismatch (input_map.NumGlobalElements(),
+ deal_ii_sparse_matrix.n()));
+
+ row_map = input_map;
+ col_map = row_map;
+
+ std::vector<int> n_entries_per_row(n_rows);
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ n_entries_per_row[(int)row] =
+ deal_ii_sparse_matrix.get_sparsity_pattern().row_length(row);
+
+ matrix = std::auto_ptr<Epetra_FECrsMatrix>
+ (new Epetra_FECrsMatrix(Copy, row_map, col_map,
+ &n_entries_per_row[0], true));
+
+ std::vector<double> values;
+ std::vector<int> row_indices;
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ {
+ values.resize (n_entries_per_row[row],0.);
+ row_indices.resize (n_entries_per_row[row],-1);
+
+ unsigned int index = 0;
+ for (dealii::SparseMatrix<double>::const_iterator
+ p = deal_ii_sparse_matrix.begin(row);
+ p != deal_ii_sparse_matrix.end(row); ++p)
+ if (std::abs(p->value()) > drop_tolerance)
+ {
+ row_indices[index] = p->column();
+ values[index] = p->value();
+ ++index;
+ }
+ const int n_row_entries = index;
+
+ matrix->InsertGlobalValues(row, n_row_entries,
+ &values[0], &row_indices[0]);
+ }
+
+ }
+
+
+
void
SparseMatrix::clear ()
{