]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Ensure that we can run iterators over Trilinos matrices that are not stored locally.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 9 Mar 2015 05:11:44 +0000 (00:11 -0500)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 9 Mar 2015 05:11:44 +0000 (00:11 -0500)
doc/news/changes.h
include/deal.II/lac/trilinos_sparse_matrix.h
source/lac/trilinos_sparse_matrix.cc
tests/mpi/trilinos_sparse_matrix_03.cc [new file with mode: 0644]
tests/mpi/trilinos_sparse_matrix_03.with_trilinos=true.mpirun=2.output [new file with mode: 0644]
tests/mpi/trilinos_sparse_matrix_04.cc [new file with mode: 0644]
tests/mpi/trilinos_sparse_matrix_04.with_trilinos=true.mpirun=2.output [new file with mode: 0644]

index e9db05ecdf2dcbb8cea8d4eecde1b90b4256a37c..67cbd6fd839442d9af0fb84dea2379ec4d969d2a 100644 (file)
@@ -364,6 +364,15 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: Iterating over the elements of a TrilinosWrappers::SparseMatrix
+  object previously led to errors if the matrix was in fact stored in
+  parallel across multiple MPI processes. This is now fixed: rows not
+  stored locally on the processor where you run the iteration simply look
+  like they're empty.
+  <br>
+  (Wolfgang Bangerth, 2015/03/08)
+  </li>
+
   <li> New: There is now a new macro DeclExceptionMsg that allows to
   declare an exception that does not take any run-time arguments
   yet still allows to specify an error message.
index 283a1ec37236dc7202abf97c0995af428cbb6171..10cabcf12e114414ef6dc8f4be687ff9bb12cc37 100644 (file)
@@ -1684,10 +1684,11 @@ namespace TrilinosWrappers
      * elements if you iterate over them.
      *
      * When you iterate over the elements of a parallel matrix, you
-     * will only be able to access the locally owned rows. In that
-     * case, you probably want to call the begin() function that takes
-     * the row as an argument to limit the range of elements to loop
-     * over.
+     * will only be able to access the locally owned rows. (You can
+     * access the other rows as well, but they will look empty.) In
+     * that case, you probably want to call the begin() function that
+     * takes the row as an argument to limit the range of elements to
+     * loop over.
      */
     const_iterator begin () const;
 
@@ -1728,7 +1729,8 @@ namespace TrilinosWrappers
      *
      * @note When you access the elements of a parallel matrix, you
      * can only access the elements of rows that are actually stored
-     * locally. Even then, if another processor has since written
+     * locally. (You can access the other rows as well, but they will
+     * look empty.) Even then, if another processor has since written
      * into, or added to, an element of the matrix that is stored on
      * the current processor, then you will still see the old value of
      * this entry unless you have called compress() between modifying
@@ -2121,7 +2123,9 @@ namespace TrilinosWrappers
 
           while ((accessor.a_row < accessor.matrix->m())
                  &&
-                 (accessor.matrix->row_length(accessor.a_row) == 0))
+                 ((accessor.matrix->in_local_range (accessor.a_row) == false)
+                 ||
+                 (accessor.matrix->row_length(accessor.a_row) == 0)))
             ++accessor.a_row;
 
           accessor.visit_present_row();
@@ -2210,7 +2214,7 @@ namespace TrilinosWrappers
   SparseMatrix::const_iterator
   SparseMatrix::begin() const
   {
-    return const_iterator(this, 0, 0);
+    return begin(0);
   }
 
 
