<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.
* 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;
*
* @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
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();
SparseMatrix::const_iterator
SparseMatrix::begin() const
{
- return const_iterator(this, 0, 0);
+ return begin(0);
}
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);
// 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
SparseMatrix::iterator
SparseMatrix::begin()
{
- return iterator(this, 0, 0);
+ return begin(0);
}
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);
// 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
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 ();
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)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+
+}
--- /dev/null
+
+DEAL:0::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff