Fix DyanmicSparsityPattern::iterator.

In my previous commit where they were introduced, I had failed to
take into account that DynamicSparsityPattern can take an IndexSet
that denotes which rows to store. We need to skip over rows
that are not stored by this sparsity pattern.

This, unfortunately, complicated the logic significantly.

diff --git a/include/deal.II/lac/dynamic_sparsity_pattern.h b/include/deal.II/lac/dynamic_sparsity_pattern.h
index 388622ffaae006ec2c18b343a15c82b583d70c92..a7b76e846c887ba974e8e9908843ffc6dac486cc 100644 (file)
--- a/include/deal.II/lac/dynamic_sparsity_pattern.h
+++ b/include/deal.II/lac/dynamic_sparsity_pattern.h
@@ -476,6 +476,12 @@ public:
    *
    * Note the discussion in the general documentation of this class about the
    * order in which elements are accessed.
+   *
+   * @note If the sparsity pattern has been initialized with an IndexSet
+   * that denotes which rows to store, then iterators will simply skip
+   * over rows that are not stored. In other words, they will look like
+   * empty rows, but no exception will be generated when iterating over
+   * such rows.
    */
   iterator begin () const;
 
@@ -494,6 +500,12 @@ public:
    *
    * Note also the discussion in the general documentation of this class about
    * the order in which elements are accessed.
+   *
+   * @note If the sparsity pattern has been initialized with an IndexSet
+   * that denotes which rows to store, then iterators will simply skip
+   * over rows that are not stored. In other words, they will look like
+   * empty rows, but no exception will be generated when iterating over
+   * such rows.
    */
   iterator begin (const size_type r) const;
 
@@ -640,12 +652,35 @@ namespace DynamicSparsityPatternIterators
     :
     sparsity_pattern(sparsity_pattern),
     current_row (row),
-    current_entry(sparsity_pattern->lines[current_row].entries.begin()+index_within_row)
+    current_entry(((sparsity_pattern->rowset.size()==0)
+                   ?
+                   sparsity_pattern->lines[current_row].entries.begin()
+                   :
+                   sparsity_pattern->lines[sparsity_pattern->rowset.index_within_set(current_row)].entries.begin())
+                  +
+                  index_within_row)
   {
     Assert (current_row < sparsity_pattern->n_rows(),
             ExcIndexRange (row, 0, sparsity_pattern->n_rows()));
-    Assert (index_within_row < sparsity_pattern->lines[current_row].entries.size(),
-            ExcIndexRange (index_within_row, 0, sparsity_pattern->lines[current_row].entries.size()));
+    Assert ((sparsity_pattern->rowset.size()==0)
+            ||
+            sparsity_pattern->rowset.is_element(current_row),
+            ExcMessage ("You can't create an iterator into a "
+                        "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())));
   }
 
 
@@ -689,7 +724,12 @@ namespace DynamicSparsityPatternIterators
     Assert (current_row < sparsity_pattern->n_rows(),
             ExcInternalError());
 
-    return current_entry - sparsity_pattern->lines[current_row].entries.begin();
+    return (current_entry -
+            ((sparsity_pattern->rowset.size()==0)
+             ?
+             sparsity_pattern->lines[current_row].entries.begin()
+             :
+             sparsity_pattern->lines[sparsity_pattern->rowset.index_within_set(current_row)].entries.begin()));
   }
 
 
@@ -732,19 +772,23 @@ namespace DynamicSparsityPatternIterators
     // if this moves us beyond the end of the row, go to the next row
     // if possible, or set the iterator to an invalid state if not.
     //
-    // use a while loop to ensure that we properly skip over lines
-    // that are empty
-    while (current_entry ==
-           sparsity_pattern->lines[current_row].entries.end())
+    // going to the next row is a bit complicated because we may have
+    // to skip over empty rows, and because we also have to avoid rows
+    // that aren't listed in a possibly passed IndexSet argument of
+    // the sparsity pattern. consequently, rather than trying to
+    // duplicate code here, just call the begin() function of the
+    // sparsity pattern itself
+    if (current_entry ==
+        ((sparsity_pattern->rowset.size()==0)
+         ?
+         sparsity_pattern->lines[current_row].entries.end()
+         :
+         sparsity_pattern->lines[sparsity_pattern->rowset.index_within_set(current_row)].entries.end()))
       {
-        ++current_row;
-        current_entry = sparsity_pattern->lines[current_row].entries.begin();
-
-        if (current_row == sparsity_pattern->n_rows())
-          {
-            *this = Accessor(sparsity_pattern);  // invalid object
-            return;
-          }
+        if (current_row+1 < sparsity_pattern->n_rows())
+          *this = *sparsity_pattern->begin(current_row+1);
+        else
+          *this = Accessor(sparsity_pattern);  // invalid object
       }
   }
 
