From 40a874a6b150c863b5be5577b19da5d06af66e4c Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 22 Feb 2019 15:04:18 +0100 Subject: [PATCH] Move implementation to templates file --- .../lac/trilinos_tpetra_vector.templates.h | 681 ++++++++++++++++++ source/lac/trilinos_tpetra_vector.cc | 646 +---------------- 2 files changed, 682 insertions(+), 645 deletions(-) create mode 100644 include/deal.II/lac/trilinos_tpetra_vector.templates.h diff --git a/include/deal.II/lac/trilinos_tpetra_vector.templates.h b/include/deal.II/lac/trilinos_tpetra_vector.templates.h new file mode 100644 index 0000000000..b102775201 --- /dev/null +++ b/include/deal.II/lac/trilinos_tpetra_vector.templates.h @@ -0,0 +1,681 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 - 2019 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii__trilinos_tpetra_vector_templates_h +#define dealii__trilinos_tpetra_vector_templates_h + +#include + +#include + +#ifdef DEAL_II_TRILINOS_WITH_TPETRA + +# ifdef DEAL_II_WITH_MPI + +# include + +# include + +# include + +# include +# include +# include + +# include + + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + namespace TpetraWrappers + { + template + Vector::Vector() + : Subscriptor() + , vector(new Tpetra::Vector( + Teuchos::RCP>( + new Tpetra::Map( + 0, + 0, + Utilities::Trilinos::tpetra_comm_self())))) + {} + + + + template + Vector::Vector(const Vector &V) + : Subscriptor() + , vector(new Tpetra::Vector( + V.trilinos_vector(), + Teuchos::Copy)) + {} + + + + template + Vector::Vector(const IndexSet ¶llel_partitioner, + const MPI_Comm &communicator) + : Subscriptor() + , vector(new Tpetra::Vector( + Teuchos::rcp(new Tpetra::Map( + parallel_partitioner.make_tpetra_map(communicator, false))))) + {} + + + + template + void + Vector::reinit(const IndexSet ¶llel_partitioner, + const MPI_Comm &communicator, + const bool omit_zeroing_entries) + { + Tpetra::Map input_map = + parallel_partitioner.make_tpetra_map(communicator, false); + if (vector->getMap()->isSameAs(input_map) == false) + vector = std_cxx14::make_unique< + Tpetra::Vector>(Teuchos::rcp( + new Tpetra::Map(input_map))); + else if (omit_zeroing_entries == false) + { + vector->putScalar(0.); + } + } + + + + template + void + Vector::reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries) + { + // Check that casting will work. + Assert(dynamic_cast *>(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast &>(V); + + reinit(down_V.locally_owned_elements(), + down_V.get_mpi_communicator(), + omit_zeroing_entries); + } + + + + template + Vector & + Vector::operator=(const Vector &V) + { + // Distinguish three cases: + // - First case: both vectors have the same layout. + // - Second case: both vectors have the same size but different layout. + // - Third case: the vectors have different size. + if (vector->getMap()->isSameAs(*(V.trilinos_vector().getMap()))) + *vector = V.trilinos_vector(); + else + { + if (size() == V.size()) + { + Tpetra::Import data_exchange( + vector->getMap(), V.trilinos_vector().getMap()); + + vector->doImport(V.trilinos_vector(), + data_exchange, + Tpetra::REPLACE); + } + else + vector = std_cxx14::make_unique< + Tpetra::Vector>( + V.trilinos_vector()); + } + + return *this; + } + + + + template + Vector & + Vector::operator=(const Number s) + { + Assert(s == Number(0.), + ExcMessage("Only 0 can be assigned to a vector.")); + + vector->putScalar(s); + + return *this; + } + + + + template + void + Vector::import( + const ReadWriteVector & V, + VectorOperation::values operation, + std::shared_ptr communication_pattern) + { + // If no communication pattern is given, create one. Otherwsie, use the + // one given. + if (communication_pattern == nullptr) + { + // The first time import is called, a communication pattern is + // created. Check if the communication pattern already exists and if + // it can be reused. + if ((source_stored_elements.size() != + V.get_stored_elements().size()) || + (source_stored_elements != V.get_stored_elements())) + { + const Teuchos::MpiComm *mpi_comm = + dynamic_cast *>( + vector->getMap()->getComm().get()); + Assert(mpi_comm != nullptr, ExcInternalError()); + create_tpetra_comm_pattern(V.get_stored_elements(), + *(mpi_comm->getRawMpiComm())()); + } + } + else + { + tpetra_comm_pattern = std::dynamic_pointer_cast< + const TpetraWrappers::CommunicationPattern>(communication_pattern); + AssertThrow( + tpetra_comm_pattern != nullptr, + ExcMessage( + std::string("The communication pattern is not of type ") + + "LinearAlgebra::TpetraWrappers::CommunicationPattern.")); + } + + Tpetra::Export tpetra_export( + tpetra_comm_pattern->get_tpetra_export()); + Tpetra::Vector source_vector( + tpetra_export.getSourceMap()); + + source_vector.template sync(); + auto x_2d = source_vector.template getLocalView(); + auto x_1d = Kokkos::subview(x_2d, Kokkos::ALL(), 0); + source_vector.template modify(); + const size_t localLength = source_vector.getLocalLength(); + auto values_it = V.begin(); + for (size_t k = 0; k < localLength; ++k) + x_1d(k) = *values_it++; + source_vector.template sync< + typename Tpetra::Vector:: + device_type::memory_space>(); + if (operation == VectorOperation::insert) + vector->doExport(source_vector, tpetra_export, Tpetra::REPLACE); + else if (operation == VectorOperation::add) + vector->doExport(source_vector, tpetra_export, Tpetra::ADD); + else + AssertThrow(false, ExcNotImplemented()); + } + + + + template + Vector & + Vector::operator*=(const Number factor) + { + AssertIsFinite(factor); + vector->scale(factor); + + return *this; + } + + + + template + Vector & + Vector::operator/=(const Number factor) + { + AssertIsFinite(factor); + Assert(factor != Number(0.), ExcZero()); + *this *= Number(1.) / factor; + + return *this; + } + + + + template + Vector & + Vector::operator+=(const VectorSpaceVector &V) + { + // Check that casting will work. + Assert(dynamic_cast *>(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast &>(V); + // If the maps are the same we can update right away. + if (vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap()))) + { + vector->update(1., down_V.trilinos_vector(), 1.); + } + else + { + Assert(this->size() == down_V.size(), + ExcDimensionMismatch(this->size(), down_V.size())); + + // TODO: Tpetra doesn't have a combine mode that also updates local + // elements, maybe there is a better workaround. + Tpetra::Vector dummy( + vector->getMap(), false); + Tpetra::Import data_exchange( + down_V.trilinos_vector().getMap(), dummy.getMap()); + + dummy.doImport(down_V.trilinos_vector(), + data_exchange, + Tpetra::INSERT); + + vector->update(1.0, dummy, 1.0); + } + + return *this; + } + + + + template + Vector & + Vector::operator-=(const VectorSpaceVector &V) + { + this->add(-1., V); + + return *this; + } + + + + template + Number Vector::operator*(const VectorSpaceVector &V) const + { + // Check that casting will work. + Assert(dynamic_cast *>(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast &>(V); + Assert(this->size() == down_V.size(), + ExcDimensionMismatch(this->size(), down_V.size())); + Assert(vector->getMap()->isSameAs(*down_V.trilinos_vector().getMap()), + ExcDifferentParallelPartitioning()); + + return vector->dot(down_V.trilinos_vector()); + } + + + + template + void + Vector::add(const Number a) + { + AssertIsFinite(a); + + vector->template sync(); + auto vector_2d = vector->template getLocalView(); + auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); + vector->template modify(); + const size_t localLength = vector->getLocalLength(); + for (size_t k = 0; k < localLength; ++k) + { + vector_1d(k) += a; + } + vector->template sync< + typename Tpetra::Vector:: + device_type::memory_space>(); + } + + + + template + void + Vector::add(const Number a, const VectorSpaceVector &V) + { + // Check that casting will work. + Assert(dynamic_cast *>(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast &>(V); + AssertIsFinite(a); + Assert(vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())), + ExcDifferentParallelPartitioning()); + + vector->update(a, down_V.trilinos_vector(), 1.); + } + + + + template + void + Vector::add(const Number a, + const VectorSpaceVector &V, + const Number b, + const VectorSpaceVector &W) + { + // Check that casting will work. + Assert(dynamic_cast *>(&V) != nullptr, + ExcVectorTypeNotCompatible()); + // Check that casting will work. + Assert(dynamic_cast *>(&W) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast &>(V); + // Downcast W. If fails, throws an exception. + const Vector &down_W = dynamic_cast &>(W); + Assert(vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())), + ExcDifferentParallelPartitioning()); + Assert(vector->getMap()->isSameAs(*(down_W.trilinos_vector().getMap())), + ExcDifferentParallelPartitioning()); + AssertIsFinite(a); + AssertIsFinite(b); + + vector->update( + a, down_V.trilinos_vector(), b, down_W.trilinos_vector(), 1.); + } + + + + template + void + Vector::sadd(const Number s, + const Number a, + const VectorSpaceVector &V) + { + // Check that casting will work. + Assert(dynamic_cast *>(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + *this *= s; + // Downcast V. It fails, throws an exception. + const Vector &down_V = dynamic_cast &>(V); + Vector tmp(down_V); + tmp *= a; + *this += tmp; + } + + + + template + void + Vector::scale(const VectorSpaceVector &scaling_factors) + { + // Check that casting will work. + Assert(dynamic_cast *>(&scaling_factors) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast scaling_factors. If fails, throws an exception. + const Vector &down_scaling_factors = + dynamic_cast &>(scaling_factors); + Assert(vector->getMap()->isSameAs( + *(down_scaling_factors.trilinos_vector().getMap())), + ExcDifferentParallelPartitioning()); + + vector->elementWiseMultiply(1., + *down_scaling_factors.vector, + *vector, + 0.); + } + + + + template + void + Vector::equ(const Number a, const VectorSpaceVector &V) + { + // Check that casting will work. + Assert(dynamic_cast *>(&V) != nullptr, + ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast &>(V); + // If we don't have the same map, copy. + if (vector->getMap()->isSameAs(*down_V.trilinos_vector().getMap()) == + false) + this->sadd(0., a, V); + else + { + // Otherwise, just update + vector->update(a, down_V.trilinos_vector(), 0.); + } + } + + + + template + bool + Vector::all_zero() const + { + // get a representation of the vector and + // loop over all the elements + Number * start_ptr = vector->getDataNonConst().get(); + const Number *ptr = start_ptr, + *eptr = start_ptr + vector->getLocalLength(); + unsigned int flag = 0; + while (ptr != eptr) + { + if (*ptr != Number(0)) + { + flag = 1; + break; + } + ++ptr; + } + + // Check that the vector is zero on _all_ processors. + const Teuchos::MpiComm *mpi_comm = + dynamic_cast *>( + vector->getMap()->getComm().get()); + Assert(mpi_comm != nullptr, ExcInternalError()); + unsigned int num_nonzero = + Utilities::MPI::sum(flag, *(mpi_comm->getRawMpiComm())()); + + return num_nonzero == 0; + } + + + + template + Number + Vector::mean_value() const + { + return vector->meanValue(); + } + + + + template + typename LinearAlgebra::VectorSpaceVector::real_type + Vector::l1_norm() const + { + return vector->norm1(); + } + + + + template + typename LinearAlgebra::VectorSpaceVector::real_type + Vector::l2_norm() const + { + return vector->norm2(); + } + + + + template + typename LinearAlgebra::VectorSpaceVector::real_type + Vector::linfty_norm() const + { + return vector->normInf(); + } + + + + template + Number + Vector::add_and_dot(const Number a, + const VectorSpaceVector &V, + const VectorSpaceVector &W) + { + this->add(a, V); + + return *this * W; + } + + + + template + typename Vector::size_type + Vector::size() const + { + return vector->getGlobalLength(); + } + + + + template + MPI_Comm + Vector::get_mpi_communicator() const + { + const auto tpetra_comm = dynamic_cast *>( + vector->getMap()->getComm().get()); + Assert(tpetra_comm != nullptr, ExcInternalError()); + return *(tpetra_comm->getRawMpiComm())(); + } + + + + template + ::dealii::IndexSet + Vector::locally_owned_elements() const + { + IndexSet is(size()); + + // easy case: local range is contiguous + if (vector->getMap()->isContiguous()) + { + is.add_range(vector->getMap()->getMinGlobalIndex(), + vector->getMap()->getMaxGlobalIndex() + 1); + } + else if (vector->getLocalLength() > 0) + { + const size_type n_indices = vector->getLocalLength(); + std::vector vector_indices; + vector_indices.reserve(n_indices); + for (unsigned int i = 0; i < n_indices; ++i) + vector_indices.push_back(vector->getMap()->getGlobalElement(i)); + + is.add_indices(vector_indices.data(), + vector_indices.data() + n_indices); + } + is.compress(); + + return is; + } + + + + template + const Tpetra::Vector & + Vector::trilinos_vector() const + { + return *vector; + } + + + + template + Tpetra::Vector & + Vector::trilinos_vector() + { + return *vector; + } + + + + template + void + Vector::print(std::ostream & out, + const unsigned int precision, + const bool scientific, + const bool across) const + { + AssertThrow(out, ExcIO()); + boost::io::ios_flags_saver restore_flags(out); + + // Get a representation of the vector and loop over all + // the elements + const auto val = vector->get1dView(); + + out.precision(precision); + if (scientific) + out.setf(std::ios::scientific, std::ios::floatfield); + else + out.setf(std::ios::fixed, std::ios::floatfield); + + vector->template sync(); + auto vector_2d = vector->template getLocalView(); + auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); + const size_t local_length = vector->getLocalLength(); + + if (across) + for (unsigned int i = 0; i < local_length; ++i) + out << vector_1d(i) << ' '; + else + for (unsigned int i = 0; i < local_length; ++i) + out << vector_1d(i) << std::endl; + out << std::endl; + + // restore the representation + // of the vector + AssertThrow(out, ExcIO()); + } + + + + template + std::size_t + Vector::memory_consumption() const + { + return sizeof(*this) + + vector->getLocalLength() * + (sizeof(Number) + sizeof(TrilinosWrappers::types::int_type)); + } + + + + template + void + Vector::create_tpetra_comm_pattern(const IndexSet &source_index_set, + const MPI_Comm &mpi_comm) + { + source_stored_elements = source_index_set; + tpetra_comm_pattern = + std::make_shared( + locally_owned_elements(), source_index_set, mpi_comm); + } + } // namespace TpetraWrappers +} // namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +# endif + +#endif + +#endif diff --git a/source/lac/trilinos_tpetra_vector.cc b/source/lac/trilinos_tpetra_vector.cc index ea218f8309..b269bdabcf 100644 --- a/source/lac/trilinos_tpetra_vector.cc +++ b/source/lac/trilinos_tpetra_vector.cc @@ -13,660 +13,17 @@ // // --------------------------------------------------------------------- -#include - -#include +#include #ifdef DEAL_II_TRILINOS_WITH_TPETRA - # ifdef DEAL_II_WITH_MPI -# include - -# include - -# include - -# include -# include -# include - -# include - - DEAL_II_NAMESPACE_OPEN namespace LinearAlgebra { namespace TpetraWrappers { - template - Vector::Vector() - : Subscriptor() - , vector(new Tpetra::Vector( - Teuchos::RCP>( - new Tpetra::Map( - 0, - 0, - Utilities::Trilinos::tpetra_comm_self())))) - {} - - - - template - Vector::Vector(const Vector &V) - : Subscriptor() - , vector(new Tpetra::Vector( - V.trilinos_vector(), - Teuchos::Copy)) - {} - - - - template - Vector::Vector(const IndexSet ¶llel_partitioner, - const MPI_Comm &communicator) - : Subscriptor() - , vector(new Tpetra::Vector( - Teuchos::rcp(new Tpetra::Map( - parallel_partitioner.make_tpetra_map(communicator, false))))) - {} - - - - template - void - Vector::reinit(const IndexSet ¶llel_partitioner, - const MPI_Comm &communicator, - const bool omit_zeroing_entries) - { - Tpetra::Map input_map = - parallel_partitioner.make_tpetra_map(communicator, false); - if (vector->getMap()->isSameAs(input_map) == false) - vector = std_cxx14::make_unique< - Tpetra::Vector>(Teuchos::rcp( - new Tpetra::Map(input_map))); - else if (omit_zeroing_entries == false) - { - vector->putScalar(0.); - } - } - - - - template - void - Vector::reinit(const VectorSpaceVector &V, - const bool omit_zeroing_entries) - { - // Check that casting will work. - Assert(dynamic_cast *>(&V) != nullptr, - ExcVectorTypeNotCompatible()); - - // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast &>(V); - - reinit(down_V.locally_owned_elements(), - down_V.get_mpi_communicator(), - omit_zeroing_entries); - } - - - - template - Vector & - Vector::operator=(const Vector &V) - { - // Distinguish three cases: - // - First case: both vectors have the same layout. - // - Second case: both vectors have the same size but different layout. - // - Third case: the vectors have different size. - if (vector->getMap()->isSameAs(*(V.trilinos_vector().getMap()))) - *vector = V.trilinos_vector(); - else - { - if (size() == V.size()) - { - Tpetra::Import data_exchange( - vector->getMap(), V.trilinos_vector().getMap()); - - vector->doImport(V.trilinos_vector(), - data_exchange, - Tpetra::REPLACE); - } - else - vector = std_cxx14::make_unique< - Tpetra::Vector>( - V.trilinos_vector()); - } - - return *this; - } - - - - template - Vector & - Vector::operator=(const Number s) - { - Assert(s == Number(0.), - ExcMessage("Only 0 can be assigned to a vector.")); - - vector->putScalar(s); - - return *this; - } - - - - template - void - Vector::import( - const ReadWriteVector & V, - VectorOperation::values operation, - std::shared_ptr communication_pattern) - { - // If no communication pattern is given, create one. Otherwsie, use the - // one given. - if (communication_pattern == nullptr) - { - // The first time import is called, a communication pattern is - // created. Check if the communication pattern already exists and if - // it can be reused. - if ((source_stored_elements.size() != - V.get_stored_elements().size()) || - (source_stored_elements != V.get_stored_elements())) - { - const Teuchos::MpiComm *mpi_comm = - dynamic_cast *>( - vector->getMap()->getComm().get()); - Assert(mpi_comm != nullptr, ExcInternalError()); - create_tpetra_comm_pattern(V.get_stored_elements(), - *(mpi_comm->getRawMpiComm())()); - } - } - else - { - tpetra_comm_pattern = std::dynamic_pointer_cast< - const TpetraWrappers::CommunicationPattern>(communication_pattern); - AssertThrow( - tpetra_comm_pattern != nullptr, - ExcMessage( - std::string("The communication pattern is not of type ") + - "LinearAlgebra::TpetraWrappers::CommunicationPattern.")); - } - - Tpetra::Export tpetra_export( - tpetra_comm_pattern->get_tpetra_export()); - Tpetra::Vector source_vector( - tpetra_export.getSourceMap()); - - source_vector.template sync(); - auto x_2d = source_vector.template getLocalView(); - auto x_1d = Kokkos::subview(x_2d, Kokkos::ALL(), 0); - source_vector.template modify(); - const size_t localLength = source_vector.getLocalLength(); - auto values_it = V.begin(); - for (size_t k = 0; k < localLength; ++k) - x_1d(k) = *values_it++; - source_vector.template sync< - typename Tpetra::Vector:: - device_type::memory_space>(); - if (operation == VectorOperation::insert) - vector->doExport(source_vector, tpetra_export, Tpetra::REPLACE); - else if (operation == VectorOperation::add) - vector->doExport(source_vector, tpetra_export, Tpetra::ADD); - else - AssertThrow(false, ExcNotImplemented()); - } - - - - template - Vector & - Vector::operator*=(const Number factor) - { - AssertIsFinite(factor); - vector->scale(factor); - - return *this; - } - - - - template - Vector & - Vector::operator/=(const Number factor) - { - AssertIsFinite(factor); - Assert(factor != Number(0.), ExcZero()); - *this *= Number(1.) / factor; - - return *this; - } - - - - template - Vector & - Vector::operator+=(const VectorSpaceVector &V) - { - // Check that casting will work. - Assert(dynamic_cast *>(&V) != nullptr, - ExcVectorTypeNotCompatible()); - - // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast &>(V); - // If the maps are the same we can update right away. - if (vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap()))) - { - vector->update(1., down_V.trilinos_vector(), 1.); - } - else - { - Assert(this->size() == down_V.size(), - ExcDimensionMismatch(this->size(), down_V.size())); - - // TODO: Tpetra doesn't have a combine mode that also updates local - // elements, maybe there is a better workaround. - Tpetra::Vector dummy( - vector->getMap(), false); - Tpetra::Import data_exchange( - down_V.trilinos_vector().getMap(), dummy.getMap()); - - dummy.doImport(down_V.trilinos_vector(), - data_exchange, - Tpetra::INSERT); - - vector->update(1.0, dummy, 1.0); - } - - return *this; - } - - - - template - Vector & - Vector::operator-=(const VectorSpaceVector &V) - { - this->add(-1., V); - - return *this; - } - - - - template - Number Vector::operator*(const VectorSpaceVector &V) const - { - // Check that casting will work. - Assert(dynamic_cast *>(&V) != nullptr, - ExcVectorTypeNotCompatible()); - - // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast &>(V); - Assert(this->size() == down_V.size(), - ExcDimensionMismatch(this->size(), down_V.size())); - Assert(vector->getMap()->isSameAs(*down_V.trilinos_vector().getMap()), - ExcDifferentParallelPartitioning()); - - return vector->dot(down_V.trilinos_vector()); - } - - - - template - void - Vector::add(const Number a) - { - AssertIsFinite(a); - - vector->template sync(); - auto vector_2d = vector->template getLocalView(); - auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); - vector->template modify(); - const size_t localLength = vector->getLocalLength(); - for (size_t k = 0; k < localLength; ++k) - { - vector_1d(k) += a; - } - vector->template sync< - typename Tpetra::Vector:: - device_type::memory_space>(); - } - - - - template - void - Vector::add(const Number a, const VectorSpaceVector &V) - { - // Check that casting will work. - Assert(dynamic_cast *>(&V) != nullptr, - ExcVectorTypeNotCompatible()); - - // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast &>(V); - AssertIsFinite(a); - Assert(vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())), - ExcDifferentParallelPartitioning()); - - vector->update(a, down_V.trilinos_vector(), 1.); - } - - - - template - void - Vector::add(const Number a, - const VectorSpaceVector &V, - const Number b, - const VectorSpaceVector &W) - { - // Check that casting will work. - Assert(dynamic_cast *>(&V) != nullptr, - ExcVectorTypeNotCompatible()); - // Check that casting will work. - Assert(dynamic_cast *>(&W) != nullptr, - ExcVectorTypeNotCompatible()); - - // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast &>(V); - // Downcast W. If fails, throws an exception. - const Vector &down_W = dynamic_cast &>(W); - Assert(vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())), - ExcDifferentParallelPartitioning()); - Assert(vector->getMap()->isSameAs(*(down_W.trilinos_vector().getMap())), - ExcDifferentParallelPartitioning()); - AssertIsFinite(a); - AssertIsFinite(b); - - vector->update( - a, down_V.trilinos_vector(), b, down_W.trilinos_vector(), 1.); - } - - - - template - void - Vector::sadd(const Number s, - const Number a, - const VectorSpaceVector &V) - { - // Check that casting will work. - Assert(dynamic_cast *>(&V) != nullptr, - ExcVectorTypeNotCompatible()); - - *this *= s; - // Downcast V. It fails, throws an exception. - const Vector &down_V = dynamic_cast &>(V); - Vector tmp(down_V); - tmp *= a; - *this += tmp; - } - - - - template - void - Vector::scale(const VectorSpaceVector &scaling_factors) - { - // Check that casting will work. - Assert(dynamic_cast *>(&scaling_factors) != nullptr, - ExcVectorTypeNotCompatible()); - - // Downcast scaling_factors. If fails, throws an exception. - const Vector &down_scaling_factors = - dynamic_cast &>(scaling_factors); - Assert(vector->getMap()->isSameAs( - *(down_scaling_factors.trilinos_vector().getMap())), - ExcDifferentParallelPartitioning()); - - vector->elementWiseMultiply(1., - *down_scaling_factors.vector, - *vector, - 0.); - } - - - - template - void - Vector::equ(const Number a, const VectorSpaceVector &V) - { - // Check that casting will work. - Assert(dynamic_cast *>(&V) != nullptr, - ExcVectorTypeNotCompatible()); - - // Downcast V. If fails, throws an exception. - const Vector &down_V = dynamic_cast &>(V); - // If we don't have the same map, copy. - if (vector->getMap()->isSameAs(*down_V.trilinos_vector().getMap()) == - false) - this->sadd(0., a, V); - else - { - // Otherwise, just update - vector->update(a, down_V.trilinos_vector(), 0.); - } - } - - - - template - bool - Vector::all_zero() const - { - // get a representation of the vector and - // loop over all the elements - Number * start_ptr = vector->getDataNonConst().get(); - const Number *ptr = start_ptr, - *eptr = start_ptr + vector->getLocalLength(); - unsigned int flag = 0; - while (ptr != eptr) - { - if (*ptr != Number(0)) - { - flag = 1; - break; - } - ++ptr; - } - - // Check that the vector is zero on _all_ processors. - const Teuchos::MpiComm *mpi_comm = - dynamic_cast *>( - vector->getMap()->getComm().get()); - Assert(mpi_comm != nullptr, ExcInternalError()); - unsigned int num_nonzero = - Utilities::MPI::sum(flag, *(mpi_comm->getRawMpiComm())()); - - return num_nonzero == 0; - } - - - - template - Number - Vector::mean_value() const - { - return vector->meanValue(); - } - - - - template - typename LinearAlgebra::VectorSpaceVector::real_type - Vector::l1_norm() const - { - return vector->norm1(); - } - - - - template - typename LinearAlgebra::VectorSpaceVector::real_type - Vector::l2_norm() const - { - return vector->norm2(); - } - - - - template - typename LinearAlgebra::VectorSpaceVector::real_type - Vector::linfty_norm() const - { - return vector->normInf(); - } - - - - template - Number - Vector::add_and_dot(const Number a, - const VectorSpaceVector &V, - const VectorSpaceVector &W) - { - this->add(a, V); - - return *this * W; - } - - - - template - typename Vector::size_type - Vector::size() const - { - return vector->getGlobalLength(); - } - - - - template - MPI_Comm - Vector::get_mpi_communicator() const - { - const auto tpetra_comm = dynamic_cast *>( - vector->getMap()->getComm().get()); - Assert(tpetra_comm != nullptr, ExcInternalError()); - return *(tpetra_comm->getRawMpiComm())(); - } - - - - template - ::dealii::IndexSet - Vector::locally_owned_elements() const - { - IndexSet is(size()); - - // easy case: local range is contiguous - if (vector->getMap()->isContiguous()) - { - is.add_range(vector->getMap()->getMinGlobalIndex(), - vector->getMap()->getMaxGlobalIndex() + 1); - } - else if (vector->getLocalLength() > 0) - { - const size_type n_indices = vector->getLocalLength(); - std::vector vector_indices; - vector_indices.reserve(n_indices); - for (unsigned int i = 0; i < n_indices; ++i) - vector_indices.push_back(vector->getMap()->getGlobalElement(i)); - - is.add_indices(vector_indices.data(), - vector_indices.data() + n_indices); - } - is.compress(); - - return is; - } - - - - template - const Tpetra::Vector & - Vector::trilinos_vector() const - { - return *vector; - } - - - - template - Tpetra::Vector & - Vector::trilinos_vector() - { - return *vector; - } - - - - template - void - Vector::print(std::ostream & out, - const unsigned int precision, - const bool scientific, - const bool across) const - { - AssertThrow(out, ExcIO()); - boost::io::ios_flags_saver restore_flags(out); - - // Get a representation of the vector and loop over all - // the elements - const auto val = vector->get1dView(); - - out.precision(precision); - if (scientific) - out.setf(std::ios::scientific, std::ios::floatfield); - else - out.setf(std::ios::fixed, std::ios::floatfield); - - vector->template sync(); - auto vector_2d = vector->template getLocalView(); - auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); - const size_t local_length = vector->getLocalLength(); - - if (across) - for (unsigned int i = 0; i < local_length; ++i) - out << vector_1d(i) << ' '; - else - for (unsigned int i = 0; i < local_length; ++i) - out << vector_1d(i) << std::endl; - out << std::endl; - - // restore the representation - // of the vector - AssertThrow(out, ExcIO()); - } - - - - template - std::size_t - Vector::memory_consumption() const - { - return sizeof(*this) + - vector->getLocalLength() * - (sizeof(Number) + sizeof(TrilinosWrappers::types::int_type)); - } - - - - template - void - Vector::create_tpetra_comm_pattern(const IndexSet &source_index_set, - const MPI_Comm &mpi_comm) - { - source_stored_elements = source_index_set; - tpetra_comm_pattern = - std::make_shared( - locally_owned_elements(), source_index_set, mpi_comm); - } - template class Vector; template class Vector; # ifdef DEAL_II_WITH_COMPLEX_VALUES @@ -679,5 +36,4 @@ namespace LinearAlgebra DEAL_II_NAMESPACE_CLOSE # endif - #endif -- 2.39.5