From bf2422d88b4f141d6a822617c8e0d184e14e76ac Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 5 Feb 2010 11:53:27 +0000 Subject: [PATCH] More flexibility with initializing a Trilinos sparse matrix from a deal.II sparse matrix. git-svn-id: https://svn.dealii.org/trunk@20503 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/trilinos_precondition.h | 10 +- .../lac/include/lac/trilinos_sparse_matrix.h | 203 +++++++------ deal.II/lac/source/trilinos_precondition.cc | 26 +- deal.II/lac/source/trilinos_sparse_matrix.cc | 279 ++++++++++-------- 4 files changed, 285 insertions(+), 233 deletions(-) diff --git a/deal.II/lac/include/lac/trilinos_precondition.h b/deal.II/lac/include/lac/trilinos_precondition.h index 227f2fa383..8e6df15845 100644 --- a/deal.II/lac/include/lac/trilinos_precondition.h +++ b/deal.II/lac/include/lac/trilinos_precondition.h @@ -49,6 +49,7 @@ DEAL_II_NAMESPACE_OPEN // forward declarations template class SparseMatrix; template class Vector; +class SparsityPattern; /*! @addtogroup TrilinosWrappers *@{ @@ -1465,7 +1466,8 @@ namespace TrilinosWrappers template void initialize (const ::dealii::SparseMatrix &deal_ii_sparse_matrix, const AdditionalData &additional_data = AdditionalData(), - const double drop_tolerance = 1e-13); + const double drop_tolerance = 1e-13, + const ::dealii::SparsityPattern *use_this_sparsity = 0); /** * This function can be used for a @@ -1492,6 +1494,12 @@ namespace TrilinosWrappers */ void reinit (); + /** + * Prints an estimate of the memory + * consumption of this class. + */ + unsigned int memory_consumption () const; + private: /** diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index 6fb4356d7f..2fabc8216e 100644 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -45,6 +45,7 @@ DEAL_II_NAMESPACE_OPEN template class BlockMatrixBase; template class SparseMatrix; +class SparsityPattern; namespace TrilinosWrappers { @@ -474,7 +475,8 @@ namespace TrilinosWrappers template void reinit (const ::dealii::SparseMatrix &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 @@ -485,93 +487,6 @@ namespace TrilinosWrappers */ 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 - * matrix=0, 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: - *
    - *
  • 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. - *
  • 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. - *