@@ -990,9 +1034,15 @@ DynamicSparsityPattern::begin (const size_type r) const
   Assert (r<n_rows(), ExcIndexRangeType<size_type>(r,0,n_rows()));
 
   // find the first row starting at r that has entries and return the
-  // begin iterator to it
+  // begin iterator to it. also skip rows for which we do not have
+  // store anything based on the IndexSet given to the sparsity
+  // pattern
+  //
+  // note: row_length(row) returns zero if the row is not locally stored
   unsigned int row = r;
-  while ((row<n_rows()) && (row_length(row)==0))
+  while ((row<n_rows())
+         &&
+         (row_length(row)==0))
     ++row;
 
   if (row == n_rows())
@@ -1010,15 +1060,21 @@ DynamicSparsityPattern::end (const size_type r) const
   Assert (r<n_rows(), ExcIndexRangeType<size_type>(r,0,n_rows()));
 
   // find the first row after r that has entries and return the begin
-  // iterator to it
+  // iterator to it. also skip rows for which we do not have
+  // store anything based on the IndexSet given to the sparsity
+  // pattern
   unsigned int row = r+1;
-  while ((row<n_rows()) && (row_length(row)==0))
-    ++row;
-
-  if (row == n_rows())
-    return iterator(this);
-  else
-    return iterator(this, row, 0);
+  while ((row<n_rows())
+         &&
+         ((row_length(row)==0)
+          ||
+          ((rowset.size() != 0) &&
+           (rowset.is_element(row) == false))))
+
+    if (row == n_rows())
+      return iterator(this);
+    else
+      return iterator(this, row, 0);
 }
 
 

diff --git a/tests/bits/dynamic_sparsity_pattern_iterator_03.cc b/tests/bits/dynamic_sparsity_pattern_iterator_03.cc
new file mode 100644 (file)
index 0000000..f3d67e8
--- /dev/null
+++ b/tests/bits/dynamic_sparsity_pattern_iterator_03.cc
@@ -0,0 +1,87 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2004 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test DynamicSparsityPattern::iterator with sparsity patterns that
+// have an associated IndexSet
+
+#include "../tests.h"
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test ()
+{
+  IndexSet rows(5);
+  rows.add_index(1);
+  rows.add_index(2);
+  rows.add_index(4);
+  DynamicSparsityPattern sp (5,5,rows);
+  sp.add (1,1);
+  sp.add (2,2);
+  sp.add (4,4);
+  sp.compress ();
+
+  DynamicSparsityPattern::const_iterator i = sp.begin();
+  for (; i!=sp.end(); ++i)
+    deallog << i->row() << ' ' << i->column() << std::endl;
+
+  deallog << "OK" << std::endl;
+
+  i = sp.begin(1);
+  deallog << i->row() << ' ' << i->column() << std::endl;
+  deallog << "OK" << std::endl;  
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test ();
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Exception on processing: " << std::endl
+              << exc.what() << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Unknown exception!" << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    };
+}

diff --git a/tests/bits/dynamic_sparsity_pattern_iterator_03.output b/tests/bits/dynamic_sparsity_pattern_iterator_03.output
new file mode 100644 (file)
index 0000000..91daa94
--- /dev/null
+++ b/tests/bits/dynamic_sparsity_pattern_iterator_03.output
@@ -0,0 +1,7 @@
+
+DEAL::1 1
+DEAL::2 2
+DEAL::4 4
+DEAL::OK
+DEAL::1 1
+DEAL::OK