template <typename MatrixType> class BlockMatrixBase;
template <typename number> class SparseMatrix;
+class SparsityPattern;
namespace TrilinosWrappers
{
template <typename number>
void reinit (const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
const double drop_tolerance=1e-13,
- const bool copy_values=true);
+ const bool copy_values=true,
+ const ::dealii::SparsityPattern *use_this_sparsity=0);
/**
* This reinit function takes as
*/
void reinit (const Epetra_CrsMatrix &input_matrix,
const bool copy_values = true);
-
- /**
- * 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.
- *
- * This is a collective operation
- * that needs to be called on all
- * processors in order to avoid a
- * dead lock.
- */
- void clear ();
-
- /**
- * This command does two things:
- * <ul>
- * <li> If the matrix was initialized
- * without a sparsity pattern,
- * elements have been added manually
- * using the set() command. When this
- * process is completed, a call to
- * compress() reorganizes the
- * internal data structures (aparsity
- * pattern) so that a fast access to
- * data is possible in matrix-vector
- * products.
- * <li> If the matrix structure has
- * already been fixed (either by
- * initialization with a sparsity
- * pattern or by calling compress()
- * during the setup phase), this
- * command does the %parallel
- * exchange of data. This is
- * necessary when we perform assembly
- * on more than one (MPI) process,
- * because then some non-local row
- * data will accumulate on nodes that
- * belong to the current's processor
- * element, but are actually held by
- * another. This command is usually
- * called after all elements have
- * been traversed.
- * </ul>
- *
- * In both cases, this function
- * compresses the data structures and
- * allows the resulting matrix to be
- * used in all other operations like
- * matrix-vector products. This is a
- * collective operation, i.e., it
- * needs to be run on all processors
- * when used in %parallel.
- *
- * See @ref GlossCompress "Compressing distributed objects"
- * for more information.
- */
- void compress ();
-
- /**
- * Returns the state of the matrix,
- * i.e., whether compress() needs to
- * be called after an operation
- * requiring data exchange. A call to
- * compress() is also needed when the
- * method set() has been called (even
- * when working in serial).
- */
- bool is_compressed () const;
//@}
/**
* @name Constructors and initialization using an Epetra_Map description
void reinit (const Epetra_Map ¶llel_partitioning,
const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
const double drop_tolerance=1e-13,
- const bool copy_values=true);
+ const bool copy_values=true,
+ const ::dealii::SparsityPattern *use_this_sparsity=0);
/**
* This function is similar to the
const Epetra_Map &col_parallel_partitioning,
const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
const double drop_tolerance=1e-13,
- const bool copy_values=true);
+ const bool copy_values=true,
+ const ::dealii::SparsityPattern *use_this_sparsity=0);
//@}
/**
* @name Constructors and initialization using an IndexSet description
const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
const MPI_Comm &communicator = MPI_COMM_WORLD,
const double drop_tolerance=1e-13,
- const bool copy_values=true);
+ const bool copy_values=true,
+ const ::dealii::SparsityPattern *use_this_sparsity=0);
/**
* This function is similar to the
const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
const MPI_Comm &communicator = MPI_COMM_WORLD,
const double drop_tolerance=1e-13,
- const bool copy_values=true);
+ const bool copy_values=true,
+ const ::dealii::SparsityPattern *use_this_sparsity=0);
//@}
/**
* @name Information on the matrix
*/
unsigned int row_length (const unsigned int row) const;
+ /**
+ * Returns the state of the matrix,
+ * i.e., whether compress() needs to
+ * be called after an operation
+ * requiring data exchange. A call to
+ * compress() is also needed when the
+ * method set() has been called (even
+ * when working in serial).
+ */
+ bool is_compressed () const;
+
/**
* Determine an estimate for the memory
* consumption (in bytes) of this
* @name Modifying entries
*/
//@{
+
+ /**
+ * 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.
+ *
+ * This is a collective operation
+ * that needs to be called on all
+ * processors in order to avoid a
+ * dead lock.
+ */
+ void clear ();
+
+ /**
+ * This command does two things:
+ * <ul>
+ * <li> If the matrix was initialized
+ * without a sparsity pattern,
+ * elements have been added manually
+ * using the set() command. When this
+ * process is completed, a call to
+ * compress() reorganizes the
+ * internal data structures (aparsity
+ * pattern) so that a fast access to
+ * data is possible in matrix-vector
+ * products.
+ * <li> If the matrix structure has
+ * already been fixed (either by
+ * initialization with a sparsity
+ * pattern or by calling compress()
+ * during the setup phase), this
+ * command does the %parallel
+ * exchange of data. This is
+ * necessary when we perform assembly
+ * on more than one (MPI) process,
+ * because then some non-local row
+ * data will accumulate on nodes that
+ * belong to the current's processor
+ * element, but are actually held by
+ * another. This command is usually
+ * called after all elements have
+ * been traversed.
+ * </ul>
+ *
+ * In both cases, this function
+ * compresses the data structures and
+ * allows the resulting matrix to be
+ * used in all other operations like
+ * matrix-vector products. This is a
+ * collective operation, i.e., it
+ * needs to be run on all processors
+ * when used in %parallel.
+ *
+ * See @ref GlossCompress "Compressing distributed objects"
+ * for more information.
+ */
+ void compress ();
+
/**
* Set the element (<i>i,j</i>)
* to @p value.
const ::dealii::SparseMatrix<number> &sparse_matrix,
const MPI_Comm &communicator,
const double drop_tolerance,
- const bool copy_values)
+ const bool copy_values,
+ const ::dealii::SparsityPattern *use_this_sparsity)
{
Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator, false);
- reinit (map, map, sparse_matrix, drop_tolerance, copy_values);
+ reinit (map, map, sparse_matrix, drop_tolerance, copy_values,
+ use_this_sparsity);
}
const ::dealii::SparseMatrix<number> &sparse_matrix,
const MPI_Comm &communicator,
const double drop_tolerance,
- const bool copy_values)
+ const bool copy_values,
+ const ::dealii::SparsityPattern *use_this_sparsity)
{
Epetra_Map row_map =
row_parallel_partitioning.make_trilinos_map (communicator, false);
Epetra_Map col_map =
col_parallel_partitioning.make_trilinos_map (communicator, false);
- reinit (row_map, col_map, sparse_matrix, drop_tolerance, copy_values);
+ reinit (row_map, col_map, sparse_matrix, drop_tolerance, copy_values,
+ use_this_sparsity);
}
const int row_length = sparsity_pattern.row_length(row);
row_indices.resize (row_length, -1);
- for (int col=0; col < row_length; ++col)
- row_indices[col] = sparsity_pattern.column_number (row, col);
-
- graph->InsertGlobalIndices (static_cast<int>(row),
- row_length, &row_indices[0]);
+ typename SparsityType::row_iterator col_num =
+ sparsity_pattern.row_begin (row),
+ row_end = sparsity_pattern.row_end(row);
+ for (unsigned int col = 0; col_num != row_end; ++col_num, ++col)
+ row_indices[col] = *col_num;
+ graph->Epetra_CrsGraph::InsertGlobalIndices (row, row_length,
+ &row_indices[0]);
}
// Eventually, optimize the graph
- // 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,
- const bool exchange_data)
- {
- // this function is similar to the other
- // reinit function with sparsity pattern
- // argument. See there for additional
- // comments.
- matrix.reset();
-
- const unsigned int n_rows = sparsity_pattern.n_rows();
-
- Assert (exchange_data == false, ExcInternalError());
- 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()));
- }
-
- column_space_map = std::auto_ptr<Epetra_Map> (new Epetra_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);
-
- std::auto_ptr<Epetra_CrsGraph> graph;
- if (input_row_map.Comm().NumProc() > 1)
- graph = std::auto_ptr<Epetra_CrsGraph>
- (new Epetra_CrsGraph (Copy, input_row_map,
- &n_entries_per_row[input_row_map.MinMyGID()], true));
- else
- graph = std::auto_ptr<Epetra_CrsGraph>
- (new Epetra_CrsGraph (Copy, input_row_map, input_col_map,
- &n_entries_per_row[input_row_map.MinMyGID()], true));
-
- Assert (graph->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
- ExcDimensionMismatch (graph->NumGlobalRows(),
- sparsity_pattern.n_rows()));
-
- std::vector<int> row_indices;
-
- for (unsigned int row=0; row<n_rows; ++row)
- if (input_row_map.MyGID(row))
- {
- const int row_length = sparsity_pattern.row_length(row);
- row_indices.resize (row_length, -1);
-
- 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;
-
- graph->InsertGlobalIndices (row, row_length, &row_indices[0]);
- }
-
- graph->FillComplete(input_col_map, input_row_map);
- graph->OptimizeStorage();
-
- Assert (graph->NumGlobalCols() == (int)sparsity_pattern.n_cols(),
- ExcDimensionMismatch (graph->NumGlobalCols(),
- sparsity_pattern.n_cols()));
-
- matrix = std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, *graph, false));
-
- last_action = Zero;
- compress();
- }
-
-
-
void
SparseMatrix::reinit (const SparsityPattern &sparsity_pattern)
{
void
SparseMatrix::reinit (const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
const double drop_tolerance,
- const bool copy_values)
+ const bool copy_values,
+ const ::dealii::SparsityPattern *use_this_sparsity)
{
const Epetra_Map rows (dealii_sparse_matrix.m(),
0,
const Epetra_Map columns (dealii_sparse_matrix.n(),
0,
Utilities::Trilinos::comm_self());
- reinit (rows, columns, dealii_sparse_matrix, drop_tolerance, copy_values);
+ reinit (rows, columns, dealii_sparse_matrix, drop_tolerance,
+ copy_values, use_this_sparsity);
}
SparseMatrix::reinit (const Epetra_Map &input_map,
const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
const double drop_tolerance,
- const bool copy_values)
+ const bool copy_values,
+ const ::dealii::SparsityPattern *use_this_sparsity)
{
reinit (input_map, input_map, dealii_sparse_matrix, drop_tolerance,
- copy_values);
+ copy_values, use_this_sparsity);
}
const Epetra_Map &input_col_map,
const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
const double drop_tolerance,
- const bool copy_values)
+ const bool copy_values,
+ const ::dealii::SparsityPattern *use_this_sparsity)
{
matrix.reset();
{
// in case we do not copy values, just
// call the other function.
- reinit (input_row_map, input_col_map,
- dealii_sparse_matrix.get_sparsity_pattern());
+ if (use_this_sparsity == 0)
+ reinit (input_row_map, input_col_map,
+ dealii_sparse_matrix.get_sparsity_pattern());
+ else
+ reinit (input_row_map, input_col_map,
+ *use_this_sparsity);
return;
}
column_space_map = std::auto_ptr<Epetra_Map> (new Epetra_Map (input_col_map));
+ const ::dealii::SparsityPattern & sparsity_pattern =
+ (use_this_sparsity!=0)? *use_this_sparsity :
+ dealii_sparse_matrix.get_sparsity_pattern();
+
std::vector<int> n_entries_per_row(n_rows);
+ unsigned int maximum_row_length = 0;
for (unsigned int row=0; row<n_rows; ++row)
- n_entries_per_row[(int)row] =
- dealii_sparse_matrix.get_sparsity_pattern().row_length(row);
+ {
+ const unsigned int row_length = sparsity_pattern.row_length(row);
+ if (row_length > maximum_row_length)
+ maximum_row_length = row_length;
+ n_entries_per_row[row] = row_length;
+ }
- matrix = std::auto_ptr<Epetra_FECrsMatrix>
- (new Epetra_FECrsMatrix(Copy, input_row_map,
- &n_entries_per_row[input_row_map.MinMyGID()],
- false));
+ std::auto_ptr<Epetra_CrsGraph> graph;
+ if (input_row_map.Comm().NumProc() > 1)
+ graph = std::auto_ptr<Epetra_CrsGraph>
+ (new Epetra_CrsGraph (Copy, input_row_map,
+ &n_entries_per_row[input_row_map.MinMyGID()], false));
+ else
+ graph = std::auto_ptr<Epetra_CrsGraph>
+ (new Epetra_CrsGraph (Copy, input_row_map, input_col_map,
+ &n_entries_per_row[input_row_map.MinMyGID()], false));
- std::vector<TrilinosScalar> values;
- std::vector<unsigned int> row_indices;
+ Assert (graph->NumGlobalRows() == static_cast<int>(n_rows),
+ ExcDimensionMismatch (graph->NumGlobalRows(), n_rows));
+ // now insert the indices
+ std::vector<int> row_indices (maximum_row_length);
for (unsigned int row=0; row<n_rows; ++row)
- {
- values.resize (n_entries_per_row[row],0.);
- row_indices.resize (n_entries_per_row[row],
- numbers::invalid_unsigned_int);
-
- unsigned int index = 0;
- for (typename ::dealii::SparseMatrix<number>::const_iterator
- p = dealii_sparse_matrix.begin(row);
- p != dealii_sparse_matrix.end(row); ++p)
- if (std::fabs(p->value()) > drop_tolerance)
+ if (input_row_map.MyGID(row))
+ {
+ typename ::dealii::SparsityPattern::row_iterator
+ col_num = sparsity_pattern.row_begin (row),
+ row_end = sparsity_pattern.row_end(row);
+ typename ::dealii::SparseMatrix<number>::const_iterator
+ matrix_col_num = dealii_sparse_matrix.begin(row),
+ matrix_row_end = dealii_sparse_matrix.end(row);
+
+ unsigned int col = 0;
+ if (sparsity_pattern.optimize_diagonal())
+ {
+ row_indices[col++] = row;
+ ++col_num;
+ ++ matrix_col_num;
+ }
+ while (col_num != row_end && matrix_col_num != matrix_row_end)
{
- row_indices[index] = p->column();
- values[index] = p->value();
- ++index;
+ while (*col_num < matrix_col_num->column() && col_num != row_end)
+ ++col_num;
+ while (*col_num > matrix_col_num->column() && matrix_col_num != matrix_row_end)
+ ++matrix_col_num;
+
+ if (std::fabs(matrix_col_num->value()) > drop_tolerance)
+ row_indices[col++] = *col_num;
+ ++col_num;
+ ++matrix_col_num;
}
- set (row, index, &row_indices[0], &values[0], false);
- }
+ graph->Epetra_CrsGraph::InsertGlobalIndices (row, col,
+ &row_indices[0]);
+ }
+ graph->FillComplete(input_col_map, input_row_map);
+ graph->OptimizeStorage();
+
+ matrix = std::auto_ptr<Epetra_FECrsMatrix>
+ (new Epetra_FECrsMatrix(Copy, *graph, false));
+ compress();
+
+ std::vector<TrilinosScalar> values (maximum_row_length);
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ if (input_row_map.MyGID(row))
+ {
+ typename ::dealii::SparsityPattern::row_iterator
+ col_num = sparsity_pattern.row_begin (row),
+ row_end = sparsity_pattern.row_end(row);
+ typename ::dealii::SparseMatrix<number>::const_iterator
+ matrix_col_num = dealii_sparse_matrix.begin(row),
+ matrix_row_end = dealii_sparse_matrix.end(row);
+
+ unsigned int col = 0;
+ if (sparsity_pattern.optimize_diagonal())
+ {
+ values[col] = matrix_col_num->value();
+ row_indices[col++] = row;
+ ++col_num;
+ ++ matrix_col_num;
+ }
+ while (col_num != row_end && matrix_col_num != matrix_row_end)
+ {
+ while (*col_num < matrix_col_num->column() && col_num != row_end)
+ ++col_num;
+ while (*col_num > matrix_col_num->column() && matrix_col_num != matrix_row_end)
+ ++matrix_col_num;
+
+ if (std::fabs(matrix_col_num->value()) > drop_tolerance)
+ {
+ values[col] = matrix_col_num->value();
+ row_indices[col++] = *col_num;
+ }
+ ++col_num;
+ ++matrix_col_num;
+ }
+ set (row, col, reinterpret_cast<unsigned int*>(&row_indices[0]),
+ &values[0], false);
+ }
compress();
}
const Epetra_Map &,
const CompressedSparsityPattern &,
const bool);
+ template void
+ SparseMatrix::reinit (const Epetra_Map &,
+ const Epetra_Map &,
+ const CompressedSimpleSparsityPattern &,
+ const bool);
+ template void
+ SparseMatrix::reinit (const Epetra_Map &,
+ const Epetra_Map &,
+ const CompressedSetSparsityPattern &,
+ const bool);
template void
SparseMatrix::reinit (const dealii::SparseMatrix<float> &,
const double,
- const bool);
+ const bool,
+ const dealii::SparsityPattern *);
template void
SparseMatrix::reinit (const dealii::SparseMatrix<double> &,
const double,
- const bool);
+ const bool,
+ const dealii::SparsityPattern *);
template void
SparseMatrix::reinit (const dealii::SparseMatrix<long double> &,
const double,
- const bool);
+ const bool,
+ const dealii::SparsityPattern *);
template void
SparseMatrix::reinit (const Epetra_Map &,
const dealii::SparseMatrix<float> &,
const double,
- const bool);
+ const bool,
+ const dealii::SparsityPattern *);
template void
SparseMatrix::reinit (const Epetra_Map &,
const dealii::SparseMatrix<double> &,
const double,
- const bool);
+ const bool,
+ const dealii::SparsityPattern *);
template void
SparseMatrix::reinit (const Epetra_Map &,
const dealii::SparseMatrix<long double> &,
const double,
- const bool);
+ const bool,
+ const dealii::SparsityPattern *);
template void
SparseMatrix::reinit (const Epetra_Map &,
const Epetra_Map &,
const dealii::SparseMatrix<float> &,
const double,
- const bool);
+ const bool,
+ const dealii::SparsityPattern *);
template void
SparseMatrix::reinit (const Epetra_Map &,
const Epetra_Map &,
const dealii::SparseMatrix<double> &,
const double,
- const bool);
+ const bool,
+ const dealii::SparsityPattern *);
template void
SparseMatrix::reinit (const Epetra_Map &,
const Epetra_Map &,
const dealii::SparseMatrix<long double> &,
const double,
- const bool);
+ const bool,
+ const dealii::SparsityPattern *);
}