- * - * 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 @@ -787,7 +702,8 @@ namespace TrilinosWrappers void reinit (const Epetra_Map ¶llel_partitioning, const ::dealii::SparseMatrix &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 @@ -815,7 +731,8 @@ namespace TrilinosWrappers const Epetra_Map &col_parallel_partitioning, const ::dealii::SparseMatrix &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 @@ -1039,7 +956,8 @@ namespace TrilinosWrappers const ::dealii::SparseMatrix &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 @@ -1068,7 +986,8 @@ namespace TrilinosWrappers const ::dealii::SparseMatrix &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 @@ -1141,6 +1060,17 @@ namespace TrilinosWrappers */ 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 @@ -1156,6 +1086,83 @@ namespace TrilinosWrappers * @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 + * matrix=0, 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: + *
    + *
  • 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. + *
  • 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. + *
+ * + * 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,j) * to @p value. @@ -3091,10 +3098,12 @@ namespace TrilinosWrappers const ::dealii::SparseMatrix &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); } @@ -3106,13 +3115,15 @@ namespace TrilinosWrappers const ::dealii::SparseMatrix &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); } diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index 25c18e22f4..39aeb5217d 100644 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -645,7 +645,8 @@ namespace TrilinosWrappers PreconditionAMG:: initialize (const ::dealii::SparseMatrix &deal_ii_sparse_matrix, const AdditionalData &additional_data, - const double drop_tolerance) + const double drop_tolerance, + const ::dealii::SparsityPattern *use_this_sparsity) { const unsigned int n_rows = deal_ii_sparse_matrix.m(); @@ -658,26 +659,41 @@ namespace TrilinosWrappers Matrix.reset(); Matrix = std_cxx1x::shared_ptr (new SparseMatrix()); - Matrix->reinit (*map, deal_ii_sparse_matrix, drop_tolerance); - Matrix->compress(); + Matrix->reinit (*map, *map, deal_ii_sparse_matrix, drop_tolerance, true, + use_this_sparsity); initialize (*Matrix, additional_data); } + void PreconditionAMG::reinit () { multilevel_operator->ReComputePreconditioner(); } + unsigned int PreconditionAMG::memory_consumption() const + { + unsigned int memory = sizeof(this); + + // todo: find a way to read out ML's data + // sizes + if (Matrix != 0) + memory += Matrix->memory_consumption(); + return memory; + } + + // explicit instantiations template void PreconditionAMG::initialize (const ::dealii::SparseMatrix &, - const AdditionalData &, const double); + const AdditionalData &, const double, + const ::dealii::SparsityPattern*); template void PreconditionAMG::initialize (const ::dealii::SparseMatrix &, - const AdditionalData &, const double); + const AdditionalData &, const double, + const ::dealii::SparsityPattern*); } diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 363c8a1788..cc952399c2 100644 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -465,12 +465,14 @@ namespace TrilinosWrappers 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(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 @@ -514,93 +516,6 @@ namespace TrilinosWrappers - // 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 (new Epetra_Map (input_col_map)); - - std::vector n_entries_per_row(n_rows); - - for (unsigned int row=0; row graph; - if (input_row_map.Comm().NumProc() > 1) - graph = std::auto_ptr - (new Epetra_CrsGraph (Copy, input_row_map, - &n_entries_per_row[input_row_map.MinMyGID()], true)); - else - graph = std::auto_ptr - (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 row_indices; - - for (unsigned int row=0; rowInsertGlobalIndices (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 - (new Epetra_FECrsMatrix(Copy, *graph, false)); - - last_action = Zero; - compress(); - } - - - void SparseMatrix::reinit (const SparsityPattern &sparsity_pattern) { @@ -639,7 +554,8 @@ namespace TrilinosWrappers void SparseMatrix::reinit (const ::dealii::SparseMatrix &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, @@ -647,7 +563,8 @@ namespace TrilinosWrappers 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); } @@ -657,10 +574,11 @@ namespace TrilinosWrappers SparseMatrix::reinit (const Epetra_Map &input_map, const ::dealii::SparseMatrix &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); } @@ -671,7 +589,8 @@ namespace TrilinosWrappers const Epetra_Map &input_col_map, const ::dealii::SparseMatrix &dealii_sparse_matrix, const double drop_tolerance, - const bool copy_values) + const bool copy_values, + const ::dealii::SparsityPattern *use_this_sparsity) { matrix.reset(); @@ -679,8 +598,12 @@ namespace TrilinosWrappers { // 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; } @@ -695,39 +618,114 @@ namespace TrilinosWrappers column_space_map = std::auto_ptr (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 n_entries_per_row(n_rows); + unsigned int maximum_row_length = 0; for (unsigned int row=0; row maximum_row_length) + maximum_row_length = row_length; + n_entries_per_row[row] = row_length; + } - matrix = std::auto_ptr - (new Epetra_FECrsMatrix(Copy, input_row_map, - &n_entries_per_row[input_row_map.MinMyGID()], - false)); + std::auto_ptr graph; + if (input_row_map.Comm().NumProc() > 1) + graph = std::auto_ptr + (new Epetra_CrsGraph (Copy, input_row_map, + &n_entries_per_row[input_row_map.MinMyGID()], false)); + else + graph = std::auto_ptr + (new Epetra_CrsGraph (Copy, input_row_map, input_col_map, + &n_entries_per_row[input_row_map.MinMyGID()], false)); - std::vector values; - std::vector row_indices; + Assert (graph->NumGlobalRows() == static_cast(n_rows), + ExcDimensionMismatch (graph->NumGlobalRows(), n_rows)); + // now insert the indices + std::vector row_indices (maximum_row_length); for (unsigned int row=0; row::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::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 + (new Epetra_FECrsMatrix(Copy, *graph, false)); + compress(); + + std::vector values (maximum_row_length); + + for (unsigned int row=0; row::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(&row_indices[0]), + &values[0], false); + } compress(); } @@ -1467,54 +1465,73 @@ namespace TrilinosWrappers 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 &, const double, - const bool); + const bool, + const dealii::SparsityPattern *); template void SparseMatrix::reinit (const dealii::SparseMatrix &, const double, - const bool); + const bool, + const dealii::SparsityPattern *); template void SparseMatrix::reinit (const dealii::SparseMatrix &, const double, - const bool); + const bool, + const dealii::SparsityPattern *); template void SparseMatrix::reinit (const Epetra_Map &, const dealii::SparseMatrix &, const double, - const bool); + const bool, + const dealii::SparsityPattern *); template void SparseMatrix::reinit (const Epetra_Map &, const dealii::SparseMatrix &, const double, - const bool); + const bool, + const dealii::SparsityPattern *); template void SparseMatrix::reinit (const Epetra_Map &, const dealii::SparseMatrix &, const double, - const bool); + const bool, + const dealii::SparsityPattern *); template void SparseMatrix::reinit (const Epetra_Map &, const Epetra_Map &, const dealii::SparseMatrix &, const double, - const bool); + const bool, + const dealii::SparsityPattern *); template void SparseMatrix::reinit (const Epetra_Map &, const Epetra_Map &, const dealii::SparseMatrix &, const double, - const bool); + const bool, + const dealii::SparsityPattern *); template void SparseMatrix::reinit (const Epetra_Map &, const Epetra_Map &, const dealii::SparseMatrix &, const double, - const bool); + const bool, + const dealii::SparsityPattern *); } -- 2.39.5