From: Martin Kronbichler Date: Mon, 21 Sep 2015 17:26:17 +0000 (+0200) Subject: Avoid using DynamicSparsityPattern::iterator past the end of the row X-Git-Tag: v8.4.0-rc2~384^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e0d8ca889ed6df70a61f9d7bb95e6e1d262864a6;p=dealii.git Avoid using DynamicSparsityPattern::iterator past the end of the row It is way too slow for the parallel case because one needs to go over all indices, not only locally owned ones Also fix several 64 bit bugs. Improve performance of nonlocal graph of Trilinos matrix by not setting dummy entries on each processor --- diff --git a/include/deal.II/lac/dynamic_sparsity_pattern.h b/include/deal.II/lac/dynamic_sparsity_pattern.h index 458c4cba1f..ee85de1820 100644 --- a/include/deal.II/lac/dynamic_sparsity_pattern.h +++ b/include/deal.II/lac/dynamic_sparsity_pattern.h @@ -69,7 +69,7 @@ namespace DynamicSparsityPatternIterators * Constructor. */ Accessor (const DynamicSparsityPattern *sparsity_pattern, - const unsigned int row, + const size_type row, const unsigned int index_within_row); /** @@ -115,7 +115,7 @@ namespace DynamicSparsityPatternIterators /** * The row we currently point into. */ - unsigned int current_row; + size_type current_row; /** * A pointer to the element within the current row that we currently point @@ -177,7 +177,7 @@ namespace DynamicSparsityPatternIterators * the zeroth row). */ Iterator (const DynamicSparsityPattern *sp, - const unsigned int row, + const size_type row, const unsigned int index_within_row); /** @@ -635,7 +635,7 @@ namespace DynamicSparsityPatternIterators inline Accessor:: Accessor (const DynamicSparsityPattern *sparsity_pattern, - const unsigned int row, + const size_type row, const unsigned int index_within_row) : sparsity_pattern(sparsity_pattern), @@ -653,8 +653,7 @@ namespace DynamicSparsityPatternIterators : sparsity_pattern->lines[sparsity_pattern->rowset.index_within_set(current_row)].entries.end()) { - Assert (current_row < sparsity_pattern->n_rows(), - ExcIndexRange (row, 0, sparsity_pattern->n_rows())); + AssertIndexRange(current_row, sparsity_pattern->n_rows()); Assert ((sparsity_pattern->rowset.size()==0) || sparsity_pattern->rowset.is_element(current_row), @@ -662,18 +661,12 @@ namespace DynamicSparsityPatternIterators "DynamicSparsityPattern's row that is not " "actually stored by that sparsity pattern " "based on the IndexSet argument to it.")); - Assert (index_within_row < - ((sparsity_pattern->rowset.size()==0) - ? - sparsity_pattern->lines[current_row].entries.size() - : - sparsity_pattern->lines[sparsity_pattern->rowset.index_within_set(current_row)].entries.size()), - ExcIndexRange (index_within_row, 0, - ((sparsity_pattern->rowset.size()==0) - ? - sparsity_pattern->lines[current_row].entries.size() - : - sparsity_pattern->lines[sparsity_pattern->rowset.index_within_set(current_row)].entries.size()))); + AssertIndexRange(index_within_row, + ((sparsity_pattern->rowset.size()==0) + ? + sparsity_pattern->lines[current_row].entries.size() + : + sparsity_pattern->lines[sparsity_pattern->rowset.index_within_set(current_row)].entries.size())); } @@ -682,7 +675,7 @@ namespace DynamicSparsityPatternIterators Accessor (const DynamicSparsityPattern *sparsity_pattern) : sparsity_pattern(sparsity_pattern), - current_row(numbers::invalid_unsigned_int), + current_row(numbers::invalid_size_type), current_entry(), end_of_row() {} @@ -739,7 +732,7 @@ namespace DynamicSparsityPatternIterators // current_entry field may not point to a deterministic location return (sparsity_pattern == other.sparsity_pattern && current_row == other.current_row && - ((current_row == numbers::invalid_unsigned_int) + ((current_row == numbers::invalid_size_type) || (current_entry == other.current_entry))); } @@ -753,14 +746,14 @@ namespace DynamicSparsityPatternIterators ExcInternalError()); // if *this is past-the-end, then it is less than no one - if (current_row == numbers::invalid_unsigned_int) + if (current_row == numbers::invalid_size_type) return (false); // now *this should be an valid value Assert (current_row < sparsity_pattern->n_rows(), ExcInternalError()); // if other is past-the-end - if (other.current_row == numbers::invalid_unsigned_int) + if (other.current_row == numbers::invalid_size_type) return (true); // now other should be an valid value Assert (other.current_row < sparsity_pattern->n_rows(), @@ -805,7 +798,7 @@ namespace DynamicSparsityPatternIterators inline Iterator::Iterator (const DynamicSparsityPattern *sparsity_pattern, - const unsigned int row, + const size_type row, const unsigned int index_within_row) : accessor(sparsity_pattern, row, index_within_row) @@ -1048,7 +1041,10 @@ DynamicSparsityPattern::begin (const size_type r) const // pattern // // note: row_length(row) returns zero if the row is not locally stored - unsigned int row = r; + // + // TODO: this is way too slow when used in parallel, so do not use it on + // non-owned rows + size_type row = r; while ((row5000 processors because that gets very + // slow. Thus, we set a flag in Epetra_CrsGraph that sets the correct + // flag. Since it is protected, we need to expose this information by + // deriving a class from Epetra_CrsGraph for the purpose of creating the + // data structure + class Epetra_CrsGraphMod : public Epetra_CrsGraph + { + public: + Epetra_CrsGraphMod (const Epetra_Map &row_map, + const int *n_entries_per_row) + : + Epetra_CrsGraph(Copy, row_map, n_entries_per_row, true) + {}; + + void SetIndicesAreGlobal() + { + this->Epetra_CrsGraph::SetIndicesAreGlobal(true); + } + }; + + // specialization for DynamicSparsityPattern which can provide us with // more information about the non-locally owned rows template <> @@ -636,28 +659,20 @@ namespace TrilinosWrappers } } - // make sure all processors create an off-processor matrix with at least - // one entry - if (have_ghost_rows == true && ghost_rows.empty() == true) - { - ghost_rows.push_back(0); - n_entries_per_ghost_row.push_back(1); - } - Epetra_Map off_processor_map(-1, ghost_rows.size(), (ghost_rows.size()>0)?(&ghost_rows[0]):NULL, 0, input_row_map.Comm()); - std_cxx11::shared_ptr graph, nonlocal_graph; + std_cxx11::shared_ptr graph; + std_cxx11::shared_ptr nonlocal_graph; if (input_row_map.Comm().NumProc() > 1) { graph.reset (new Epetra_CrsGraph (Copy, input_row_map, (n_entries_per_row.size()>0)?(&n_entries_per_row[0]):NULL, exchange_data ? false : true)); if (have_ghost_rows == true) - nonlocal_graph.reset (new Epetra_CrsGraph (Copy, off_processor_map, - &n_entries_per_ghost_row[0], - true)); + nonlocal_graph.reset (new Epetra_CrsGraphMod (off_processor_map, + &n_entries_per_ghost_row[0])); } else graph.reset (new Epetra_CrsGraph (Copy, input_row_map, input_col_map, @@ -676,11 +691,8 @@ namespace TrilinosWrappers continue; row_indices.resize (row_length, -1); - { - dealii::DynamicSparsityPattern::iterator p = sparsity_pattern.begin(global_row); - for (size_type col=0; p != sparsity_pattern.end(global_row); ++p, ++col) - row_indices[col] = p->column(); - } + for (int col=0; col < row_length; ++col) + row_indices[col] = sparsity_pattern.column_number(global_row, col); if (input_row_map.MyGID(global_row)) graph->InsertGlobalIndices (global_row, row_length, &row_indices[0]); @@ -695,14 +707,11 @@ namespace TrilinosWrappers // finalize nonlocal graph and create nonlocal matrix if (nonlocal_graph.get() != 0) { - if (nonlocal_graph->IndicesAreGlobal() == false && - nonlocal_graph->RowMap().NumMyElements() > 0) - { - // insert dummy element - TrilinosWrappers::types::int_type row = - nonlocal_graph->RowMap().MyGID(TrilinosWrappers::types::int_type(0)); - nonlocal_graph->InsertGlobalIndices(row, 1, &row); - } + // must make sure the IndicesAreGlobal flag is set on all processors + // because some processors might not call InsertGlobalIndices (and + // we do not want to insert dummy indices on all processors for + // large-scale simulations due to the bad impact on performance) + nonlocal_graph->SetIndicesAreGlobal(); Assert(nonlocal_graph->IndicesAreGlobal() == true, ExcInternalError()); nonlocal_graph->FillComplete(input_col_map, input_row_map); diff --git a/source/lac/trilinos_sparsity_pattern.cc b/source/lac/trilinos_sparsity_pattern.cc index fa50f974d3..9a3ca88d72 100644 --- a/source/lac/trilinos_sparsity_pattern.cc +++ b/source/lac/trilinos_sparsity_pattern.cc @@ -486,8 +486,14 @@ namespace TrilinosWrappers row_indices.resize (row_length, -1); { typename SparsityType::iterator p = sp.begin(row); - for (size_type col=0; p != sp.end(row); ++p, ++col) - row_indices[col] = p->column(); + // avoid incrementing p over the end of the current row because + // it is slow for DynamicSparsityPattern in parallel + for (int col=0; colcolumn(); + if (col < row_length) + ++p; + } } graph->Epetra_CrsGraph::InsertGlobalIndices (row, row_length, &row_indices[0]); @@ -502,8 +508,14 @@ namespace TrilinosWrappers row_indices.resize (row_length, -1); { typename SparsityType::iterator p = sp.begin(row); - for (size_type col=0; p != sp.end(row); ++p, ++col) - row_indices[col] = p->column(); + // avoid incrementing p over the end of the current row because + // it is slow for DynamicSparsityPattern in parallel + for (int col=0; colcolumn(); + if (col < row_length) + ++p; + } } graph->InsertGlobalIndices (1, reinterpret_cast(&row),