]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Avoid using DynamicSparsityPattern::iterator past the end of the row
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 21 Sep 2015 17:26:17 +0000 (19:26 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 22 Sep 2015 15:28:42 +0000 (17:28 +0200)
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

include/deal.II/lac/dynamic_sparsity_pattern.h
source/base/index_set.cc
source/lac/trilinos_sparse_matrix.cc
source/lac/trilinos_sparsity_pattern.cc

index 458c4cba1f18fbce76191cbc23520c35e3a8f1bb..ee85de182081946338907531835048048ceb4097 100644 (file)
@@ -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 ((row<n_rows())
          &&
          (row_length(row)==0))
index f654da01974aac9adf42be934fee1a885d19cce8..bad685f15686fa755494e632349cd34ede82ec99 100644 (file)
@@ -489,7 +489,7 @@ IndexSet::make_trilinos_map (const MPI_Comm &communicator,
 #ifdef DEBUG
   if (!overlapping)
     {
-      const unsigned int n_global_elements
+      const size_type n_global_elements
         = Utilities::MPI::sum (n_elements(), communicator);
       Assert (n_global_elements == size(),
               ExcMessage ("You are trying to create an Epetra_Map object "
index eb76168c7d30645c754b37f4ffae5536101f59fd..f1f6edd5f77a967043503c390deb9c07c0059665 100644 (file)
@@ -565,6 +565,29 @@ namespace TrilinosWrappers
 
 
 
+    // for the non-local graph, we need to circumvent the problem that some
+    // processors will not add into the non-local graph at all: We do not want
+    // to insert dummy elements on >5000 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<Epetra_CrsGraph> graph, nonlocal_graph;
+      std_cxx11::shared_ptr<Epetra_CrsGraph> graph;
+      std_cxx11::shared_ptr<Epetra_CrsGraphMod> 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);
index fa50f974d39c7f2e435d8dafebb2cf42e4ee873f..9a3ca88d72fd44f27aaa35caec3492d10d695dd7 100644 (file)
@@ -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; col<row_length; )
+                {
+                  row_indices[col++] = p->column();
+                  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; col<row_length; )
+                {
+                  row_indices[col++] = p->column();
+                  if (col < row_length)
+                    ++p;
+                }
             }
             graph->InsertGlobalIndices (1,
                                         reinterpret_cast<TrilinosWrappers::types::int_type *>(&row),

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.