#include <base/table.h>
#include <lac/block_matrix_base.h>
#include <lac/block_sparse_matrix.h>
-#include <lac/block_sparsity_pattern.h>
#include <lac/trilinos_sparse_matrix.h>
#include <lac/trilinos_block_vector.h>
#include <lac/exceptions.h>
DEAL_II_NAMESPACE_OPEN
+ // forward declarations
+class BlockSparsityPattern;
+class BlockCompressedSparsityPattern;
+class BlockCompressedSetSparsityPattern;
+class BlockCompressedSimpleSparsityPattern;
namespace TrilinosWrappers
* assumes that a quadratic block
* matrix is generated.
*/
+ template <typename BlockSparsityType>
void reinit (const std::vector<Epetra_Map> &input_maps,
- const BlockSparsityPattern &block_sparsity_pattern);
+ const BlockSparsityType &block_sparsity_pattern);
/**
* Resize the matrix and initialize it
* result is a block matrix for which
* all elements are stored locally.
*/
- void reinit (const BlockSparsityPattern &block_sparsity_pattern);
-
+ template <typename BlockSparsityType>
+ void reinit (const BlockSparsityType &block_sparsity_pattern);
+
/**
* This function initializes the
* Trilinos matrix using the deal.II
const ::dealii::BlockSparseMatrix<double> &deal_ii_sparse_matrix,
const double drop_tolerance=1e-13);
- /**
+ /**
+ * 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). Since no
+ * Epetra_Map is given, all the
+ * elements will be locally
+ * stored.
+ */
+ void reinit (const ::dealii::BlockSparseMatrix<double> &deal_ii_sparse_matrix,
+ const double drop_tolerance=1e-13);
+
+ /**
* This function calls the compress()
* command of all matrices after
* the assembly is
* block, so that the number of
* outer iterations is
* reduced. See the discussion
- * in the @step_22 step-22
- * tutorial program.
+ * in the @ref step_22
+ * "step-22" tutorial program.
*/
unsigned int Av_do_solve;
#include <base/config.h>
#include <base/subscriptor.h>
-#include <lac/sparsity_pattern.h>
#include <lac/sparse_matrix.h>
#include <lac/exceptions.h>
#include <lac/trilinos_vector_base.h>
DEAL_II_NAMESPACE_OPEN
+ // forward declarations
+class SparsityPattern;
+class CompressedSparsityPattern;
+class CompressedSetSparsityPattern;
+class CompressedSimpleSparsityPattern;
namespace TrilinosWrappers
{
* Epetra matrix know the
* position of nonzero entries
* according to the sparsity
- * pattern. Note that, when
- * using this function, the
- * matrix must already be
- * initialized with a suitable
- * Epetra_Map that describes
- * the distribution of the
- * matrix among the MPI
- * processes. Otherwise, an
- * error will be thrown. In a
- * parallel run, it is currently
- * necessary that each processor
- * holds the sparsity_pattern
- * structure because each
- * processor sets its rows.
+ * pattern. This function is
+ * meant for use in serial
+ * programs, where there is no
+ * need to specify how the
+ * matrix is going to be
+ * distributed among the
+ * processors. This function
+ * works in parallel, too, but
+ * it is recommended to
+ * manually specify the
+ * parallel partioning of the
+ * matrix using an
+ * Epetra_Map. When run in
+ * parallel, it is currently
+ * necessary that each
+ * processor holds the
+ * sparsity_pattern structure
+ * because each processor sets
+ * its rows.
*
- * This is a
- * collective operation that
- * needs to be called on all
- * processors in order to avoid a
- * dead lock.
+ * This is a collective
+ * operation that needs to be
+ * called on all processors in
+ * order to avoid a dead lock.
*/
- void reinit (const SparsityPattern &sparsity_pattern);
+ template<typename SparsityType>
+ void reinit (const SparsityType &sparsity_pattern);
/**
* This function is initializes
* processors in order to avoid a
* dead lock.
*/
- void reinit (const Epetra_Map &input_map,
- const SparsityPattern &sparsity_pattern);
+ template<typename SparsityType>
+ void reinit (const Epetra_Map &input_map,
+ const SparsityType &sparsity_pattern);
/**
* This function is similar to
* processors in order to avoid a
* dead lock.
*/
- void reinit (const Epetra_Map &input_row_map,
- const Epetra_Map &input_col_map,
- const SparsityPattern &sparsity_pattern);
+ template<typename SparsityType>
+ void reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
+ const SparsityType &sparsity_pattern);
/**
* This function copies the
* processors in order to avoid a
* dead lock.
*/
+ void reinit (const ::dealii::SparseMatrix<double> &dealii_sparse_matrix,
+ const double drop_tolerance=1e-13);
+
+ /**
+ * 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 with
+ * modulus larger than the
+ * threshold (so zeros in the
+ * deal.II matrix can be
+ * filtered away). In contrast
+ * to the other reinit function
+ * with deal.II sparse matrix
+ * argument, this function
+ * takes a parallel
+ * partitioning specified by
+ * the user instead of
+ * internally generating one.
+ *
+ * This is a
+ * collective operation that
+ * needs to be called on all
+ * processors in order to avoid a
+ * dead lock.
+ */
void reinit (const Epetra_Map &input_map,
const ::dealii::SparseMatrix<double> &dealii_sparse_matrix,
const double drop_tolerance=1e-13);
#include <lac/trilinos_block_sparse_matrix.h>
+#include <lac/block_sparsity_pattern.h>
#ifdef DEAL_II_USE_TRILINOS
+ template <typename BlockSparsityType>
void
BlockSparseMatrix::
reinit (const std::vector<Epetra_Map> &input_maps,
- const BlockSparsityPattern &block_sparsity_pattern)
+ const BlockSparsityType &block_sparsity_pattern)
{
Assert (input_maps.size() == block_sparsity_pattern.n_block_rows(),
ExcDimensionMismatch (input_maps.size(),
+ template <typename BlockSparsityType>
void
BlockSparseMatrix::
- reinit (const BlockSparsityPattern &block_sparsity_pattern)
+ reinit (const BlockSparsityType &block_sparsity_pattern)
{
Assert (block_sparsity_pattern.n_block_rows() ==
block_sparsity_pattern.n_block_cols(),
reinit (input_maps, block_sparsity_pattern);
}
-
+
void
+ void
+ BlockSparseMatrix::
+ reinit (const ::dealii::BlockSparseMatrix<double> &dealii_block_sparse_matrix,
+ const double drop_tolerance)
+ {
+ Assert (dealii_block_sparse_matrix.n_block_rows() ==
+ dealii_block_sparse_matrix.n_block_cols(),
+ ExcDimensionMismatch (dealii_block_sparse_matrix.n_block_rows(),
+ dealii_block_sparse_matrix.n_block_cols()));
+ Assert (dealii_block_sparse_matrix.m() ==
+ dealii_block_sparse_matrix.n(),
+ ExcDimensionMismatch (dealii_block_sparse_matrix.m(),
+ dealii_block_sparse_matrix.n()));
+
+ // produce a dummy local map and pass it
+ // off to the other function
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD);
+#else
+ Epetra_SerialComm trilinos_communicator;
+#endif
+
+ std::vector<Epetra_Map> input_maps;
+ for (unsigned int i=0; i<dealii_block_sparse_matrix.n_block_rows(); ++i)
+ input_maps.push_back (Epetra_Map(dealii_block_sparse_matrix.block(i,0).m(),
+ 0,
+ trilinos_communicator));
+
+ reinit (input_maps, dealii_block_sparse_matrix, drop_tolerance);
+ }
+
+
+
void
BlockSparseMatrix::compress()
{
return dst.l2_norm();
}
+
+
+
+
+ // -------------------- explicit instantiations -----------------------
+ //
+ template void
+ BlockSparseMatrix::reinit (const BlockSparsityPattern &);
+ template void
+ BlockSparseMatrix::reinit (const BlockCompressedSparsityPattern &);
+ template void
+ BlockSparseMatrix::reinit (const BlockCompressedSetSparsityPattern &);
+ template void
+ BlockSparseMatrix::reinit (const BlockCompressedSimpleSparsityPattern &);
+
+
+ template void
+ BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
+ const BlockSparsityPattern &);
+ template void
+ BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
+ const BlockCompressedSparsityPattern &);
+ template void
+ BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
+ const BlockCompressedSetSparsityPattern &);
+ template void
+ BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
+ const BlockCompressedSimpleSparsityPattern &);
+
}
#include <lac/trilinos_sparse_matrix.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/compressed_set_sparsity_pattern.h>
+#include <lac/compressed_simple_sparsity_pattern.h>
+
#ifdef DEAL_II_USE_TRILINOS
DEAL_II_NAMESPACE_OPEN
+ template <typename SparsityType>
void
- SparseMatrix::reinit (const SparsityPattern &sparsity_pattern)
+ SparseMatrix::reinit (const SparsityType &sparsity_pattern)
{
#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD);
+ template <typename SparsityType>
void
- SparseMatrix::reinit (const Epetra_Map &input_map,
- const SparsityPattern &sparsity_pattern)
+ SparseMatrix::reinit (const Epetra_Map &input_map,
+ const SparsityType &sparsity_pattern)
{
reinit (input_map, input_map, sparsity_pattern);
}
+ template <typename SparsityType>
void
- SparseMatrix::reinit (const Epetra_Map &input_row_map,
- const Epetra_Map &input_col_map,
- const SparsityPattern &sparsity_pattern)
+ SparseMatrix::reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
+ const SparsityType &sparsity_pattern)
{
matrix.reset();
+ // The CompressedSetSparsityPattern
+ // class stores the columns
+ // differently, so we need to
+ // explicitly rewrite this function
+ // in that case.
+ template<>
+ void
+ SparseMatrix::reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
+ const CompressedSetSparsityPattern &sparsity_pattern)
+ {
+ matrix.reset();
+
+ unsigned int n_rows = sparsity_pattern.n_rows();
+
+ if (input_row_map.Comm().MyPID() == 0)
+ {
+ 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[row] = sparsity_pattern.row_length(row);
+
+ // TODO: There seems to be problem
+ // in Epetra when a matrix is
+ // initialized with both row and
+ // column map. Maybe find something
+ // more out about this... It could
+ // be related to the bug #4123. For
+ // the moment, just ignore the
+ // column information and generate
+ // the matrix as if it were
+ // square. The call to
+ // GlobalAssemble later will set the
+ // correct values.
+ matrix = std::auto_ptr<Epetra_FECrsMatrix>
+ (new Epetra_FECrsMatrix(Copy, row_map,
+ &n_entries_per_row[0], false));
+
+
+ // TODO: As of now, assume that the
+ // sparsity pattern sits at the all
+ // processors (completely), and let
+ // each processor set its rows. Since
+ // this is wasteful, a better solution
+ // should be found in the future.
+ Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
+ ExcDimensionMismatch (matrix->NumGlobalRows(),
+ sparsity_pattern.n_rows()));
+
+ // Trilinos has a bug for rectangular
+ // matrices at this point, so do not
+ // check for consistent column numbers
+ // here.
+ //
+ // this bug is filed in the Sandia
+ // bugzilla under #4123 and should be
+ // fixed for version 9.0
+// Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(),
+// ExcDimensionMismatch (matrix->NumGlobalCols(),
+// sparsity_pattern.n_cols()));
+
+ std::vector<double> values;
+ std::vector<int> row_indices;
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ if (row_map.MyGID(row))
+ {
+ const int row_length = sparsity_pattern.row_length(row);
+ row_indices.resize (row_length, -1);
+ values.resize (row_length, 0.);
+
+ CompressedSetSparsityPattern::row_iterator col_num =
+ sparsity_pattern.row_begin (row);
+
+ for (unsigned int col = 0;
+ col_num != sparsity_pattern.row_end (row);
+ ++col_num, ++col)
+ row_indices[col] = *col_num;
+
+ matrix->InsertGlobalValues (row, row_length,
+ &values[0], &row_indices[0]);
+ }
+
+ last_action = Zero;
+
+ // In the end, the matrix needs to
+ // be compressed in order to be
+ // really ready.
+ compress();
+ }
+
+
+
void
SparseMatrix::reinit (const SparseMatrix &sparse_matrix)
{
+ void
+ SparseMatrix::reinit (const ::dealii::SparseMatrix<double> &dealii_sparse_matrix,
+ const double drop_tolerance)
+ {
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD);
+#else
+ Epetra_SerialComm trilinos_communicator;
+#endif
+
+ const Epetra_Map rows (dealii_sparse_matrix.m(),
+ 0,
+ trilinos_communicator);
+ const Epetra_Map columns (dealii_sparse_matrix.n(),
+ 0,
+ trilinos_communicator);
+ reinit (rows, columns, dealii_sparse_matrix, drop_tolerance);
+ }
+
+
+
void
SparseMatrix::reinit (const Epetra_Map &input_map,
const ::dealii::SparseMatrix<double> &dealii_sparse_matrix,
// As of now, no particularly neat
// ouput is generated in case of
// multiple processors.
- void SparseMatrix::print (std::ostream &out) const
+ void
+ SparseMatrix::print (std::ostream &out) const
{
double * values;
int * indices;
AssertThrow (out, ExcIO());
}
+
+
+
+ // explicit instantiations
+ //
+ template void
+ SparseMatrix::reinit (const SparsityPattern &);
+ template void
+ SparseMatrix::reinit (const CompressedSparsityPattern &);
+ template void
+ SparseMatrix::reinit (const CompressedSetSparsityPattern &);
+ template void
+ SparseMatrix::reinit (const CompressedSimpleSparsityPattern &);
+
+
+ template void
+ SparseMatrix::reinit (const Epetra_Map &,
+ const SparsityPattern &);
+ template void
+ SparseMatrix::reinit (const Epetra_Map &,
+ const CompressedSparsityPattern &);
+ template void
+ SparseMatrix::reinit (const Epetra_Map &,
+ const CompressedSetSparsityPattern &);
+ template void
+ SparseMatrix::reinit (const Epetra_Map &,
+ const CompressedSimpleSparsityPattern &);
+
+
+ template void
+ SparseMatrix::reinit (const Epetra_Map &,
+ const Epetra_Map &,
+ const SparsityPattern &);
+ template void
+ SparseMatrix::reinit (const Epetra_Map &,
+ const Epetra_Map &,
+ const CompressedSparsityPattern &);
+ template void
+ SparseMatrix::reinit (const Epetra_Map &,
+ const Epetra_Map &,
+ const CompressedSimpleSparsityPattern &);
+
}
DEAL_II_NAMESPACE_CLOSE