From: Wolfgang Bangerth Date: Mon, 9 Mar 2015 05:11:44 +0000 (-0500) Subject: Ensure that we can run iterators over Trilinos matrices that are not stored locally. X-Git-Tag: v8.3.0-rc1~380^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f275de690eaa384672962874eab969a84bf3dcfa;p=dealii.git Ensure that we can run iterators over Trilinos matrices that are not stored locally. --- diff --git a/doc/news/changes.h b/doc/news/changes.h index e9db05ecdf..67cbd6fd83 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -364,6 +364,15 @@ inconvenience this causes.

Specific improvements

    +
  1. 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. +
    + (Wolfgang Bangerth, 2015/03/08) +
    2. +
  2. 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. diff --git a/include/deal.II/lac/trilinos_sparse_matrix.h b/include/deal.II/lac/trilinos_sparse_matrix.h index 283a1ec372..10cabcf12e 100644 --- a/include/deal.II/lac/trilinos_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_sparse_matrix.h @@ -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 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 0) + if (in_local_range (i) + && + (row_length(i) > 0)) return iterator(this, i, 0); // if there is no such line, then take the diff --git a/source/lac/trilinos_sparse_matrix.cc b/source/lac/trilinos_sparse_matrix.cc index abf80a79f4..4d808c7a3e 100644 --- a/source/lac/trilinos_sparse_matrix.cc +++ b/source/lac/trilinos_sparse_matrix.cc @@ -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 index 0000000000..b5ae003844 --- /dev/null +++ b/tests/mpi/trilinos_sparse_matrix_03.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include + + +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 index 0000000000..be8d055f86 --- /dev/null +++ b/tests/mpi/trilinos_sparse_matrix_03.with_trilinos=true.mpirun=2.output @@ -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 index 0000000000..2424813027 --- /dev/null +++ b/tests/mpi/trilinos_sparse_matrix_04.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include + + +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 index 0000000000..bd2018a681 --- /dev/null +++ b/tests/mpi/trilinos_sparse_matrix_04.with_trilinos=true.mpirun=2.output @@ -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