From 004d08bf15470f84ae59aa75a665dc89ef2f04be Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Wed, 5 Mar 2014 11:48:30 +0000 Subject: [PATCH] Temporary fix for Issue 191, allow sadd for Trilinos vectors git-svn-id: https://svn.dealii.org/trunk@32612 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/lac/trilinos_vector_base.h | 120 ++++++++++++------ deal.II/source/lac/trilinos_vector_base.cc | 17 ++- ...add_01.with_trilinos=true.mpirun=2.output} | 0 tests/mpi/trilinos_sadd_02.cc | 111 ++++++++++++++++ ...sadd_02.with_trilinos=true.mpirun=1.output | 45 +++++++ ...sadd_02.with_trilinos=true.mpirun=2.output | 85 +++++++++++++ 6 files changed, 336 insertions(+), 42 deletions(-) rename tests/mpi/{trilinos_sadd_01.mpirun=2.output => trilinos_sadd_01.with_trilinos=true.mpirun=2.output} (100%) create mode 100644 tests/mpi/trilinos_sadd_02.cc create mode 100644 tests/mpi/trilinos_sadd_02.with_trilinos=true.mpirun=1.output create mode 100644 tests/mpi/trilinos_sadd_02.with_trilinos=true.mpirun=2.output diff --git a/deal.II/include/deal.II/lac/trilinos_vector_base.h b/deal.II/include/deal.II/lac/trilinos_vector_base.h index dfb77a229e..84845480eb 100644 --- a/deal.II/include/deal.II/lac/trilinos_vector_base.h +++ b/deal.II/include/deal.II/lac/trilinos_vector_base.h @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2008 - 2014 by the deal.II authors +// Copyright (C) 2008 - 2013 by the deal.II authors // // This file is part of the deal.II library. // @@ -1383,7 +1383,7 @@ namespace TrilinosWrappers // use pre-allocated vector for non-local entries if it exists for // addition operation const TrilinosWrappers::types::int_type my_row = nonlocal_vector->Map().LID(static_cast(row)); - Assert(my_row != -1, + Assert(my_row != -1, ExcMessage("Attempted to write into off-processor vector entry " "that has not be specified as being writable upon " "initialization")); @@ -1704,14 +1704,22 @@ namespace TrilinosWrappers // if we have ghost values, do not allow // writing to this vector at all. Assert (!has_ghost_elements(), ExcGhostsPresent()); - Assert (local_size() == v.local_size(), - ExcDimensionMismatch(local_size(), v.local_size())); + Assert (size() == v.size(), + ExcDimensionMismatch (size(), v.size())); Assert (numbers::is_finite(s), ExcNumberNotFinite()); - const int ierr = vector->Update(1., *(v.vector), s); - - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + if(local_size() == v.local_size()) + { + const int ierr = vector->Update(1., *(v.vector), s); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + else + { + VectorBase tmp = v; + tmp *= s; + this->add(tmp, true); + } } @@ -1725,15 +1733,23 @@ namespace TrilinosWrappers // if we have ghost values, do not allow // writing to this vector at all. Assert (!has_ghost_elements(), ExcGhostsPresent()); - Assert (local_size() == v.local_size(), - ExcDimensionMismatch(local_size(), v.local_size())); - + Assert (size() == v.size(), + ExcDimensionMismatch (size(), v.size())); Assert (numbers::is_finite(s), ExcNumberNotFinite()); Assert (numbers::is_finite(a), ExcNumberNotFinite()); - const int ierr = vector->Update(a, *(v.vector), s); - - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + if(local_size() == v.local_size()) + { + const int ierr = vector->Update(a, *(v.vector), s); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + else + { + (*this)*=s; + VectorBase tmp = v; + tmp *= a; + this->add(tmp, true); + } } @@ -1749,18 +1765,29 @@ namespace TrilinosWrappers // if we have ghost values, do not allow // writing to this vector at all. Assert (!has_ghost_elements(), ExcGhostsPresent()); - Assert (local_size() == v.local_size(), - ExcDimensionMismatch(local_size(), v.local_size())); - Assert (local_size() == w.local_size(), - ExcDimensionMismatch(local_size(), w.local_size())); - + Assert (size() == v.size(), + ExcDimensionMismatch (size(), v.size())); + Assert (size() == w.size(), + ExcDimensionMismatch (size(), w.size())); Assert (numbers::is_finite(s), ExcNumberNotFinite()); Assert (numbers::is_finite(a), ExcNumberNotFinite()); Assert (numbers::is_finite(b), ExcNumberNotFinite()); - - const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s); - - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + if(local_size() == v.local_size() && local_size() == w.local_size()) + { + const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + else + { + (*this)*=s; + VectorBase tmp = v; + tmp *= a; + this->add(tmp, true); + tmp = w; + tmp *= b; + this->add(tmp, true); + } } @@ -1778,27 +1805,44 @@ namespace TrilinosWrappers // if we have ghost values, do not allow // writing to this vector at all. Assert (!has_ghost_elements(), ExcGhostsPresent()); - Assert (local_size() == v.local_size(), - ExcDimensionMismatch(local_size(), v.local_size())); - Assert (local_size() == w.local_size(), - ExcDimensionMismatch(local_size(), w.local_size())); - Assert (local_size() == x.local_size(), - ExcDimensionMismatch(local_size(), x.local_size())); - + Assert (size() == v.size(), + ExcDimensionMismatch (size(), v.size())); + Assert (size() == w.size(), + ExcDimensionMismatch (size(), w.size())); + Assert (size() == x.size(), + ExcDimensionMismatch (size(), x.size())); Assert (numbers::is_finite(s), ExcNumberNotFinite()); Assert (numbers::is_finite(a), ExcNumberNotFinite()); Assert (numbers::is_finite(b), ExcNumberNotFinite()); Assert (numbers::is_finite(c), ExcNumberNotFinite()); - // Update member can only - // input two other vectors so - // do it in two steps - const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - - const int jerr = vector->Update(c, *(x.vector), 1.); - Assert (jerr == 0, ExcTrilinosError(jerr)); - (void)jerr; // removes -Wunused-parameter warning in optimized mode + if(local_size() == v.local_size() + && local_size() == w.local_size() + && local_size() == x.local_size()) + { + // Update member can only + // input two other vectors so + // do it in two steps + const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + const int jerr = vector->Update(c, *(x.vector), 1.); + Assert (jerr == 0, ExcTrilinosError(jerr)); + (void)jerr; // removes -Wunused-parameter warning in optimized mode + } + else + { + (*this)*=s; + VectorBase tmp = v; + tmp *= a; + this->add(tmp, true); + tmp = w; + tmp *= b; + this->add(tmp, true); + tmp = x; + tmp *= c; + this->add(tmp, true); + } } diff --git a/deal.II/source/lac/trilinos_vector_base.cc b/deal.II/source/lac/trilinos_vector_base.cc index 95a9b1387b..da1e491f16 100644 --- a/deal.II/source/lac/trilinos_vector_base.cc +++ b/deal.II/source/lac/trilinos_vector_base.cc @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2008 - 2013 by the deal.II authors +// Copyright (C) 2008 - 2014 by the deal.II authors // // This file is part of the deal.II library. // @@ -291,15 +291,24 @@ namespace TrilinosWrappers *this += v; else { + Assert (!has_ghost_elements(), ExcGhostsPresent()); AssertThrow (size() == v.size(), ExcDimensionMismatch (size(), v.size())); - Epetra_Import data_exchange (vector->Map(), v.vector->Map()); + //HACK: For some reason the following doesn't work properly, so do it manually + // Epetra_Import data_exchange (vector->Map(), v.vector->Map()); + // int ierr = vector->Import(*v.vector, data_exchange, Add); + // AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + // last_action = Add; - int ierr = vector->Import(*v.vector, data_exchange, Add); + Epetra_MultiVector dummy(vector->Map(), 1, false); + Epetra_Import data_exchange (dummy.Map(), v.vector->Map()); + + int ierr = dummy.Import(*v.vector, data_exchange, Insert); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - last_action = Insert; + ierr = vector->Update (1.0, dummy, 1.0); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } } diff --git a/tests/mpi/trilinos_sadd_01.mpirun=2.output b/tests/mpi/trilinos_sadd_01.with_trilinos=true.mpirun=2.output similarity index 100% rename from tests/mpi/trilinos_sadd_01.mpirun=2.output rename to tests/mpi/trilinos_sadd_01.with_trilinos=true.mpirun=2.output diff --git a/tests/mpi/trilinos_sadd_02.cc b/tests/mpi/trilinos_sadd_02.cc new file mode 100644 index 0000000000..a5723bedee --- /dev/null +++ b/tests/mpi/trilinos_sadd_02.cc @@ -0,0 +1,111 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2011 - 2013 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. +// +// --------------------------------------------------------------------- + + + +// check correct behaviour of sadd of Trilinos vectors +// if they have different Epetra maps + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include + +void test () +{ + TrilinosWrappers::MPI::Vector ghosted, distributed; + //All processes should own 10 entries + const int entries_per_process = 10; + const int n_proc = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + IndexSet locally_owned(entries_per_process*n_proc); + IndexSet locally_relevant(entries_per_process*n_proc); + const int my_id = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + const int begin_index = my_id*entries_per_process; + const int end_index = (my_id+1)*entries_per_process; + const int local_begin = std::max(0, begin_index-entries_per_process/2); + const int local_end = entries_per_process*n_proc; + + locally_owned.add_range(begin_index, end_index); + locally_relevant.add_range + (local_begin, local_end); + + distributed.reinit(locally_owned, MPI_COMM_WORLD); + ghosted.reinit (locally_owned, locally_relevant, MPI_COMM_WORLD); + + distributed=1.; + ghosted=distributed; + distributed.sadd (2., ghosted); + ghosted = distributed; + deallog << "sadd (s, v)" <