@@ -2229,7 +2233,9 @@ namespace TrilinosWrappers
   SparseMatrix::begin(const size_type r) const
   {
     Assert (r < m(), ExcIndexRange(r, 0, m()));
-    if (row_length(r) > 0)
+    if (in_local_range (r)
+       &&
+       (row_length(r) > 0))
       return const_iterator(this, r, 0);
     else
       return end (r);
@@ -2247,7 +2253,9 @@ namespace TrilinosWrappers
     // past this line, or at the end of the
     // matrix
     for (size_type i=r+1; i<m(); ++i)
-      if (row_length(i) > 0)
+      if (in_local_range (i)
+         &&
+         (row_length(i) > 0))
         return const_iterator(this, i, 0);
 
     // if there is no such line, then take the
@@ -2261,7 +2269,7 @@ namespace TrilinosWrappers
   SparseMatrix::iterator
   SparseMatrix::begin()
   {
-    return iterator(this, 0, 0);
+    return begin(0);
   }
 
 
@@ -2280,7 +2288,9 @@ namespace TrilinosWrappers
   SparseMatrix::begin(const size_type r)
   {
     Assert (r < m(), ExcIndexRange(r, 0, m()));
-    if (row_length(r) > 0)
+    if (in_local_range (r)
+       &&
+       (row_length(r) > 0))
       return iterator(this, r, 0);
     else
       return end (r);
@@ -2298,7 +2308,9 @@ namespace TrilinosWrappers
     // past this line, or at the end of the
     // matrix
     for (size_type i=r+1; i<m(); ++i)
-      if (row_length(i) > 0)
+      if (in_local_range (i)
+         &&
+         (row_length(i) > 0))
         return iterator(this, i, 0);
 
     // if there is no such line, then take the
index abf80a79f4d92f70ffb36ecbd243d99e92f83003..4d808c7a3e61caf2c704bd38a89772eebc09feb5 100644 (file)
@@ -114,11 +114,15 @@ namespace TrilinosWrappers
     void
     AccessorBase::visit_present_row ()
     {
-      // if we are asked to visit the
-      // past-the-end line, then simply
-      // release all our caches and go on
-      // with life
-      if (this->a_row == matrix->m())
+      // if we are asked to visit the past-the-end line, then simply
+      // release all our caches and go on with life.
+      //
+      // do the same if the row we're supposed to visit is not locally
+      // owned. this is simply going to make non-locally owned rows
+      // look like they're empty
+      if ((this->a_row == matrix->m())
+         ||
+         (matrix->in_local_range (this->a_row) == false))
         {
           colnum_cache.reset ();
           value_cache.reset ();
@@ -126,8 +130,7 @@ namespace TrilinosWrappers
           return;
         }
 
-      // get a representation of the present
-      // row
+      // get a representation of the present row
       int ncols;
       TrilinosWrappers::types::int_type colnums = matrix->n();
       if (value_cache.get() == 0)
diff --git a/tests/mpi/trilinos_sparse_matrix_03.cc b/tests/mpi/trilinos_sparse_matrix_03.cc
new file mode 100644 (file)
index 0000000..b5ae003
--- /dev/null
@@ -0,0 +1,230 @@
+// ---------------------------------------------------------------------
+//
+// 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 TrilinosWrappers::SparseMatrix::iterator semantics
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test ()
+{
+  const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  const unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  const unsigned int n_rows = 3;
+  const unsigned int n_cols = 4;
+
+  IndexSet row_partitioning (n_rows);
+  IndexSet col_partitioning (n_cols);
+
+  if (n_procs == 1)
+    {
+      row_partitioning.add_range(0, n_rows);
+      col_partitioning.add_range(0, n_cols);
+    }
+  else if (n_procs == 2)
+    {
+      // row_partitioning should be { [0, 2), [2, n_rows) }
+      // col_partitioning should be { [0, 2), [2, n_cols) }
+      // col_relevant_set should be { [0, 3), [1, n_cols) }
+      if (my_id == 0)
+        {
+          row_partitioning.add_range(0, 2);
+          col_partitioning.add_range(0, 2);
+        }
+      else if (my_id == 1)
+        {
+          row_partitioning.add_range(2, n_rows);
+          col_partitioning.add_range(2, n_cols);
+        }
+    }
+  else
+    Assert (false, ExcNotImplemented());
+
+  TrilinosWrappers::SparsityPattern sp (row_partitioning,
+                                        col_partitioning, MPI_COMM_WORLD);
+  if (my_id == 0)
+    {
+      sp.add (0, 0);
+      sp.add (0, 2);
+    }
+  if ((n_procs == 1) || (my_id == 1))
+    sp.add(2,3);
+  sp.compress();
+
+  TrilinosWrappers::SparseMatrix A;
+  A.reinit (sp);
+  if (my_id==0)
+    {
+      A.set(0, 0, 0.1);
+      A.set(0, 2, 0.2);
+    }
+  if ((n_procs == 1) || (my_id == 1))
+    {
+      A.set(0, 0, 0.1);
+      A.set(0, 2, 0.2);
+      A.set(2,3, 0.3);
+    }
+
+  A.compress(VectorOperation::insert);
+
+  // now access elements by iterator. we know what should be in the
+  // matrix, just make sure
+  if (my_id == 0)
+    {
+      for (TrilinosWrappers::SparseMatrix::iterator p=A.begin(0);
+          p != A.end(0);
+          ++p)
+       {
+         AssertThrow (p->row() == 0, ExcInternalError());
+         if (p->column() == 0)
+           AssertThrow (p->value() == 0.1, ExcInternalError())
+         else if (p->column() == 2)
+           AssertThrow (p->value() == 0.2, ExcInternalError())
+         else
+           // well, we didn't write here, so the only thing that
+           // should be in there is a zero
+           AssertThrow (p->value() == 0.0, ExcInternalError());
+       }
+    }
+  else
+    {
+      for (TrilinosWrappers::SparseMatrix::iterator p=A.begin(2);
+          p != A.end(2);
+          ++p)
+       {
+         AssertThrow (p->row() == 2, ExcInternalError());
+         if (p->column() == 3)
+           AssertThrow (p->value() == 0.3, ExcInternalError())
+         else
+           // well, we didn't write here, so the only thing that
+           // should be in there is a zero
+           AssertThrow (p->value() == 0.0, ExcInternalError());
+       }
+    }
+
+  // now let each of the processors write something into the memory of
+  // the other. the values we can locally access are of course the
+  // same ones as before because we haven't communicated the values
+  // yet:
+  if (my_id == 0)
+    {
+      A.set (2, 3, 42.);
+      for (TrilinosWrappers::SparseMatrix::iterator p=A.begin(0);
+          p != A.end(0);
+          ++p)
+       {
+         AssertThrow (p->row() == 0, ExcInternalError());
+         if (p->column() == 0)
+           AssertThrow (p->value() == 0.1, ExcInternalError())
+         else if (p->column() == 2)
+           AssertThrow (p->value() == 0.2, ExcInternalError())
+         else
+           // well, we didn't write here, so the only thing that
+           // should be in there is a zero
+           AssertThrow (p->value() == 0.0, ExcInternalError());
+       }
+    }
+  else
+    {
+      A.set (0, 0, 108.);
+      for (TrilinosWrappers::SparseMatrix::iterator p=A.begin(2);
+          p != A.end(2);
+          ++p)
+       {
+         AssertThrow (p->row() == 2, ExcInternalError());
+         if (p->column() == 3)
+           AssertThrow (p->value() == 0.3, ExcInternalError())
+         else
+           // well, we didn't write here, so the only thing that
+           // should be in there is a zero
+           AssertThrow (p->value() == 0.0, ExcInternalError());
+       }
+    }
+
+  // then call compress() and ensure that we get the correct values
+  A.compress (VectorOperation::insert);
+  if (my_id == 0)
+    {
+      for (TrilinosWrappers::SparseMatrix::iterator p=A.begin(0);
+          p != A.end(0);
+          ++p)
+       {
+         AssertThrow (p->row() == 0, ExcInternalError());
+         if (p->column() == 0)
+           AssertThrow (p->value() == 108, ExcInternalError())
+         else if (p->column() == 2)
+           AssertThrow (p->value() == 0.2, ExcInternalError())
+         else
+           // well, we didn't write here, so the only thing that
+           // should be in there is a zero
+           AssertThrow (p->value() == 0.0, ExcInternalError());
+       }
+    }
+  else
+    {
+      for (TrilinosWrappers::SparseMatrix::iterator p=A.begin(2);
+          p != A.end(2);
+          ++p)
+       {
+         AssertThrow (p->row() == 2, ExcInternalError());
+         if (p->column() == 3)
+           AssertThrow (p->value() == 42, ExcInternalError())
+         else
+           // well, we didn't write here, so the only thing that
+           // should be in there is a zero
+           AssertThrow (p->value() == 0.0, ExcInternalError());
+       }
+    }  
+  
+  if (my_id == 0)
+    deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
+
+  const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  deallog.push(Utilities::int_to_string(myid));
+
+  if (myid == 0)
+    {
+      std::ofstream logfile("output");
+      deallog.attach(logfile);
+      deallog << std::setprecision(4);
+      deallog.depth_console(0);
+      deallog.threshold_double(1.e-10);
+
+      test();
+    }
+  else
+    test();
+
+}
diff --git a/tests/mpi/trilinos_sparse_matrix_03.with_trilinos=true.mpirun=2.output b/tests/mpi/trilinos_sparse_matrix_03.with_trilinos=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..be8d055
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL:0::OK
diff --git a/tests/mpi/trilinos_sparse_matrix_04.cc b/tests/mpi/trilinos_sparse_matrix_04.cc
new file mode 100644 (file)
index 0000000..2424813
--- /dev/null
@@ -0,0 +1,141 @@
+// ---------------------------------------------------------------------
+//
+// 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 TrilinosWrappers::SparseMatrix::iterator semantics. make sure
+// that rows not stored locally look like they're empty
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/parallel_vector.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test ()
+{
+  const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  const unsigned int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  const unsigned int n_rows = 3;
+  const unsigned int n_cols = 4;
+
+  IndexSet row_partitioning (n_rows);
+  IndexSet col_partitioning (n_cols);
+
+  if (n_procs == 1)
+    {
+      row_partitioning.add_range(0, n_rows);
+      col_partitioning.add_range(0, n_cols);
+    }
+  else if (n_procs == 2)
+    {
+      // row_partitioning should be { [0, 2), [2, n_rows) }
+      // col_partitioning should be { [0, 2), [2, n_cols) }
+      // col_relevant_set should be { [0, 3), [1, n_cols) }
+      if (my_id == 0)
+        {
+          row_partitioning.add_range(0, 2);
+          col_partitioning.add_range(0, 2);
+        }
+      else if (my_id == 1)
+        {
+          row_partitioning.add_range(2, n_rows);
+          col_partitioning.add_range(2, n_cols);
+        }
+    }
+  else
+    Assert (false, ExcNotImplemented());
+
+  TrilinosWrappers::SparsityPattern sp (row_partitioning,
+                                        col_partitioning, MPI_COMM_WORLD);
+  if (my_id == 0)
+    {
+      sp.add (0, 0);
+      sp.add (0, 2);
+    }
+  if ((n_procs == 1) || (my_id == 1))
+    sp.add(2,3);
+  sp.compress();
+
+  TrilinosWrappers::SparseMatrix A;
+  A.reinit (sp);
+  if (my_id==0)
+    {
+      A.set(0, 0, 0.1);
+      A.set(0, 2, 0.2);
+    }
+  if ((n_procs == 1) || (my_id == 1))
+    {
+      A.set(0, 0, 0.1);
+      A.set(0, 2, 0.2);
+      A.set(2,3, 0.3);
+    }
+
+  A.compress(VectorOperation::insert);
+
+  // now access elements by iterator. ensure that we can iterate over
+  // all rows but that iterators into rows not stored locally just
+  // look empty
+  for (TrilinosWrappers::SparseMatrix::iterator p=A.begin();
+       p != A.end();
+       ++p)
+    if (my_id == 0)
+      {
+       deallog << "Looking at entry (" << p->row() << ','
+               << p->column() << ") with value "
+               << p->value()
+               << std::endl;
+       AssertThrow (p->row() == 0, ExcInternalError());
+      }
+    else
+      {
+       AssertThrow (p->row() == 2, ExcInternalError());
+      }  
+  
+  if (my_id == 0)
+    deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, numbers::invalid_unsigned_int);
+
+  const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  unsigned int myid = Utilities::MPI::this_mpi_process (MPI_COMM_WORLD);
+  deallog.push(Utilities::int_to_string(myid));
+
+  if (myid == 0)
+    {
+      std::ofstream logfile("output");
+      deallog.attach(logfile);
+      deallog << std::setprecision(4);
+      deallog.depth_console(0);
+      deallog.threshold_double(1.e-10);
+
+      test();
+    }
+  else
+    test();
+
+}
diff --git a/tests/mpi/trilinos_sparse_matrix_04.with_trilinos=true.mpirun=2.output b/tests/mpi/trilinos_sparse_matrix_04.with_trilinos=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..bd2018a
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL:0::Looking at entry (0,0) with value 0.1000
+DEAL:0::Looking at entry (0,2) with value 0.2000
+DEAL:0::OK

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.