CompressedSetSparsityPattern &,
const bool,
const Table<2,bool> &) const;
-
template void ConstraintMatrix::
add_entries_local_to_global<CompressedSimpleSparsityPattern> (const std::vector<unsigned int> &,
CompressedSimpleSparsityPattern &,
const bool,
const Table<2,bool> &) const;
+#ifdef DEAL_II_USE_TRILINOS
+template void ConstraintMatrix::
+add_entries_local_to_global<TrilinosWrappers::SparsityPattern>
+ (const std::vector<unsigned int> &,
+ TrilinosWrappers::SparsityPattern &,
+ const bool,
+ const Table<2,bool> &) const;
+#endif
template void ConstraintMatrix::
add_entries_local_to_global<BlockSparsityPattern> (const std::vector<unsigned int> &,
- BlockSparsityPattern &,
+ BlockSparsityPattern &,
const bool,
const Table<2,bool> &) const;
template void ConstraintMatrix::
const Table<2,bool> &) const;
template void ConstraintMatrix::
add_entries_local_to_global<BlockCompressedSetSparsityPattern> (const std::vector<unsigned int> &,
- BlockCompressedSetSparsityPattern &,
+ BlockCompressedSetSparsityPattern &,
const bool,
const Table<2,bool> &) const;
-
template void ConstraintMatrix::
add_entries_local_to_global<BlockCompressedSimpleSparsityPattern> (const std::vector<unsigned int> &,
- BlockCompressedSimpleSparsityPattern &,
+ BlockCompressedSimpleSparsityPattern &,
const bool,
const Table<2,bool> &) const;
+#ifdef DEAL_II_USE_TRILINOS
+template void ConstraintMatrix::
+add_entries_local_to_global<TrilinosWrappers::BlockSparsityPattern>
+ (const std::vector<unsigned int> &,
+ TrilinosWrappers::BlockSparsityPattern &,
+ const bool,
+ const Table<2,bool> &) const;
+#endif
DEAL_II_NAMESPACE_CLOSE
#include <base/quadrature_lib.h>
#include <base/table.h>
#include <base/template_constraints.h>
+#include <base/utilities.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
#include <grid/intergrid_map.h>
#include <lac/compressed_sparsity_pattern.h>
#include <lac/compressed_set_sparsity_pattern.h>
#include <lac/compressed_simple_sparsity_pattern.h>
+#include <lac/trilinos_sparsity_pattern.h>
#include <lac/block_sparsity_pattern.h>
#include <lac/vector.h>
#include <numerics/vectors.h>
dofs_on_this_cell.reserve (max_dofs_per_cell(dof));
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
+
+ // In case we work with a distributed
+ // sparsity pattern of Trilinos type, we
+ // only have to do the work if the
+ // current cell is owned by the calling
+ // processor. Otherwise, just continue.
for (; cell!=endc; ++cell)
+#ifdef DEAL_II_USE_TRILINOS
+ if ((types_are_equal<SparsityPattern,TrilinosWrappers::SparsityPattern>::value
+ ||
+ types_are_equal<SparsityPattern,TrilinosWrappers::BlockSparsityPattern>::value)
+ &&
+ cell->subdomain_id() !=
+ Utilities::Trilinos::get_this_mpi_process(Utilities::Trilinos::comm_world()))
+ continue;
+ else
+
+#endif
{
const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
dofs_on_this_cell.resize (dofs_per_cell);
std::vector<unsigned int> dofs_on_this_cell(fe_collection.max_dofs_per_cell());
typename DH::active_cell_iterator cell = dof.begin_active(),
endc = dof.end();
- for (; cell!=endc; ++cell)
+
+ // In case we work with a distributed
+ // sparsity pattern of Trilinos type, we
+ // only have to do the work if the
+ // current cell is owned by the calling
+ // processor. Otherwise, just continue.
+ for (; cell!=endc; ++cell)
+#ifdef DEAL_II_USE_TRILINOS
+ if ((types_are_equal<SparsityPattern,TrilinosWrappers::SparsityPattern>::value
+ ||
+ types_are_equal<SparsityPattern,TrilinosWrappers::BlockSparsityPattern>::value)
+ &&
+ cell->subdomain_id() !=
+ Utilities::Trilinos::get_this_mpi_process(Utilities::Trilinos::comm_world()))
+ continue;
+ else
+
+#endif
{
const unsigned int fe_index
= cell->active_fe_index();
CompressedSimpleSparsityPattern &sparsity,
const ConstraintMatrix &,
const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
+ TrilinosWrappers::SparsityPattern>
+(const DoFHandler<deal_II_dimension> &dof,
+ TrilinosWrappers::SparsityPattern &sparsity,
+ const ConstraintMatrix &,
+ const bool);
+#endif
+
template void
DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
BlockSparsityPattern>
BlockCompressedSimpleSparsityPattern &sparsity,
const ConstraintMatrix &,
const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
+ TrilinosWrappers::BlockSparsityPattern>
+(const DoFHandler<deal_II_dimension> &dof,
+ TrilinosWrappers::BlockSparsityPattern &sparsity,
+ const ConstraintMatrix &,
+ const bool);
+#endif
template void
DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
SparsityPattern &sparsity,
const ConstraintMatrix &,
const bool);
-
-
template void
DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
CompressedSparsityPattern>
CompressedSparsityPattern &sparsity,
const ConstraintMatrix &,
const bool);
-
template void
DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
CompressedSetSparsityPattern>
CompressedSetSparsityPattern &sparsity,
const ConstraintMatrix &,
const bool);
-
template void
DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
CompressedSimpleSparsityPattern>
CompressedSimpleSparsityPattern &sparsity,
const ConstraintMatrix &,
const bool);
-
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
+ TrilinosWrappers::SparsityPattern>
+(const hp::DoFHandler<deal_II_dimension> &dof,
+ TrilinosWrappers::SparsityPattern &sparsity,
+ const ConstraintMatrix &,
+ const bool);
+#endif
template void
DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
BlockCompressedSimpleSparsityPattern &sparsity,
const ConstraintMatrix &,
const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
+ TrilinosWrappers::BlockSparsityPattern>
+(const hp::DoFHandler<deal_II_dimension> &dof,
+ TrilinosWrappers::BlockSparsityPattern &sparsity,
+ const ConstraintMatrix &,
+ const bool);
+#endif
+
template void
CompressedSimpleSparsityPattern&,
const ConstraintMatrix &,
const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
+ TrilinosWrappers::SparsityPattern>
+(const DoFHandler<deal_II_dimension>&,
+ const Table<2,Coupling>&,
+ TrilinosWrappers::SparsityPattern&,
+ const ConstraintMatrix &,
+ const bool);
+#endif
+
template void
DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
BlockSparsityPattern>
BlockCompressedSimpleSparsityPattern&,
const ConstraintMatrix &,
const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
+ TrilinosWrappers::BlockSparsityPattern>
+(const DoFHandler<deal_II_dimension>&,
+ const Table<2,Coupling>&,
+ TrilinosWrappers::BlockSparsityPattern&,
+ const ConstraintMatrix &,
+ const bool);
+#endif
template void
DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
CompressedSimpleSparsityPattern&,
const ConstraintMatrix &,
const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
+ TrilinosWrappers::SparsityPattern>
+(const hp::DoFHandler<deal_II_dimension>&,
+ const Table<2,Coupling>&,
+ TrilinosWrappers::SparsityPattern&,
+ const ConstraintMatrix &,
+ const bool);
+#endif
+
template void
DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
BlockSparsityPattern>
BlockCompressedSimpleSparsityPattern&,
const ConstraintMatrix &,
const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
+ TrilinosWrappers::BlockSparsityPattern>
+(const hp::DoFHandler<deal_II_dimension>&,
+ const Table<2,Coupling>&,
+ TrilinosWrappers::BlockSparsityPattern&,
+ const ConstraintMatrix &,
+ const bool);
+#endif
+
template void
#include <base/subscriptor.h>
#include <base/smartpointer.h>
#include <lac/sparsity_pattern.h>
+#include <lac/trilinos_sparsity_pattern.h>
#include <lac/compressed_sparsity_pattern.h>
#include <lac/compressed_set_sparsity_pattern.h>
#include <lac/compressed_simple_sparsity_pattern.h>
class BlockCompressedSparsityPattern;
class BlockCompressedSimpleSparsityPattern;
class BlockCompressedSetSparsityPattern;
+#ifdef DEAL_II_USE_TRILINOS
+namespace TrilinosWrappers
+{
+ class BlockSparsityPattern;
+}
+#endif
/*! @addtogroup Sparsity
*@{
* There are several, exchangeable variations of this class, see @ref Sparsity,
* section 'Dynamic block sparsity patterns' for more information.
*
- * This class is used in @ref step_22 "step-22".
- *
* @author Wolfgang Bangerth, 2007
*/
class BlockCompressedSetSparsityPattern : public BlockSparsityPatternBase<CompressedSetSparsityPattern>
* There are several, exchangeable variations of this class, see @ref Sparsity,
* section 'Dynamic block sparsity patterns' for more information.
*
+ * This class is used in @ref step_22 "step-22" and @ref step_31 "step-31".
+ *
* @author Timo Heister, 2008
*/
class BlockCompressedSimpleSparsityPattern : public BlockSparsityPatternBase<CompressedSimpleSparsityPattern>
* you assign them sizes.
*/
BlockCompressedSimpleSparsityPattern (const unsigned int n_rows,
- const unsigned int n_columns);
+ const unsigned int n_columns);
/**
* Initialize the pattern with
+
+#ifdef DEAL_II_USE_TRILINOS
+
+
+/**
+ * This class extends the base class to implement an array of Trilinos
+ * sparsity patterns that can be used to initialize Trilinos block sparse
+ * matrices that can be distributed among different processors. It is used
+ * in the same way as the BlockSparsityPattern except that it builds upon
+ * the TrilinosWrappers::SparsityPattern instead of the
+ * dealii::SparsityPattern. See the documentation of the
+ * BlockSparsityPattern for examples.
+ *
+ * This class is has properties of the "dynamic" type of @ref Sparsity (in
+ * the sense that it can extend the memory if too little elements were
+ * allocated), but otherwise is more like the basic deal.II SparsityPattern
+ * (in the sense that the method compress() needs to be called before the
+ * pattern can be used).
+ *
+ * This class is used in @ref step_32 "step-32".
+ *
+ * @author Martin Kronbichler, 2008
+ */
+namespace TrilinosWrappers
+{
+ class BlockSparsityPattern :
+ public dealii::BlockSparsityPatternBase<SparsityPattern>
+ {
+ public:
+
+ /**
+ * Initialize the matrix empty,
+ * that is with no memory
+ * allocated. This is useful if
+ * you want such objects as
+ * member variables in other
+ * classes. You can make the
+ * structure usable by calling
+ * the reinit() function.
+ */
+ BlockSparsityPattern ();
+
+ /**
+ * Initialize the matrix with the
+ * given number of block rows and
+ * columns. The blocks themselves
+ * are still empty, and you have
+ * to call collect_sizes() after
+ * you assign them sizes.
+ */
+ BlockSparsityPattern (const unsigned int n_rows,
+ const unsigned int n_columns);
+
+ /**
+ * Initialize the pattern with
+ * two BlockIndices for the block
+ * structures of matrix rows and
+ * columns. This function is
+ * equivalent to calling the
+ * previous constructor with the
+ * length of the two index vector
+ * and then entering the index
+ * values.
+ */
+ BlockSparsityPattern (const std::vector<unsigned int>& row_block_sizes,
+ const std::vector<unsigned int>& col_block_sizes);
+
+ /**
+ * Initialize the pattern with an array
+ * Epetra_Map that specifies both rows
+ * and columns of the matrix (so the
+ * final matrix will be a square
+ * matrix), where the Epetra_Map
+ * specifies the parallel distribution
+ * of the degrees of freedom on the
+ * individual block. This function is
+ * equivalent to calling the second
+ * constructor with the length of the
+ * mapping vector and then entering the
+ * index values.
+ */
+ BlockSparsityPattern (const std::vector<Epetra_Map>& input_maps);
+
+ /**
+ * Resize the matrix to a tensor
+ * product of matrices with
+ * dimensions defined by the
+ * arguments.
+ *
+ * The matrix will have as many
+ * block rows and columns as
+ * there are entries in the two
+ * arguments. The block at
+ * position (<i>i,j</i>) will
+ * have the dimensions
+ * <tt>row_block_sizes[i]</tt>
+ * times <tt>col_block_sizes[j]</tt>.
+ */
+ void reinit (const std::vector< unsigned int > &row_block_sizes,
+ const std::vector< unsigned int > &col_block_sizes);
+
+ /**
+ * Resize the matrix to a square tensor
+ * product of matrices with parallel
+ * distribution according to the
+ * specifications in the array of
+ * Epetra_Maps.
+ */
+ void reinit (const std::vector<Epetra_Map>& input_maps);
+
+
+ /**
+ * Allow the use of the reinit
+ * functions of the base class as
+ * well.
+ */
+ using BlockSparsityPatternBase<SparsityPattern>::reinit;
+ };
+}
+
+#endif
+
+
/*@}*/
/*---------------------- Template functions -----------------------------------*/
DEAL_II_NAMESPACE_OPEN
- // forward declarations
-class SparsityPattern;
-class CompressedSparsityPattern;
-class CompressedSetSparsityPattern;
-class CompressedSimpleSparsityPattern;
namespace TrilinosWrappers
{
// forward declarations
class VectorBase;
class SparseMatrix;
+ class SparsityPattern;
namespace MatrixIterators
{
* over the elements of Trilinos matrices. The implementation of this
* class is similar to the one for PETSc matrices.
*
- * Note that Trilinos does not give any guarantees as to the order of
- * elements within each row. Note also that accessing the elements of
- * a full matrix surprisingly only shows the nonzero elements of the
- * matrix, not all elements.
+ * Note that Trilinos stores the elements within each row in ascending
+ * order. This is opposed to the deal.II sparse matrix style where the
+ * diagonal element (if it exists) is stored before all other values, and
+ * the PETSc sparse matrices, where one can't guarantee a certain order of
+ * the elements.
*
* @ingroup TrilinosWrappers
* @author Martin Kronbichler, Wolfgang Bangerth, 2008
* access, a const matrix
* pointer is sufficient.
*/
- Accessor (const SparseMatrix *matrix,
- const unsigned int row,
- const unsigned int index);
+ Accessor (const SparseMatrix *matrix,
+ const unsigned int row,
+ const unsigned int index);
/**
* Row number of the element
* matrix for the given row and
* the index within it.
*/
- const_iterator (const SparseMatrix *matrix,
+ const_iterator (const SparseMatrix *matrix,
const unsigned int row,
const unsigned int index);
SparseMatrix ();
/**
- * Constructor using an
- * Epetra_Map and a maximum
- * number of nonzero matrix
- * entries. Note that this
- * number does not need to be
- * exact, and it is even
- * allowed that the actual
- * matrix structure has more
- * nonzero entries than
- * specified in the
- * constructor. However it is
- * still advantageous to
- * provide good estimates here
- * since this will considerably
- * increase the performance of
- * the matrix.
+ * Constructor using an Epetra_Map
+ * and a maximum number of nonzero
+ * matrix entries. Note that this
+ * number does not need to be exact,
+ * and it is even allowed that the
+ * actual matrix structure has more
+ * nonzero entries than specified in
+ * the constructor. However it is
+ * still advantageous to provide good
+ * estimates here since this will
+ * considerably increase the
+ * performance of the matrix
+ * setup. However, there should be no
+ * effect in the performance of
+ * matrix-vector products, since
+ * Trilinos wants to reorganize the
+ * matrix memory prior to use.
*/
SparseMatrix (const Epetra_Map &InputMap,
const unsigned int n_max_entries_per_row);
const std::vector<unsigned int> &n_entries_per_row);
/**
- * This constructor is similar to
- * the one above, but it now
- * takes two different Epetra
- * maps for rows and
- * columns. This interface is
- * meant to be used for
- * generating rectangular
- * matrices, where one map takes
- * care of the columns and the
- * other one of the rows. Note
- * that there is no real
- * parallelism along the columns
- * – the processor that
- * owns a certain row always owns
- * all the column elements, no
- * matter how far they might be
- * spread out. The second
- * Epetra_Map is only used to
- * specify the number of columns
- * and for internal arragements
- * when doing matrix-vector
- * products with vectors based on
- * that column map.
+ * This constructor is similar to the
+ * one above, but it now takes two
+ * different Epetra maps for rows and
+ * columns. This interface is meant
+ * to be used for generating
+ * rectangular matrices, where one
+ * map describes the parallel
+ * partitioning of the dofs
+ * associated with the matrix rows
+ * and the other one the partitioning
+ * of dofs in the matrix
+ * columns. Note that there is no
+ * real parallelism along the columns
+ * – the processor that owns a
+ * certain row always owns all the
+ * column elements, no matter how far
+ * they might be spread out. The
+ * second Epetra_Map is only used to
+ * specify the number of columns and
+ * for internal arragements when
+ * doing matrix-vector products with
+ * vectors based on that column map.
*
* The number of columns entries
* per row is specified as the
* to be used for generating
* rectangular matrices, where one
* map specifies the parallel
- * distribution of rows and the
- * second one specifies the number of
- * columns in the total matrix. It
- * also provides information for the
+ * distribution of degrees of freedom
+ * associated with matrix rows and
+ * the second one specifies the
+ * parallel distribution the dofs
+ * associated with columns in the
+ * matrix. The second map also
+ * provides information for the
* internal arrangement in matrix
- * vector products, but is not used
- * for the distribution of the
- * columns – rather, all column
- * elements of a row are stored on
- * the same processor. The vector
+ * vector products (i.e., the
+ * distribution of vector this matrix
+ * is to be multiplied with), but is
+ * not used for the distribution of
+ * the columns – rather, all
+ * column elements of a row are
+ * stored on the same processor in
+ * any case. The vector
* <tt>n_entries_per_row</tt>
* specifies the number of entries in
* each row of the newly generated
const unsigned int n,
const std::vector<unsigned int> &n_entries_per_row);
+ /**
+ * Generate a matrix from a Trilinos
+ * sparsity pattern object.
+ */
+ SparseMatrix (const SparsityPattern &InputSparsityPattern);
+
/**
* Copy constructor. Sets the
* calling matrix to be the same
virtual ~SparseMatrix ();
/**
- * This function initializes
- * the Trilinos matrix with a
- * deal.II sparsity pattern,
- * i.e. it makes the Trilinos
- * Epetra matrix know the
- * position of nonzero entries
+ * This function initializes the
+ * Trilinos matrix with a deal.II
+ * sparsity pattern, i.e. it makes
+ * the Trilinos Epetra matrix know
+ * the position of nonzero entries
* according to the sparsity
- * 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.
+ * 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 different
+ * 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.
*/
template<typename SparsityType>
void reinit (const SparsityType &sparsity_pattern);
/**
- * This function is initializes
- * the Trilinos Epetra matrix
- * according to the specified
- * sparsity_pattern, and also
- * reassigns the matrix rows to
- * different processes
- * according to a user-supplied
- * Epetra map. In programs
- * following the style of the
- * tutorial programs, this
- * function (and the respective
- * call for a rectangular matrix)
- * are the natural way to
- * initialize the matrix size,
- * its distribution among the MPI
- * processes (if run in parallel)
- * as well as the locatoin of
- * non-zero elements. Trilinos
+ * This function is initializes the
+ * Trilinos Epetra matrix according
+ * to the specified sparsity_pattern,
+ * and also reassigns the matrix rows
+ * to different processes according
+ * to a user-supplied Epetra map. In
+ * programs following the style of
+ * the tutorial programs, this
+ * function (and the respective call
+ * for a rectangular matrix) are the
+ * natural way to initialize the
+ * matrix size, its distribution
+ * among the MPI processes (if run in
+ * parallel) as well as the locatoin
+ * of non-zero elements. Trilinos
* stores the sparsity pattern
- * internally, so it won't be
- * needed any more after this
- * call, in contrast to the
- * deal.II own object. In a
- * parallel run, it is currently
- * necessary that each processor
- * holds the sparsity_pattern
- * structure because each
- * processor sets its
- * rows.
+ * internally, so it won't be needed
+ * any more after this call, in
+ * contrast to the deal.II own
+ * object. In a parallel run, 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
+ * This is a collective operation
+ * that needs to be called on all
* processors in order to avoid a
* dead lock.
*/
const SparsityType &sparsity_pattern);
/**
- * This function is similar to
- * the other initialization
- * function above, but now also
- * reassigns the matrix rows
- * and columns according to two
- * user-supplied Epetra maps.
+ * This function is similar to the
+ * other initialization function
+ * above, but now also reassigns the
+ * matrix rows and columns according
+ * to two user-supplied Epetra maps.
* To be used for rectangular
* matrices.
*
- * This is a
- * collective operation that
- * needs to be called on all
+ * This is a collective operation
+ * that needs to be called on all
* processors in order to avoid a
* dead lock.
*/
const SparsityType &sparsity_pattern);
/**
- * This function copies the
- * content in
- * <tt>sparse_matrix</tt> to
- * the calling matrix.
+ * This function reinitializes the
+ * Trilinos sparse matrix from a
+ * (possibly distributed) Trilinos
+ * sparsity pattern.
+ *
+ * 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);
+
+ /**
+ * This function copies the content
+ * in <tt>sparse_matrix</tt> to the
+ * calling matrix.
*
- * This is a
- * collective operation that
- * needs to be called on all
+ * This is a collective operation
+ * that needs to be called on all
* processors in order to avoid a
* dead lock.
*/
void reinit (const SparseMatrix &sparse_matrix);
/**
- * 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).
+ * 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).
*
- * This is a
- * collective operation that
- * needs to be called on all
+ * This is a collective operation
+ * that needs to be called on all
* processors in order to avoid a
* dead lock.
*/
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
+ * 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
+ * This is a collective operation
+ * that needs to be called on all
* processors in order to avoid a
* dead lock.
*/
const double drop_tolerance=1e-13);
/**
- * This function is similar to
- * the other initialization
- * function with deal.II sparse
- * matrix input above, but now
- * takes Epetra maps for both
- * the rows and the columns of
- * the matrix. Chosen for
- * rectangular matrices.
+ * This function is similar to the
+ * other initialization function with
+ * deal.II sparse matrix input above,
+ * but now takes Epetra maps for both
+ * the rows and the columns of the
+ * matrix. Chosen for rectangular
+ * matrices.
*
- * This is a
- * collective operation that
- * needs to be called on all
+ * This is a collective operation
+ * that needs to be called on all
* processors in order to avoid a
* dead lock.
*/
void reinit (const Epetra_CrsMatrix &input_matrix);
/**
- * This operator assigns a scalar
- * to a matrix. Since this does
- * usually not make much sense
- * (should we set all matrix
- * entries to this value? Only
- * the nonzero entries of the
- * sparsity pattern?), this
- * operation is only allowed if
- * the actual value to be
- * assigned is zero. This
- * operator only exists to allow
- * for the obvious notation
- * <tt>matrix=0</tt>, which sets
- * all elements of the matrix to
- * zero, but keeps the sparsity
- * pattern previously used.
+ * This operator assigns a scalar to
+ * a matrix. Since this does usually
+ * not make much sense (should we set
+ * all matrix entries to this value?
+ * Only the nonzero entries of the
+ * sparsity pattern?), this operation
+ * is only allowed if the actual
+ * value to be assigned is zero. This
+ * operator only exists to allow for
+ * the obvious notation
+ * <tt>matrix=0</tt>, which sets all
+ * elements of the matrix to zero,
+ * but keeps the sparsity pattern
+ * previously used.
*/
SparseMatrix &
operator = (const double d);
/**
- * Release all memory and
- * return to a state just like
- * after having called the
- * default constructor.
+ * Release all memory and return to a
+ * state just like after having
+ * called the default constructor.
*
- * This is a
- * collective operation that
- * needs to be called on all
+ * This is a collective operation
+ * that needs to be called on all
* processors in order to avoid a
* dead lock.
*/
void clear ();
/**
- * Trilinos matrices store their
- * own sparsity patterns. So, in
- * analogy to our own
- * SparsityPattern class, this
- * function compresses the
- * sparsity pattern and allows
- * the resulting matrix to be
- * used in all other operations
- * where before only assembly
- * functions were allowed. This
- * function must therefore be
- * called once you have assembled
- * the matrix. This is a
+ * Trilinos matrices store their own
+ * sparsity patterns. So, in analogy
+ * to our own SparsityPattern class,
+ * this function compresses the
+ * sparsity pattern and allows the
+ * resulting matrix to be used in all
+ * other operations where before only
+ * assembly functions were
+ * allowed. This function must
+ * therefore be called once you have
+ * assembled the matrix. This is a
* collective operation, i.e., it
- * needs to be run on all
- * processors when used in
- * parallel.
+ * needs to be run on all processors
+ * when used in parallel.
*/
void compress ();
* i.e., whether compress() needs to
* be called after an operation
* requiring data exchange. A call to
- * compress() is also after the
- * method set() is called (even when
- * working in serial).
+ * compress() is also needed when the
+ * method set() has been called (even
+ * when working in serial).
*/
bool is_compressed () const;
//@}
* to the individual processes.
*/
Epetra_Map col_map;
-
+
/**
* Trilinos doesn't allow to
* mix additions to matrix
#endif // DEAL_II_USE_TRILINOS
-/*---------------------------- trilinos_sparse_matrix.h ---------------------------*/
+/*----------------------- trilinos_sparse_matrix.h --------------------*/
#endif
-/*---------------------------- trilinos_sparse_matrix.h ---------------------------*/
+/*----------------------- trilinos_sparse_matrix.h --------------------*/
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2008 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__trilinos_sparsity_pattern_h
+#define __deal2__trilinos_sparsity_pattern_h
+
+
+#include <base/config.h>
+#include <base/subscriptor.h>
+#include <lac/exceptions.h>
+
+#include <boost/shared_ptr.hpp>
+#include <vector>
+#include <cmath>
+#include <memory>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+# include <Epetra_FECrsGraph.h>
+# include <Epetra_Map.h>
+# ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+# include <Epetra_MpiComm.h>
+# include "mpi.h"
+# else
+# include "Epetra_SerialComm.h"
+# endif
+
+DEAL_II_NAMESPACE_OPEN
+
+ // forward declarations
+class SparsityPattern;
+class CompressedSparsityPattern;
+class CompressedSetSparsityPattern;
+class CompressedSimpleSparsityPattern;
+
+namespace TrilinosWrappers
+{
+ // forward declarations
+ class SparsityPattern;
+
+ namespace SparsityPatternIterators
+ {
+/**
+ * STL conforming iterator. This class acts as an iterator walking
+ * over the elements of Trilinos sparsity pattern.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ */
+ class const_iterator
+ {
+ private:
+ /**
+ * Accessor class for iterators
+ */
+ class Accessor
+ {
+ public:
+ /**
+ * Constructor. Since we use
+ * accessors only for read
+ * access, a const matrix
+ * pointer is sufficient.
+ */
+ Accessor (const SparsityPattern *sparsity_pattern,
+ const unsigned int row,
+ const unsigned int index);
+
+ /**
+ * Row number of the element
+ * represented by this object.
+ */
+ unsigned int row() const;
+
+ /**
+ * Index in row of the element
+ * represented by this object.
+ */
+ unsigned int index() const;
+
+ /**
+ * Column number of the element
+ * represented by this object.
+ */
+ unsigned int column() const;
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcBeyondEndOfSparsityPattern);
+
+ /**
+ * Exception
+ */
+ DeclException3 (ExcAccessToNonlocalRow,
+ int, int, int,
+ << "You tried to access row " << arg1
+ << " of a distributed sparsity pattern, "
+ << " but only rows " << arg2 << " through " << arg3
+ << " are stored locally and can be accessed.");
+
+ private:
+ /**
+ * The matrix accessed.
+ */
+ mutable SparsityPattern *sparsity_pattern;
+
+ /**
+ * Current row number.
+ */
+ unsigned int a_row;
+
+ /**
+ * Current index in row.
+ */
+ unsigned int a_index;
+
+ /**
+ * Cache where we store the
+ * column indices of the
+ * present row. This is
+ * necessary, since Trilinos
+ * makes access to the elements
+ * of its matrices rather hard,
+ * and it is much more
+ * efficient to copy all column
+ * entries of a row once when
+ * we enter it than repeatedly
+ * asking Trilinos for
+ * individual ones. This also
+ * makes some sense since it is
+ * likely that we will access
+ * them sequentially anyway.
+ *
+ * In order to make copying of
+ * iterators/accessor of
+ * acceptable performance, we
+ * keep a shared pointer to
+ * these entries so that more
+ * than one accessor can access
+ * this data if necessary.
+ */
+ boost::shared_ptr<const std::vector<unsigned int> > colnum_cache;
+
+ /**
+ * Discard the old row caches
+ * (they may still be used by
+ * other accessors) and
+ * generate new ones for the
+ * row pointed to presently by
+ * this accessor.
+ */
+ void visit_present_row ();
+
+ /**
+ * Make enclosing class a
+ * friend.
+ */
+ friend class const_iterator;
+ };
+
+ public:
+
+ /**
+ * Constructor. Create an
+ * iterator into the matrix @p
+ * matrix for the given row and
+ * the index within it.
+ */
+ const_iterator (const SparsityPattern *sparsity_pattern,
+ const unsigned int row,
+ const unsigned int index);
+
+ /**
+ * Prefix increment.
+ */
+ const_iterator& operator++ ();
+
+ /**
+ * Postfix increment.
+ */
+ const_iterator operator++ (int);
+
+ /**
+ * Dereferencing operator.
+ */
+ const Accessor& operator* () const;
+
+ /**
+ * Dereferencing operator.
+ */
+ const Accessor* operator-> () const;
+
+ /**
+ * Comparison. True, if both
+ * iterators point to the same
+ * matrix position.
+ */
+ bool operator == (const const_iterator&) const;
+
+ /**
+ * Inverse of <tt>==</tt>.
+ */
+ bool operator != (const const_iterator&) const;
+
+ /**
+ * Comparison operator. Result
+ * is true if either the first
+ * row number is smaller or if
+ * the row numbers are equal
+ * and the first index is
+ * smaller.
+ */
+ bool operator < (const const_iterator&) const;
+
+ /**
+ * Exception
+ */
+ DeclException2 (ExcInvalidIndexWithinRow,
+ int, int,
+ << "Attempt to access element " << arg2
+ << " of row " << arg1
+ << " which doesn't have that many elements.");
+
+ private:
+ /**
+ * Store an object of the
+ * accessor class.
+ */
+ Accessor accessor;
+ };
+
+ }
+
+
+/**
+ * This class implements a wrapper class to use the Trilinos distributed
+ * sparsity pattern class Epetra_FECrsGraph. This class is designed to be
+ * used for construction of parallel Trilinos matrices. The functionality of
+ * this class is modeled after the existing sparsity pattern classes, with
+ * the difference that this class can work fully in parallel according to a
+ * partitioning of the sparsity pattern rows.
+ *
+ * This class has many similarities to the compressed sparsity pattern
+ * classes of deal.II (i.e., the classes CompressedSparsityPattern,
+ * CompressedSetSparsityPattern, and CompressedSimpleSparsityPattern), since
+ * it can dynamically add elements to the pattern without any memory being
+ * previously reserved for it. However, it also has a method
+ * SparsityPattern::compress(), that finalizes the pattern and enables its
+ * use with Trilinos sparse matrices.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Sparsity
+ * @author Martin Kronbichler, 2008
+ */
+ class SparsityPattern : public Subscriptor
+ {
+ public:
+
+ /**
+ * Declare a typedef for the
+ * iterator class.
+ */
+ typedef SparsityPatternIterators::const_iterator const_iterator;
+
+/**
+ * @name Constructors and initalization.
+ */
+//@{
+ /**
+ * Default constructor. Generates an
+ * empty (zero-size) sparsity
+ * pattern.
+ */
+ SparsityPattern ();
+
+ /**
+ * Constructor for a square sparsity
+ * pattern using an Epetra_Map and
+ * the number of nonzero entries in
+ * the rows of the sparsity
+ * pattern. Note that this number
+ * does not need to be exact, and it
+ * is even allowed that the actual
+ * sparsity structure has more
+ * nonzero entries than specified in
+ * the constructor. However it is
+ * still advantageous to provide good
+ * estimates here since this will
+ * considerably increase the
+ * performance when creating the
+ * sparsity pattern.
+ */
+ SparsityPattern (const Epetra_Map &InputMap,
+ const unsigned int n_entries_per_row = 1);
+
+ /**
+ * Same as before, but now use the
+ * exact number of nonzeros in each m
+ * row. Since we know the number of
+ * elements in the sparsity pattern
+ * exactly in this case, we can
+ * already allocate the right amount
+ * of memory, which makes the
+ * creation process by the respective
+ * SparsityPattern::reinit call
+ * considerably faster. However, this
+ * is a rather unusual situation,
+ * since knowing the number of
+ * entries in each row is usually
+ * connected to knowing the indices
+ * of nonzero entries, which the
+ * sparsity pattern is designed to
+ * describe.
+ */
+ SparsityPattern (const Epetra_Map &InputMap,
+ const std::vector<unsigned int> &n_entries_per_row);
+
+ /**
+ * This constructor is similar to the
+ * one above, but it now takes two
+ * different Epetra maps for rows and
+ * columns. This interface is meant
+ * to be used for generating
+ * rectangular sparsity pattern,
+ * where one map describes the
+ * parallel partitioning of the dofs
+ * associated with the sparsity
+ * pattern rows and the other one of
+ * the sparsity pattern columns. Note
+ * that there is no real parallelism
+ * along the columns – the
+ * processor that owns a certain row
+ * always owns all the column
+ * elements, no matter how far they
+ * might be spread out. The second
+ * Epetra_Map is only used to specify
+ * the number of columns and for
+ * internal arragements when doing
+ * matrix-vector products with
+ * vectors based on that column map.
+ *
+ * The number of columns entries
+ * per row is specified as the
+ * maximum number of entries
+ * argument.
+ */
+ SparsityPattern (const Epetra_Map &InputRowMap,
+ const Epetra_Map &InputColMap,
+ const unsigned int n_entries_per_row = 1);
+
+ /**
+ * This constructor is similar to the
+ * one above, but it now takes two
+ * different Epetra maps for rows and
+ * columns. This interface is meant
+ * to be used for generating
+ * rectangular matrices, where one
+ * map specifies the parallel
+ * distribution of rows and the
+ * second one specifies the
+ * distribution of degrees of freedom
+ * associated with matrix
+ * columns. This second map is
+ * however not used for the
+ * distribution of the columns
+ * themselves – rather, all
+ * column elements of a row are
+ * stored on the same processor. The
+ * vector <tt>n_entries_per_row</tt>
+ * specifies the number of entries in
+ * each row of the newly generated
+ * matrix.
+ */
+ SparsityPattern (const Epetra_Map &InputRowMap,
+ const Epetra_Map &InputColMap,
+ const std::vector<unsigned int> &n_entries_per_row);
+
+ /**
+ * Generate a sparsity pattern that
+ * is completely stored locally,
+ * having #m rows and #n columns. The
+ * resulting matrix will be
+ * completely stored locally.
+ *
+ * The number of columns entries per
+ * row is specified as the maximum
+ * number of entries argument. As
+ * above, this does not need to be an
+ * accurate number since the entries
+ * are allocated dynamically in a
+ * similar manner as for the deal.II
+ * CompressedSparsityPattern classes,
+ * but a good estimate will reduce
+ * the setup time of the sparsity
+ * pattern.
+ */
+ SparsityPattern (const unsigned int m,
+ const unsigned int n,
+ const unsigned int n_entries_per_row = 1);
+
+ /**
+ * Generate a sparsity pattern that
+ * is completely stored locally,
+ * having #m rows and #n columns. The
+ * resulting matrix will be
+ * completely stored locally.
+ *
+ * The vector
+ * <tt>n_entries_per_row</tt>
+ * specifies the number of entries in
+ * each row.
+ */
+ SparsityPattern (const unsigned int m,
+ const unsigned int n,
+ const std::vector<unsigned int> &n_entries_per_row);
+
+ /**
+ * Copy constructor. Sets the calling
+ * sparsity pattern to be the same as
+ * the input sparsity pattern.
+ */
+ SparsityPattern (const SparsityPattern &SP);
+
+ /**
+ * Destructor. Made virtual so that
+ * one can use pointers to this
+ * class.
+ */
+ virtual ~SparsityPattern ();
+
+ /**
+ * Reinitialization function for
+ * generating a square sparsity
+ * pattern using an Epetra_Map and
+ * the number of nonzero entries in
+ * the rows of the sparsity
+ * pattern. Note that this number
+ * does not need to be exact, and it
+ * is even allowed that the actual
+ * sparsity structure has more
+ * nonzero entries than specified in
+ * the constructor. However it is
+ * still advantageous to provide good
+ * estimates here since this will
+ * considerably increase the
+ * performance when creating the
+ * sparsity pattern.
+ *
+ * This function does not create any
+ * entries by itself, but provides
+ * the correct data structures that
+ * can be used by the respective
+ * add() function.
+ */
+ void
+ reinit (const Epetra_Map &InputMap,
+ const unsigned int n_entries_per_row = 1);
+
+ /**
+ * Same as before, but now use the
+ * exact number of nonzeros in each m
+ * row. Since we know the number of
+ * elements in the sparsity pattern
+ * exactly in this case, we can
+ * already allocate the right amount
+ * of memory, which makes process of
+ * adding entries to the sparsity
+ * pattern considerably
+ * faster. However, this is a rather
+ * unusual situation, since knowing
+ * the number of entries in each row
+ * is usually connected to knowing
+ * the indices of nonzero entries,
+ * which the sparsity pattern is
+ * designed to describe.
+ */
+ void
+ reinit (const Epetra_Map &InputMap,
+ const std::vector<unsigned int> &n_entries_per_row);
+
+ /**
+ * This reinit function is similar to
+ * the one above, but it now takes
+ * two different Epetra maps for rows
+ * and columns. This interface is
+ * meant to be used for generating
+ * rectangular sparsity pattern,
+ * where one map describes the
+ * parallel partitioning of the dofs
+ * associated with the sparsity
+ * pattern rows and the other one of
+ * the sparsity pattern columns. Note
+ * that there is no real parallelism
+ * along the columns – the
+ * processor that owns a certain row
+ * always owns all the column
+ * elements, no matter how far they
+ * might be spread out. The second
+ * Epetra_Map is only used to specify
+ * the number of columns and for
+ * internal arragements when doing
+ * matrix-vector products with
+ * vectors based on that column map.
+ *
+ * The number of columns entries per
+ * row is specified by the argument
+ * <tt>n_entries_per_row</tt>.
+ */
+ void
+ reinit (const Epetra_Map &InputRowMap,
+ const Epetra_Map &InputColMap,
+ const unsigned int n_entries_per_row = 1);
+
+ /**
+ * This reinit function is similar to
+ * the one above, but it now takes
+ * two different Epetra maps for rows
+ * and columns. This interface is
+ * meant to be used for generating
+ * rectangular matrices, where one
+ * map specifies the parallel
+ * distribution of rows and the
+ * second one specifies the
+ * distribution of degrees of freedom
+ * associated with matrix
+ * columns. This second map is
+ * however not used for the
+ * distribution of the columns
+ * themselves – rather, all
+ * column elements of a row are
+ * stored on the same processor. The
+ * vector <tt>n_entries_per_row</tt>
+ * specifies the number of entries in
+ * each row of the newly generated
+ * matrix.
+ */
+ void
+ reinit (const Epetra_Map &InputRowMap,
+ const Epetra_Map &InputColMap,
+ const std::vector<unsigned int> &n_entries_per_row);
+
+ /**
+ * Initialize a sparsity pattern that
+ * is completely stored locally,
+ * having #m rows and #n columns. The
+ * resulting matrix will be
+ * completely stored locally.
+ *
+ * The number of columns entries per
+ * row is specified as the maximum
+ * number of entries argument. As
+ * above, this does not need to be an
+ * accurate number since the entries
+ * are allocated dynamically in a
+ * similar manner as for the deal.II
+ * CompressedSparsityPattern classes,
+ * but a good estimate will reduce
+ * the setup time of the sparsity
+ * pattern.
+ */
+ void
+ reinit (const unsigned int m,
+ const unsigned int n,
+ const unsigned int n_entries_per_row = 1);
+
+ /**
+ * Initialize a sparsity pattern that
+ * is completely stored locally,
+ * having #m rows and #n columns. The
+ * resulting matrix will be
+ * completely stored locally.
+ *
+ * The vector
+ * <tt>n_entries_per_row</tt>
+ * specifies the number of entries in
+ * each row.
+ */
+ void
+ reinit (const unsigned int m,
+ const unsigned int n,
+ const std::vector<unsigned int> &n_entries_per_row);
+
+ /**
+ * Reinit function. Takes one of the
+ * deal.II sparsity patterns and a
+ * parallel partitioning of the rows
+ * and columns for initializing the
+ * current Trilinos sparsity pattern.
+ */
+ template<typename SparsityType>
+ void
+ reinit (const Epetra_Map &InputRowMap,
+ const Epetra_Map &InputColMap,
+ const SparsityType &SP);
+
+ /**
+ * Reinit function. Takes one of the
+ * deal.II sparsity patterns and a
+ * parallel partitioning of the rows
+ * and columns for initializing the
+ * current Trilinos sparsity pattern.
+ */
+ template<typename SparsityType>
+ void
+ reinit (const Epetra_Map &InputMap,
+ const SparsityType &SP);
+
+ /**
+ * Copy function. Sets the calling
+ * sparsity pattern to be the same as
+ * the input sparsity pattern.
+ */
+ void
+ copy_from (const SparsityPattern &SP);
+
+ /**
+ * Copy function from one of the
+ * deal.II sparsity patterns. If used
+ * in parallel, this function uses an
+ * ad-hoc partitioning of the rows
+ * and columns.
+ */
+ template<typename SparsityType>
+ void
+ copy_from (const SparsityType &SP);
+
+ /**
+ * Release all memory and
+ * return to a state just like
+ * after having called the
+ * default constructor.
+ *
+ * This is a
+ * collective operation that
+ * needs to be called on all
+ * processors in order to avoid a
+ * dead lock.
+ */
+ void clear ();
+
+ /**
+ * In analogy to our own
+ * SparsityPattern class, this
+ * function compresses the sparsity
+ * pattern and allows the resulting
+ * pattern to be used for actually
+ * generating a matrix. This function
+ * also exchanges non-local data that
+ * might have accumulated during the
+ * addition of new elements. This
+ * function must therefore be called
+ * once the structure is fixed. This
+ * is a collective operation, i.e.,
+ * it needs to be run on all
+ * processors when used in parallel.
+ */
+ void compress ();
+//@}
+/**
+ * @name Information on the sparsity pattern
+ */
+//@{
+
+ /**
+ * Returns the state of the sparsity
+ * pattern, i.e., whether compress()
+ * needs to be called after an
+ * operation requiring data
+ * exchange.
+ */
+ bool is_compressed () const;
+
+ /**
+ * Gives the maximum number of
+ * entries per row on the current
+ * processor.
+ */
+ unsigned int max_entries_per_row () const;
+
+ /**
+ * Return the number of rows in this
+ * sparsity pattern.
+ */
+ unsigned int n_rows () const;
+
+ /**
+ * Return the number of columns in
+ * this sparsity pattern.
+ */
+ unsigned int n_cols () const;
+
+ /**
+ * Return the local dimension of the
+ * sparsity pattern, i.e. the number
+ * of rows stored on the present MPI
+ * process. In the sequential case,
+ * this number is the same as
+ * n_rows(), but for parallel
+ * matrices it may be smaller.
+ *
+ * To figure out which elements
+ * exactly are stored locally,
+ * use local_range().
+ */
+ unsigned int local_size () const;
+
+ /**
+ * Return a pair of indices
+ * indicating which rows of this
+ * sparsity pattern are stored
+ * locally. The first number is the
+ * index of the first row stored, the
+ * second the index of the one past
+ * the last one that is stored
+ * locally. If this is a sequential
+ * matrix, then the result will be
+ * the pair (0,n_rows()), otherwise
+ * it will be a pair (i,i+n), where
+ * <tt>n=local_size()</tt>.
+ */
+ std::pair<unsigned int, unsigned int>
+ local_range () const;
+
+ /**
+ * Return whether @p index is
+ * in the local range or not,
+ * see also local_range().
+ */
+ bool in_local_range (const unsigned int index) const;
+
+ /**
+ * Return the number of nonzero
+ * elements of this sparsity pattern.
+ */
+ unsigned int n_nonzero_elements () const;
+
+ /**
+ * Number of entries in a
+ * specific row.
+ */
+ unsigned int row_length (const unsigned int row) const;
+
+ /**
+ * Return whether the object is
+ * empty. It is empty if no memory is
+ * allocated, which is the same as
+ * when both dimensions are zero.
+ */
+ bool empty () const;
+
+ /**
+ * Return whether the index
+ * (<i>i,j</i>) exists in the
+ * sparsity pattern (i.e., it may be
+ * non-zero) or not.
+ */
+ bool exists (const unsigned int i,
+ const unsigned int j) const;
+
+ /**
+ * Determine an estimate for the
+ * memory consumption (in bytes)
+ * of this object. Currently not
+ * implemented for this class.
+ */
+ unsigned int memory_consumption () const;
+
+//@}
+/**
+ * @name Adding entries
+ */
+//@{
+ /**
+ * Add the element (<i>i,j</i>) to
+ * the sparsity pattern.
+ */
+ void add (const unsigned int i,
+ const unsigned int j);
+
+
+ /**
+ * Add several elements in one row to
+ * the sparsity pattern.
+ */
+ void add (const unsigned int row,
+ const unsigned int n_cols,
+ const unsigned int *col_indices);
+//@}
+/**
+ * @name Iterators
+ */
+//@{
+
+ /**
+ * STL-like iterator with the
+ * first entry.
+ */
+ const_iterator begin () const;
+
+ /**
+ * Final iterator.
+ */
+ const_iterator end () const;
+
+ /**
+ * STL-like iterator with the
+ * first entry of row @p r.
+ *
+ * Note that if the given row
+ * is empty, i.e. does not
+ * contain any nonzero entries,
+ * then the iterator returned
+ * by this function equals
+ * <tt>end(r)</tt>. Note also
+ * that the iterator may not be
+ * dereferencable in that case.
+ */
+ const_iterator begin (const unsigned int r) const;
+
+ /**
+ * Final iterator of row
+ * <tt>r</tt>. It points to the
+ * first element past the end
+ * of line @p r, or past the
+ * end of the entire sparsity
+ * pattern.
+ *
+ * Note that the end iterator
+ * is not necessarily
+ * dereferencable. This is in
+ * particular the case if it is
+ * the end iterator for the
+ * last row of a matrix.
+ */
+ const_iterator end (const unsigned int r) const;
+
+//@}
+/**
+ * @name Input/Output
+ */
+//@{
+
+ /**
+ * Abstract Trilinos object
+ * that helps view in ASCII
+ * other Trilinos
+ * objects. Currently this
+ * function is not
+ * implemented. TODO: Not
+ * implemented.
+ */
+ void write_ascii ();
+
+ /**
+ * Print the sparsity pattern to the
+ * given stream, using the format
+ * <tt>(line,col)</tt>.
+ */
+ void print (std::ostream &out) const;
+
+ // TODO: Write an overloading
+ // of the operator << for output.
+ // Since the underlying Trilinos
+ // object supports it, this should
+ // be very easy.
+
+//@}
+ /** @addtogroup Exceptions
+ * @{ */
+ /**
+ * Exception
+ */
+ DeclException1 (ExcTrilinosError,
+ int,
+ << "An error with error number " << arg1
+ << " occured while calling a Trilinos function");
+
+ /**
+ * Exception
+ */
+ DeclException2 (ExcInvalidIndex,
+ int, int,
+ << "The entry with index <" << arg1 << ',' << arg2
+ << "> does not exist.");
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcSourceEqualsDestination);
+
+ /**
+ * Exception
+ */
+ DeclException4 (ExcAccessToNonLocalElement,
+ int, int, int, int,
+ << "You tried to access element (" << arg1
+ << "/" << arg2 << ")"
+ << " of a distributed matrix, but only rows "
+ << arg3 << " through " << arg4
+ << " are stored locally and can be accessed.");
+
+ /**
+ * Exception
+ */
+ DeclException2 (ExcAccessToNonPresentElement,
+ int, int,
+ << "You tried to access element (" << arg1
+ << "/" << arg2 << ")"
+ << " of a sparse matrix, but it appears to not"
+ << " exist in the Trilinos sparsity pattern.");
+ //@}
+ private:
+ /**
+ * Epetra Trilinos
+ * mapping of the matrix rows
+ * that assigns parts of the
+ * matrix to the individual
+ * processes. This map is
+ * provided either via the
+ * constructor or in a reinit
+ * function.
+ */
+ Epetra_Map row_map;
+
+ /**
+ * Pointer to the user-supplied
+ * Epetra Trilinos mapping of
+ * the matrix columns that
+ * assigns parts of the matrix
+ * to the individual processes.
+ */
+ Epetra_Map col_map;
+
+ /**
+ * A boolean variable to hold
+ * information on whether the
+ * vector is compressed or not.
+ */
+ bool compressed;
+
+ /**
+ * A sparsity pattern object in
+ * Trilinos to be used for finite
+ * element based problems which
+ * allows for adding non-local
+ * elements to the pattern.
+ */
+ std::auto_ptr<Epetra_FECrsGraph> graph;
+
+ friend class SparseMatrix;
+ friend class SparsityPatternIterators::const_iterator;
+ };
+
+
+
+// -------------------------- inline and template functions ----------------------
+
+
+#ifndef DOXYGEN
+
+ namespace SparsityPatternIterators
+ {
+
+ inline
+ const_iterator::Accessor::
+ Accessor (const SparsityPattern *sp,
+ const unsigned int row,
+ const unsigned int index)
+ :
+ sparsity_pattern(const_cast<SparsityPattern*>(sp)),
+ a_row(row),
+ a_index(index)
+ {
+ visit_present_row ();
+ }
+
+
+ inline
+ unsigned int
+ const_iterator::Accessor::row() const
+ {
+ Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern());
+ return a_row;
+ }
+
+
+
+ inline
+ unsigned int
+ const_iterator::Accessor::column() const
+ {
+ Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern());
+ return (*colnum_cache)[a_index];
+ }
+
+
+
+ inline
+ unsigned int
+ const_iterator::Accessor::index() const
+ {
+ Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern());
+ return a_index;
+ }
+
+
+
+ inline
+ const_iterator::
+ const_iterator(const SparsityPattern *sp,
+ const unsigned int row,
+ const unsigned int index)
+ :
+ accessor(sp, row, index)
+ {}
+
+
+
+ inline
+ const_iterator &
+ const_iterator::operator++ ()
+ {
+ Assert (accessor.a_row < accessor.sparsity_pattern->n_rows(),
+ ExcIteratorPastEnd());
+
+ ++accessor.a_index;
+
+ // If at end of line: do one
+ // step, then cycle until we
+ // find a row with a nonzero
+ // number of entries.
+ if (accessor.a_index >= accessor.colnum_cache->size())
+ {
+ accessor.a_index = 0;
+ ++accessor.a_row;
+
+ while ((accessor.a_row < accessor.sparsity_pattern->n_rows())
+ &&
+ (accessor.sparsity_pattern->row_length(accessor.a_row) == 0))
+ ++accessor.a_row;
+
+ accessor.visit_present_row();
+ }
+ return *this;
+ }
+
+
+
+ inline
+ const_iterator
+ const_iterator::operator++ (int)
+ {
+ const const_iterator old_state = *this;
+ ++(*this);
+ return old_state;
+ }
+
+
+
+ inline
+ const const_iterator::Accessor &
+ const_iterator::operator* () const
+ {
+ return accessor;
+ }
+
+
+
+ inline
+ const const_iterator::Accessor *
+ const_iterator::operator-> () const
+ {
+ return &accessor;
+ }
+
+
+
+ inline
+ bool
+ const_iterator::
+ operator == (const const_iterator& other) const
+ {
+ return (accessor.a_row == other.accessor.a_row &&
+ accessor.a_index == other.accessor.a_index);
+ }
+
+
+
+ inline
+ bool
+ const_iterator::
+ operator != (const const_iterator& other) const
+ {
+ return ! (*this == other);
+ }
+
+
+
+ inline
+ bool
+ const_iterator::
+ operator < (const const_iterator& other) const
+ {
+ return (accessor.row() < other.accessor.row() ||
+ (accessor.row() == other.accessor.row() &&
+ accessor.index() < other.accessor.index()));
+ }
+
+ }
+
+
+
+ inline
+ SparsityPattern::const_iterator
+ SparsityPattern::begin() const
+ {
+ return const_iterator(this, 0, 0);
+ }
+
+
+
+ inline
+ SparsityPattern::const_iterator
+ SparsityPattern::end() const
+ {
+ return const_iterator(this, n_rows(), 0);
+ }
+
+
+
+ inline
+ SparsityPattern::const_iterator
+ SparsityPattern::begin(const unsigned int r) const
+ {
+ Assert (r < n_rows(), ExcIndexRange(r, 0, n_rows()));
+ if (row_length(r) > 0)
+ return const_iterator(this, r, 0);
+ else
+ return end (r);
+ }
+
+
+
+ inline
+ SparsityPattern::const_iterator
+ SparsityPattern::end(const unsigned int r) const
+ {
+ Assert (r < n_rows(), ExcIndexRange(r, 0, n_rows()));
+
+ // place the iterator on the first entry
+ // past this line, or at the end of the
+ // matrix
+ for (unsigned int i=r+1; i<n_rows(); ++i)
+ if (row_length(i) > 0)
+ return const_iterator(this, i, 0);
+
+ // if there is no such line, then take the
+ // end iterator of the matrix
+ return end();
+ }
+
+
+
+ inline
+ bool
+ SparsityPattern::in_local_range (const unsigned int index) const
+ {
+ int begin, end;
+ begin = graph->RowMap().MinMyGID();
+ end = graph->RowMap().MaxMyGID()+1;
+
+ return ((index >= static_cast<unsigned int>(begin)) &&
+ (index < static_cast<unsigned int>(end)));
+ }
+
+
+
+ inline
+ bool
+ SparsityPattern::is_compressed () const
+ {
+ return compressed;
+ }
+
+
+
+ inline
+ bool
+ SparsityPattern::empty () const
+ {
+ return ((n_rows() == 0) && (n_cols() == 0));
+ }
+
+
+
+ inline
+ void
+ SparsityPattern::add (const unsigned int i,
+ const unsigned int j)
+ {
+ add (i, 1, &j);
+ }
+
+
+
+ inline
+ void
+ SparsityPattern::add (const unsigned int row,
+ const unsigned int n_cols,
+ const unsigned int *col_indices)
+ {
+ int * col_index_ptr = (int*)col_indices;
+ compressed = false;
+
+ int ierr;
+
+ // If the calling sparsity pattern owns
+ // the row to which we want to add
+ // values, we can directly call the
+ // Epetra_CrsGraph input function, which
+ // is much faster than the
+ // Epetra_FECrsGraph function.
+ if (row_map.MyGID(row) == true)
+ ierr = graph->Epetra_CrsGraph::InsertGlobalIndices(row,
+ n_cols,
+ col_index_ptr);
+ else
+ {
+ // When we're at off-processor data, we
+ // have to stick with the standard
+ // SumIntoGlobalValues
+ // function. Nevertheless, the way we
+ // call it is the fastest one (any other
+ // will lead to repeated allocation and
+ // deallocation of memory in order to
+ // call the function we already use,
+ // which is very unefficient if writing
+ // one element at a time).
+
+ ierr = graph->InsertGlobalIndices (1, (int*)&row, n_cols,
+ col_index_ptr);
+ }
+
+ //Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
+ AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
+ }
+
+
+#endif // DOXYGEN
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif // DEAL_II_USE_TRILINOS
+
+
+/*-------------------- trilinos_sparsity_pattern.h --------------------*/
+
+#endif
+/*-------------------- trilinos_sparsity_pattern.h --------------------*/
template class BlockSparsityPatternBase<CompressedSparsityPattern>;
template class BlockSparsityPatternBase<CompressedSimpleSparsityPattern>;
template class BlockSparsityPatternBase<CompressedSetSparsityPattern>;
-
+#ifdef DEAL_II_USE_TRILINOS
+template class BlockSparsityPatternBase<TrilinosWrappers::SparsityPattern>;
+#endif
this->collect_sizes();
}
-BlockCompressedSimpleSparsityPattern::BlockCompressedSimpleSparsityPattern ()
+
+
+BlockCompressedSetSparsityPattern::BlockCompressedSetSparsityPattern ()
{}
-BlockCompressedSimpleSparsityPattern::
-BlockCompressedSimpleSparsityPattern (const unsigned int n_rows,
+BlockCompressedSetSparsityPattern::
+BlockCompressedSetSparsityPattern (const unsigned int n_rows,
const unsigned int n_columns)
:
- BlockSparsityPatternBase<CompressedSimpleSparsityPattern>(n_rows,
+ BlockSparsityPatternBase<CompressedSetSparsityPattern>(n_rows,
n_columns)
{}
-BlockCompressedSimpleSparsityPattern::
-BlockCompressedSimpleSparsityPattern (const std::vector<unsigned int>& row_indices,
+BlockCompressedSetSparsityPattern::
+BlockCompressedSetSparsityPattern (const std::vector<unsigned int>& row_indices,
const std::vector<unsigned int>& col_indices)
:
- BlockSparsityPatternBase<CompressedSimpleSparsityPattern>(row_indices.size(),
+ BlockSparsityPatternBase<CompressedSetSparsityPattern>(row_indices.size(),
col_indices.size())
{
for (unsigned int i=0;i<row_indices.size();++i)
void
-BlockCompressedSimpleSparsityPattern::reinit (
+BlockCompressedSetSparsityPattern::reinit (
const std::vector< unsigned int > &row_block_sizes,
const std::vector< unsigned int > &col_block_sizes)
{
- BlockSparsityPatternBase<CompressedSimpleSparsityPattern>::reinit(row_block_sizes.size(), col_block_sizes.size());
+ BlockSparsityPatternBase<CompressedSetSparsityPattern>::reinit(row_block_sizes.size(), col_block_sizes.size());
for (unsigned int i=0;i<row_block_sizes.size();++i)
for (unsigned int j=0;j<col_block_sizes.size();++j)
this->block(i,j).reinit(row_block_sizes[i],col_block_sizes[j]);
-BlockCompressedSetSparsityPattern::BlockCompressedSetSparsityPattern ()
+BlockCompressedSimpleSparsityPattern::BlockCompressedSimpleSparsityPattern ()
{}
-BlockCompressedSetSparsityPattern::
-BlockCompressedSetSparsityPattern (const unsigned int n_rows,
+BlockCompressedSimpleSparsityPattern::
+BlockCompressedSimpleSparsityPattern (const unsigned int n_rows,
const unsigned int n_columns)
:
- BlockSparsityPatternBase<CompressedSetSparsityPattern>(n_rows,
+ BlockSparsityPatternBase<CompressedSimpleSparsityPattern>(n_rows,
n_columns)
{}
-BlockCompressedSetSparsityPattern::
-BlockCompressedSetSparsityPattern (const std::vector<unsigned int>& row_indices,
+
+BlockCompressedSimpleSparsityPattern::
+BlockCompressedSimpleSparsityPattern (const std::vector<unsigned int>& row_indices,
const std::vector<unsigned int>& col_indices)
:
- BlockSparsityPatternBase<CompressedSetSparsityPattern>(row_indices.size(),
+ BlockSparsityPatternBase<CompressedSimpleSparsityPattern>(row_indices.size(),
col_indices.size())
{
for (unsigned int i=0;i<row_indices.size();++i)
}
+
void
-BlockCompressedSetSparsityPattern::reinit (
+BlockCompressedSimpleSparsityPattern::reinit (
const std::vector< unsigned int > &row_block_sizes,
const std::vector< unsigned int > &col_block_sizes)
{
- BlockSparsityPatternBase<CompressedSetSparsityPattern>::reinit(row_block_sizes.size(), col_block_sizes.size());
+ BlockSparsityPatternBase<CompressedSimpleSparsityPattern>::reinit(row_block_sizes.size(), col_block_sizes.size());
for (unsigned int i=0;i<row_block_sizes.size();++i)
for (unsigned int j=0;j<col_block_sizes.size();++j)
this->block(i,j).reinit(row_block_sizes[i],col_block_sizes[j]);
+#ifdef DEAL_II_USE_TRILINOS
+namespace TrilinosWrappers
+{
+
+ BlockSparsityPattern::BlockSparsityPattern ()
+ {}
+
+
+
+ BlockSparsityPattern::
+ BlockSparsityPattern (const unsigned int n_rows,
+ const unsigned int n_columns)
+ :
+ dealii::BlockSparsityPatternBase<SparsityPattern>(n_rows,
+ n_columns)
+ {}
+
+
+
+ BlockSparsityPattern::
+ BlockSparsityPattern (const std::vector<unsigned int>& row_indices,
+ const std::vector<unsigned int>& col_indices)
+ :
+ BlockSparsityPatternBase<SparsityPattern>(row_indices.size(),
+ col_indices.size())
+ {
+ for (unsigned int i=0;i<row_indices.size();++i)
+ for (unsigned int j=0;j<col_indices.size();++j)
+ this->block(i,j).reinit(row_indices[i],col_indices[j]);
+ this->collect_sizes();
+ }
+
+
+
+ BlockSparsityPattern::
+ BlockSparsityPattern (const std::vector<Epetra_Map>& input_maps)
+ :
+ BlockSparsityPatternBase<SparsityPattern>(input_maps.size(),
+ input_maps.size())
+ {
+ for (unsigned int i=0;i<input_maps.size();++i)
+ for (unsigned int j=0;j<input_maps.size();++j)
+ this->block(i,j).reinit(input_maps[i],input_maps[j]);
+ this->collect_sizes();
+ }
+
+
+
+ void
+ BlockSparsityPattern::reinit (const std::vector<unsigned int> &row_block_sizes,
+ const std::vector<unsigned int> &col_block_sizes)
+ {
+ dealii::BlockSparsityPatternBase<SparsityPattern>::reinit(row_block_sizes.size(), col_block_sizes.size());
+ for (unsigned int i=0;i<row_block_sizes.size();++i)
+ for (unsigned int j=0;j<col_block_sizes.size();++j)
+ this->block(i,j).reinit(row_block_sizes[i],col_block_sizes[j]);
+ this->collect_sizes();
+ }
+
+
+
+ void
+ BlockSparsityPattern::reinit (const std::vector<Epetra_Map> &input_maps)
+ {
+ dealii::BlockSparsityPatternBase<SparsityPattern>::reinit(input_maps.size(),
+ input_maps.size());
+ for (unsigned int i=0;i<input_maps.size();++i)
+ for (unsigned int j=0;j<input_maps.size();++j)
+ this->block(i,j).reinit(input_maps[i],input_maps[j]);
+ this->collect_sizes();
+ }
+
+}
+
+#endif
// Remark: The explicit instantiations for "BlockSparsityPatternBase" were moved
// to the top of this source file. The reason is a slightly buggy version
+ template <>
+ void
+ BlockSparseMatrix::
+ reinit (const BlockSparsityPattern &block_sparsity_pattern)
+ {
+
+ // Call the other basic reinit function, ...
+ reinit (block_sparsity_pattern.n_block_rows(),
+ block_sparsity_pattern.n_block_cols());
+
+ // ... set the correct sizes, ...
+ this->row_block_indices = block_sparsity_pattern.get_row_indices();
+ this->column_block_indices = block_sparsity_pattern.get_column_indices();
+
+ // ... and then assign the correct
+ // data to the blocks.
+ for (unsigned int r=0; r<this->n_block_rows(); ++r)
+ for (unsigned int c=0; c<this->n_block_cols(); ++c)
+ {
+ this->sub_objects[r][c]->reinit (block_sparsity_pattern.block(r,c));
+ }
+ }
+
+
+
void
BlockSparseMatrix::
reinit (const std::vector<Epetra_Map> &input_maps,
// -------------------- explicit instantiations -----------------------
//
template void
- BlockSparseMatrix::reinit (const BlockSparsityPattern &);
+ BlockSparseMatrix::reinit (const dealii::BlockSparsityPattern &);
template void
BlockSparseMatrix::reinit (const BlockCompressedSparsityPattern &);
template void
template void
BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
- const BlockSparsityPattern &);
+ const dealii::BlockSparsityPattern &);
template void
BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
const BlockCompressedSparsityPattern &);
break;
case gmres:
aztecBlockParams->sublist("Forward Solve")
- .sublist("AztecOO Settings").set("Aztec Solver", "GMRES");
+ .sublist("AztecOO Settings").set("Aztec Solver", "MinRes");
aztecBlockParams->sublist("Forward Solve")
.sublist("AztecOO Settings").set("Size of Krylov Subspace",
(int)additional_data.gmres_restart_parameter);
#include <lac/trilinos_sparse_matrix.h>
+#include <lac/trilinos_sparsity_pattern.h>
#include <lac/sparsity_pattern.h>
#include <lac/compressed_sparsity_pattern.h>
#include <lac/compressed_set_sparsity_pattern.h>
matrix->FillComplete();
}
+
+
SparseMatrix::SparseMatrix (const Epetra_Map &InputMap,
const unsigned int n_max_entries_per_row)
:
int(n_max_entries_per_row), false)))
{}
+
+
SparseMatrix::SparseMatrix (const Epetra_Map &InputMap,
const std::vector<unsigned int> &n_entries_per_row)
:
false)))
{}
+
+
SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap,
const Epetra_Map &InputColMap,
const unsigned int n_max_entries_per_row)
int(n_max_entries_per_row), false)))
{}
+
+
SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap,
const Epetra_Map &InputColMap,
const std::vector<unsigned int> &n_entries_per_row)
false)))
{}
+
+
SparseMatrix::SparseMatrix (const unsigned int m,
const unsigned int n,
const unsigned int n_max_entries_per_row)
int(n_max_entries_per_row), false)))
{}
+
+
SparseMatrix::SparseMatrix (const unsigned int m,
const unsigned int n,
const std::vector<unsigned int> &n_entries_per_row)
false)))
{}
+
+
+ SparseMatrix::SparseMatrix (const SparsityPattern &InputSP)
+ :
+ Subscriptor(),
+ row_map (InputSP.row_map),
+ col_map (InputSP.col_map),
+ last_action (Zero),
+ compressed (true),
+ matrix (std::auto_ptr<Epetra_FECrsMatrix>
+ (new Epetra_FECrsMatrix(Copy, *InputSP.graph, false)))
+ {
+ Assert(InputSP.graph->Filled() == true,
+ ExcMessage("The Trilinos sparsity pattern has not been compressed."));
+ compress();
+ }
+
+
+
SparseMatrix::SparseMatrix (const SparseMatrix &InputMatrix)
:
Subscriptor(),
+ void
+ SparseMatrix::reinit (const SparsityPattern &sparsity_pattern)
+ {
+ matrix.reset();
+
+ row_map = sparsity_pattern.row_map;
+ col_map = sparsity_pattern.col_map;
+
+ Assert (sparsity_pattern.graph->Filled() == true,
+ ExcMessage("The Trilinos sparsity pattern has not been compressed"));
+
+ matrix = std::auto_ptr<Epetra_FECrsMatrix>
+ (new Epetra_FECrsMatrix(Copy, *sparsity_pattern.graph, false));
+
+ compress();
+ }
+
+
+
void
SparseMatrix::reinit (const SparseMatrix &sparse_matrix)
{
// If the data is not on the
// present processor, we throw
- // an exception. This is on of
+ // an exception. This is one of
// the two tiny differences to
// the el(i,j) call, which does
// not throw any assertions.
- if ((trilinos_i == -1 ) || (trilinos_j == -1))
+ if (trilinos_i == -1)
{
Assert (false, ExcAccessToNonLocalElement(i, j, local_range().first,
local_range().second));
}
else
{
- // Check whether the matrix
- // already is transformed to
- // local indices.
+ // Check whether the matrix has
+ // already been transformed to local
+ // indices.
if (matrix->Filled() == false)
matrix->GlobalAssemble(col_map, row_map, true);
unsigned int
SparseMatrix::n () const
{
- int n_cols = matrix -> NumGlobalCols();
+ unsigned int n_cols = matrix -> NumGlobalCols();
return n_cols;
}
unsigned int
SparseMatrix::local_size () const
{
- int n_rows = matrix -> NumMyRows();
+ unsigned int n_rows = matrix -> NumMyRows();
return n_rows;
}
std::pair<unsigned int, unsigned int>
SparseMatrix::local_range () const
{
- int begin, end;
+ unsigned int begin, end;
begin = matrix -> RowMap().MinMyGID();
end = matrix -> RowMap().MaxMyGID()+1;
unsigned int
SparseMatrix::n_nonzero_elements () const
{
- int nnz = matrix->NumGlobalNonzeros();
+ unsigned int nnz = matrix->NumGlobalNonzeros();
- return static_cast<unsigned int>(nnz);
+ return nnz;
}
// get a representation of the
// present row
int ncols = -1;
- int local_row = matrix->RowMap().LID(row);
+ int local_row = matrix->LRID(row);
// on the processor who owns this
// row, we'll have a non-negative
}
+
TrilinosScalar
SparseMatrix::l1_norm () const
{
// explicit instantiations
//
template void
- SparseMatrix::reinit (const SparsityPattern &);
+ SparseMatrix::reinit (const dealii::SparsityPattern &);
template void
SparseMatrix::reinit (const CompressedSparsityPattern &);
template void
template void
SparseMatrix::reinit (const Epetra_Map &,
- const SparsityPattern &);
+ const dealii::SparsityPattern &);
template void
SparseMatrix::reinit (const Epetra_Map &,
const CompressedSparsityPattern &);
template void
SparseMatrix::reinit (const Epetra_Map &,
const Epetra_Map &,
- const SparsityPattern &);
+ const dealii::SparsityPattern &);
template void
SparseMatrix::reinit (const Epetra_Map &,
const Epetra_Map &,
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2008 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <lac/trilinos_sparsity_pattern.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
+
+namespace TrilinosWrappers
+{
+ namespace SparsityPatternIterators
+ {
+ void
+ SparsityPattern::const_iterator::Accessor::
+ visit_present_row ()
+ {
+ // if we are asked to visit the
+ // past-the-end line, then simply
+ // release all our caches and go on
+ // with life
+ if (this->a_row == sparsity_pattern->n_rows())
+ {
+ colnum_cache.reset ();
+
+ return;
+ }
+
+ // otherwise first flush Trilinos caches
+ sparsity_pattern->compress ();
+
+ // get a representation of the present
+ // row
+ int ncols;
+ int colnums = sparsity_pattern->n_cols();
+
+ int ierr;
+ ierr = sparsity_pattern->graph->ExtractGlobalRowCopy((int)this->a_row,
+ colnums,
+ ncols,
+ (int*)&(*colnum_cache)[0]);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ // copy it into our caches if the
+ // line isn't empty. if it is, then
+ // we've done something wrong, since
+ // we shouldn't have initialized an
+ // iterator for an empty line (what
+ // would it point to?)
+ Assert (ncols != 0, ExcInternalError());
+ colnum_cache.reset (new std::vector<unsigned int> (colnums,
+ colnums+ncols));
+ }
+ }
+
+
+ // The constructor is actually the
+ // only point where we have to check
+ // whether we build a serial or a
+ // parallel Trilinos matrix.
+ // Actually, it does not even matter
+ // how many threads there are, but
+ // only if we use an MPI compiler or
+ // a standard compiler. So, even one
+ // thread on a configuration with
+ // MPI will still get a parallel
+ // interface.
+ SparsityPattern::SparsityPattern ()
+ :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ row_map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+#else
+ row_map (0, 0, Epetra_SerialComm()),
+#endif
+ col_map (row_map),
+ compressed (true),
+ graph (std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(View, row_map, 0)))
+ {
+ graph->FillComplete();
+ }
+
+ SparsityPattern::SparsityPattern (const Epetra_Map &InputMap,
+ const unsigned int n_entries_per_row)
+ :
+ row_map (InputMap),
+ col_map (row_map),
+ compressed (false),
+ graph (std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(Copy, row_map,
+ int(n_entries_per_row), false)))
+ {}
+
+ SparsityPattern::SparsityPattern (const Epetra_Map &InputMap,
+ const std::vector<unsigned int> &n_entries_per_row)
+ :
+ row_map (InputMap),
+ col_map (row_map),
+ compressed (false),
+ graph (std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(Copy, row_map,
+ (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
+ false)))
+ {}
+
+ SparsityPattern::SparsityPattern (const Epetra_Map &InputRowMap,
+ const Epetra_Map &InputColMap,
+ const unsigned int n_entries_per_row)
+ :
+ row_map (InputRowMap),
+ col_map (InputColMap),
+ compressed (false),
+ graph (std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(Copy, row_map,
+ int(n_entries_per_row), false)))
+ {}
+
+ SparsityPattern::SparsityPattern (const Epetra_Map &InputRowMap,
+ const Epetra_Map &InputColMap,
+ const std::vector<unsigned int> &n_entries_per_row)
+ :
+ row_map (InputRowMap),
+ col_map (InputColMap),
+ compressed (false),
+ graph (std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(Copy, row_map,
+ (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
+ false)))
+ {}
+
+ SparsityPattern::SparsityPattern (const unsigned int m,
+ const unsigned int n,
+ const unsigned int n_entries_per_row)
+ :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ row_map (m, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+ col_map (n, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+#else
+ row_map (m, 0, Epetra_SerialComm()),
+ col_map (n, 0, Epetra_SerialComm()),
+#endif
+ compressed (false),
+ graph (std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(Copy, row_map,
+ int(n_entries_per_row), false)))
+ {}
+
+ SparsityPattern::SparsityPattern (const unsigned int m,
+ const unsigned int n,
+ const std::vector<unsigned int> &n_entries_per_row)
+ :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ row_map (m, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+ col_map (n, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+#else
+ row_map (m, 0, Epetra_SerialComm()),
+ col_map (n, 0, Epetra_SerialComm()),
+#endif
+ compressed (false),
+ graph (std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(Copy, row_map,
+ (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
+ false)))
+ {}
+
+ // Copy function is currently not working
+ // because the Trilinos Epetra_FECrsGraph
+ // does not implement a reinit function
+ // from another graph.
+ /*
+ SparsityPattern::SparsityPattern (const SparsityPattern &InputSP)
+ :
+ Subscriptor(),
+ row_map (InputSP.row_map),
+ col_map (InputSP.col_map),
+ compressed (false),
+ graph (std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(*InputSP.graph)))
+ {}
+ */
+
+
+
+ SparsityPattern::~SparsityPattern ()
+ {}
+
+
+
+ void
+ SparsityPattern::reinit (const Epetra_Map &input_map,
+ const unsigned int n_entries_per_row)
+ {
+ reinit (input_map, input_map, n_entries_per_row);
+ }
+
+
+ void
+ SparsityPattern::reinit (const unsigned int m,
+ const unsigned int n,
+ const unsigned int n_entries_per_row)
+ {
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD);
+#else
+ Epetra_SerialComm trilinos_communicator;
+#endif
+
+ const Epetra_Map rows (m, 0, trilinos_communicator);
+ const Epetra_Map columns (n, 0, trilinos_communicator);
+
+ reinit (rows, columns, n_entries_per_row);
+ }
+
+
+ void
+ SparsityPattern::reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
+ const unsigned int n_entries_per_row)
+ {
+ graph.reset();
+
+ row_map = input_row_map;
+ col_map = input_col_map;
+
+ graph = std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(Copy, row_map, n_entries_per_row, false));
+ }
+
+
+
+ void
+ SparsityPattern::reinit (const Epetra_Map &input_map,
+ const std::vector<unsigned int> &n_entries_per_row)
+ {
+ reinit (input_map, input_map, n_entries_per_row);
+ }
+
+
+
+ void
+ SparsityPattern::reinit (const unsigned int m,
+ const unsigned int n,
+ const std::vector<unsigned int> &n_entries_per_row)
+ {
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD);
+#else
+ Epetra_SerialComm trilinos_communicator;
+#endif
+
+ const Epetra_Map rows (m, 0, trilinos_communicator);
+ const Epetra_Map columns (n, 0, trilinos_communicator);
+
+ reinit (rows, columns, n_entries_per_row);
+ }
+
+
+
+ void
+ SparsityPattern::reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
+ const std::vector<unsigned int> &n_entries_per_row)
+ {
+ graph.reset();
+
+ Assert (n_entries_per_row.size() ==
+ static_cast<unsigned int>(input_row_map.NumGlobalElements()),
+ ExcDimensionMismatch (n_entries_per_row.size(),
+ input_row_map.NumGlobalElements()));
+ row_map = input_row_map;
+ col_map = input_col_map;
+
+ graph = std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(Copy, row_map,
+ n_entries_per_row[input_row_map.MinMyGID()],
+ false));
+ }
+
+
+
+ template <typename SparsityType>
+ void
+ SparsityPattern::reinit (const Epetra_Map &input_map,
+ const SparsityType &sp)
+ {
+ reinit (input_map, input_map, sp);
+ }
+
+
+
+ template <typename SparsityType>
+ void
+ SparsityPattern::reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
+ const SparsityType &sp)
+ {
+ graph.reset();
+
+ Assert (sp.n_rows() ==
+ static_cast<unsigned int>(input_row_map.NumGlobalElements()),
+ ExcDimensionMismatch (sp.n_rows(),
+ input_row_map.NumGlobalElements()));
+ Assert (sp.n_cols() ==
+ static_cast<unsigned int>(input_col_map.NumGlobalElements()),
+ ExcDimensionMismatch (sp.n_cols(),
+ input_col_map.NumGlobalElements()));
+
+ row_map = input_row_map;
+ col_map = input_col_map;
+
+ const unsigned int n_rows = sp.n_rows();
+
+ std::vector<int> n_entries_per_row(n_rows);
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ n_entries_per_row[row] = sp.row_length(row);
+
+ graph = std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(Copy, row_map,
+ n_entries_per_row[input_row_map.MinMyGID()],
+ false));
+
+ Assert (graph->NumGlobalRows() == (int)sp.n_rows(),
+ ExcDimensionMismatch (graph->NumGlobalRows(),
+ sp.n_rows()));
+
+
+ std::vector<int> row_indices;
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ if (row_map.MyGID(row))
+ {
+ const int row_length = sp.row_length(row);
+ row_indices.resize (row_length, -1);
+
+ for (int col=0; col < row_length; ++col)
+ row_indices[col] = sp.column_number (row, col);
+
+ graph->Epetra_CrsGraph::InsertGlobalIndices (row, row_length,
+ &row_indices[0]);
+ }
+
+ compress();
+ }
+
+
+
+ template<>
+ void
+ SparsityPattern::reinit (const Epetra_Map &input_row_map,
+ const Epetra_Map &input_col_map,
+ const CompressedSetSparsityPattern &sp)
+ {
+ graph.reset();
+
+ Assert (sp.n_rows() ==
+ static_cast<unsigned int>(input_row_map.NumGlobalElements()),
+ ExcDimensionMismatch (sp.n_rows(),
+ input_row_map.NumGlobalElements()));
+ Assert (sp.n_cols() ==
+ static_cast<unsigned int>(input_col_map.NumGlobalElements()),
+ ExcDimensionMismatch (sp.n_cols(),
+ input_col_map.NumGlobalElements()));
+
+ row_map = input_row_map;
+ col_map = input_col_map;
+
+ const unsigned int n_rows = sp.n_rows();
+
+ std::vector<int> n_entries_per_row(n_rows);
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ n_entries_per_row[row] = sp.row_length(row);
+
+ graph = std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(Copy, row_map,
+ n_entries_per_row[input_row_map.MinMyGID()],
+ false));
+
+ Assert (graph->NumGlobalRows() == (int)sp.n_rows(),
+ ExcDimensionMismatch (graph->NumGlobalRows(),
+ sp.n_rows()));
+
+
+ std::vector<int> row_indices;
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ if (row_map.MyGID(row))
+ {
+ const int row_length = sp.row_length(row);
+ row_indices.resize (row_length, -1);
+
+ CompressedSetSparsityPattern::row_iterator col_num =
+ sp.row_begin (row);
+
+ for (unsigned int col = 0;
+ col_num != sp.row_end (row);
+ ++col_num, ++col)
+ row_indices[col] = *col_num;
+
+ graph->Epetra_CrsGraph::InsertGlobalIndices (row, row_length,
+ &row_indices[0]);
+ }
+
+ compress();
+ }
+
+
+
+ /* void
+ SparsityPattern::copy_from (const SparsityPattern &sp)
+ {
+ graph.reset();
+ row_map = sp.row_map;
+ col_map = sp.col_map;
+
+ graph = std::auto_ptr<Epetra_FECrsGraph> (new Epetra_FECrsGraph(*sp.graph));
+ }
+ */
+
+
+ template <typename SparsityType>
+ void
+ SparsityPattern::copy_from (const SparsityType &sp)
+ {
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ Epetra_MpiComm trilinos_communicator (MPI_COMM_WORLD);
+#else
+ Epetra_SerialComm trilinos_communicator;
+#endif
+
+ const Epetra_Map rows (sp.n_rows(), 0, trilinos_communicator);
+ const Epetra_Map columns (sp.n_cols(), 0, trilinos_communicator);
+
+ reinit (rows, columns, sp);
+ }
+
+
+
+ void
+ SparsityPattern::clear ()
+ {
+ // When we clear the matrix, reset
+ // the pointer and generate an
+ // empty matrix.
+ graph.reset();
+
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ row_map = Epetra_Map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD));
+#else
+ row_map = Epetra_Map (0, 0, Epetra_SerialComm());
+#endif
+
+ col_map = row_map;
+
+ graph = std::auto_ptr<Epetra_FECrsGraph>
+ (new Epetra_FECrsGraph(View, row_map, 0));
+
+ graph->FillComplete();
+
+ compressed = true;
+ }
+
+
+
+ void
+ SparsityPattern::compress ()
+ {
+ // flush buffers
+ int ierr;
+ ierr = graph->GlobalAssemble (col_map, row_map, true);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ ierr = graph->OptimizeStorage ();
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ compressed = true;
+ }
+
+
+
+ bool
+ SparsityPattern::exists (const unsigned int i,
+ const unsigned int j) const
+ {
+ // Extract local indices in
+ // the matrix.
+ int trilinos_i = graph->LRID(i), trilinos_j = graph->LRID(j);
+
+ // If the data is not on the
+ // present processor, we throw
+ // an exception. This is on of
+ // the two tiny differences to
+ // the el(i,j) call, which does
+ // not throw any assertions.
+ if (trilinos_i == -1)
+ {
+ return false;
+ }
+ else
+ {
+ // Check whether the matrix
+ // already is transformed to
+ // local indices.
+ if (graph->Filled() == false)
+ {
+ int nnz_present = graph->NumGlobalIndices(i);
+ int nnz_extracted;
+ int *col_indices;
+
+ // Generate the view and make
+ // sure that we have not generated
+ // an error.
+ int ierr = graph->ExtractGlobalRowView(trilinos_i, nnz_extracted,
+ col_indices);
+ Assert (ierr==0, ExcTrilinosError(ierr));
+ Assert (nnz_present == nnz_extracted,
+ ExcDimensionMismatch(nnz_present, nnz_extracted));
+
+ // Search the index
+ int* el_find = std::find(col_indices, col_indices + nnz_present,
+ trilinos_j);
+
+ int local_col_index = (int)(el_find - col_indices);
+
+ if (local_col_index == nnz_present)
+ return false;
+ }
+ else
+ {
+ // Prepare pointers for extraction
+ // of a view of the row.
+ int nnz_present = graph->NumGlobalIndices(i);
+ int nnz_extracted;
+ int *col_indices;
+
+ // Generate the view and make
+ // sure that we have not generated
+ // an error.
+ int ierr = graph->ExtractMyRowView(trilinos_i, nnz_extracted,
+ col_indices);
+ Assert (ierr==0, ExcTrilinosError(ierr));
+
+ Assert (nnz_present == nnz_extracted,
+ ExcDimensionMismatch(nnz_present, nnz_extracted));
+
+ // Search the index
+ int* el_find = std::find(col_indices, col_indices + nnz_present,
+ trilinos_j);
+
+ int local_col_index = (int)(el_find - col_indices);
+
+ if (local_col_index == nnz_present)
+ return false;
+ }
+ }
+
+ return true;
+ }
+
+
+
+ unsigned int
+ SparsityPattern::n_rows () const
+ {
+ int n_rows = graph -> NumGlobalRows();
+
+ return n_rows;
+ }
+
+
+
+ unsigned int
+ SparsityPattern::n_cols () const
+ {
+ int n_cols = graph -> NumGlobalCols();
+ return n_cols;
+ }
+
+
+
+ unsigned int
+ SparsityPattern::local_size () const
+ {
+ int n_rows = graph -> NumMyRows();
+
+ return n_rows;
+ }
+
+
+
+ std::pair<unsigned int, unsigned int>
+ SparsityPattern::local_range () const
+ {
+ unsigned int begin, end;
+ begin = graph -> RowMap().MinMyGID();
+ end = graph -> RowMap().MaxMyGID()+1;
+
+ return std::make_pair (begin, end);
+ }
+
+
+
+ unsigned int
+ SparsityPattern::n_nonzero_elements () const
+ {
+ int nnz = graph->NumGlobalEntries();
+
+ return static_cast<unsigned int>(nnz);
+ }
+
+
+
+ unsigned int
+ SparsityPattern::max_entries_per_row () const
+ {
+ int nnz = graph->MaxRowDim();
+
+ return static_cast<unsigned int>(nnz);
+ }
+
+
+
+ unsigned int
+ SparsityPattern::row_length (const unsigned int row) const
+ {
+ Assert (row < n_rows(), ExcInternalError());
+
+ // get a representation of the
+ // present row
+ int ncols = -1;
+ int local_row = graph->LRID(row);
+
+ // on the processor who owns this
+ // row, we'll have a non-negative
+ // value.
+ if (local_row >= 0)
+ ncols = graph->NumMyIndices (local_row);
+
+ return static_cast<unsigned int>(ncols);
+ }
+
+
+
+ void
+ SparsityPattern::write_ascii ()
+ {
+ Assert (false, ExcNotImplemented());
+ }
+
+
+
+ // As of now, no particularly neat
+ // ouput is generated in case of
+ // multiple processors.
+ void
+ SparsityPattern::print (std::ostream &out) const
+ {
+ int * indices;
+ int num_entries;
+
+ for (int i=0; i<graph->NumMyRows(); ++i)
+ {
+ graph->ExtractMyRowView (i, num_entries, indices);
+ for (int j=0; j<num_entries; ++j)
+ out << "(" << i << "," << indices[graph->GRID(j)] << ") "
+ << std::endl;
+ }
+
+ AssertThrow (out, ExcIO());
+ }
+
+
+
+
+ // explicit instantiations
+ //
+ template void
+ SparsityPattern::copy_from (const dealii::SparsityPattern &);
+ template void
+ SparsityPattern::copy_from (const CompressedSparsityPattern &);
+ template void
+ SparsityPattern::copy_from (const CompressedSetSparsityPattern &);
+ template void
+ SparsityPattern::copy_from (const CompressedSimpleSparsityPattern &);
+
+
+ template void
+ SparsityPattern::reinit (const Epetra_Map &,
+ const dealii::SparsityPattern &);
+ template void
+ SparsityPattern::reinit (const Epetra_Map &,
+ const CompressedSparsityPattern &);
+ template void
+ SparsityPattern::reinit (const Epetra_Map &,
+ const CompressedSetSparsityPattern &);
+ template void
+ SparsityPattern::reinit (const Epetra_Map &,
+ const CompressedSimpleSparsityPattern &);
+
+
+ template void
+ SparsityPattern::reinit (const Epetra_Map &,
+ const Epetra_Map &,
+ const dealii::SparsityPattern &);
+ template void
+ SparsityPattern::reinit (const Epetra_Map &,
+ const Epetra_Map &,
+ const CompressedSparsityPattern &);
+ template void
+ SparsityPattern::reinit (const Epetra_Map &,
+ const Epetra_Map &,
+ const CompressedSimpleSparsityPattern &);
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS