From 3d57ef7106d83ebc9cc6e2a50cd8654fb13972e5 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Fri, 10 Feb 2017 09:26:33 -0500 Subject: [PATCH] Augment VectorSpaceVector interface. --- include/deal.II/lac/cuda_vector.h | 13 +++++ .../deal.II/lac/la_parallel_block_vector.h | 7 +++ .../lac/la_parallel_block_vector.templates.h | 12 +++++ include/deal.II/lac/la_parallel_vector.h | 7 +++ .../lac/la_parallel_vector.templates.h | 14 ++++++ include/deal.II/lac/la_vector.h | 43 ++++++++++++++++ include/deal.II/lac/la_vector.templates.h | 49 ++++++++++++++++++- include/deal.II/lac/read_write_vector.h | 17 ++----- include/deal.II/lac/trilinos_epetra_vector.h | 39 ++++++++++++++- include/deal.II/lac/trilinos_sparse_matrix.h | 1 + include/deal.II/lac/vector_space_vector.h | 12 +++++ source/lac/cuda_vector.cu | 25 +++++++++- source/lac/trilinos_epetra_vector.cc | 27 ++++++++++ tests/cuda/cuda_vector_01.cc | 8 +++ 14 files changed, 257 insertions(+), 17 deletions(-) diff --git a/include/deal.II/lac/cuda_vector.h b/include/deal.II/lac/cuda_vector.h index 17b57dc28d..c10ac4129f 100644 --- a/include/deal.II/lac/cuda_vector.h +++ b/include/deal.II/lac/cuda_vector.h @@ -83,6 +83,13 @@ namespace LinearAlgebra void reinit(const size_type n, const bool omit_zeroing_entries = false); + /** + * Change the dimension to that of the vector V. The elements of V are not + * copied. + */ + virtual void reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries = false) override; + /** * Import all the element from the input vector @p V. * VectorOperation::values @p operation is used to decide if the @@ -96,6 +103,12 @@ namespace LinearAlgebra std::shared_ptr communication_pattern = std::shared_ptr ()) override; + /** + * Sets all elements of the vector to the scalar @p s. This operation is + * only allowed if @p s is equal to zero. + */ + virtual Vector &operator= (const Number s) override; + /** * Multiply the entive vector by a fixed factor. */ diff --git a/include/deal.II/lac/la_parallel_block_vector.h b/include/deal.II/lac/la_parallel_block_vector.h index 7bc0135691..7dbccdf7bd 100644 --- a/include/deal.II/lac/la_parallel_block_vector.h +++ b/include/deal.II/lac/la_parallel_block_vector.h @@ -405,6 +405,13 @@ namespace LinearAlgebra */ //@{ + /** + * Change the dimension to that of the vector V. The elements of V are not + * copied. + */ + virtual void reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries = false); + /** * Multiply the entire vector by a fixed factor. */ diff --git a/include/deal.II/lac/la_parallel_block_vector.templates.h b/include/deal.II/lac/la_parallel_block_vector.templates.h index cdce69c73a..bb7e08add4 100644 --- a/include/deal.II/lac/la_parallel_block_vector.templates.h +++ b/include/deal.II/lac/la_parallel_block_vector.templates.h @@ -290,6 +290,18 @@ namespace LinearAlgebra + template + void + BlockVector::reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries) + { + Assert(dynamic_cast *>(&V)!=NULL, + ExcVectorTypeNotCompatible()); + const BlockVector &down_V = dynamic_cast &>(V); + reinit(down_V, omit_zeroing_entries); + } + + template BlockVector & BlockVector::operator *= (const Number factor) diff --git a/include/deal.II/lac/la_parallel_vector.h b/include/deal.II/lac/la_parallel_vector.h index f75a71332e..5c135c9119 100644 --- a/include/deal.II/lac/la_parallel_vector.h +++ b/include/deal.II/lac/la_parallel_vector.h @@ -534,6 +534,13 @@ namespace LinearAlgebra */ //@{ + /** + * Change the dimension to that of the vector V. The elements of V are not + * copied. + */ + virtual void reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries = false); + /** * Multiply the entire vector by a fixed factor. */ diff --git a/include/deal.II/lac/la_parallel_vector.templates.h b/include/deal.II/lac/la_parallel_vector.templates.h index fedc272467..a10d0e8a33 100644 --- a/include/deal.II/lac/la_parallel_vector.templates.h +++ b/include/deal.II/lac/la_parallel_vector.templates.h @@ -939,6 +939,20 @@ namespace LinearAlgebra + template + void Vector::reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries) + { + // Downcast. Throws an exception if invalid. + Assert(dynamic_cast *>(&V)!=NULL, + ExcVectorTypeNotCompatible()); + const Vector &down_V = dynamic_cast &>(V); + + reinit(down_V, omit_zeroing_entries); + } + + + template Vector & Vector::operator += (const VectorSpaceVector &vv) diff --git a/include/deal.II/lac/la_vector.h b/include/deal.II/lac/la_vector.h index fb64b714e8..609207b07e 100644 --- a/include/deal.II/lac/la_vector.h +++ b/include/deal.II/lac/la_vector.h @@ -87,6 +87,49 @@ namespace LinearAlgebra template Vector(const InputIterator first, const InputIterator last); + /** + * Set the global size of the vector to @p size. The stored elements have + * their index in [0,size). + * + * If the flag @p omit_zeroing_entries is set to false, the memory will be + * initialized with zero, otherwise the memory will be untouched (and the + * user must make sure to fill it with reasonable data before using it). + */ + virtual void reinit(const size_type size, + const bool omit_zeroing_entries = false); + + /** + * Uses the same IndexSet as the one of the input vector @p in_vector and + * allocates memory for this vector. + * + * If the flag @p omit_zeroing_entries is set to false, the memory will be + * initialized with zero, otherwise the memory will be untouched (and the + * user must make sure to fill it with reasonable data before using it). + */ + template + void reinit(const ReadWriteVector &in_vector, + const bool omit_zeroing_entries = false); + + /** + * Initializes the vector. The indices are specified by @p + * locally_stored_indices. + * + * If the flag @p omit_zeroing_entries is set to false, the memory will be + * initialized with zero, otherwise the memory will be untouched (and the + * user must make sure to fill it with reasonable data before using it). + * locally_stored_indices. + */ + virtual void reinit(const IndexSet &locally_stored_indices, + const bool omit_zeroing_entries = false); + + + /** + * Change the dimension to that of the vector V. The elements of V are not + * copied. + */ + virtual void reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries = false); + /** * Copies the data of the input vector @p in_vector. */ diff --git a/include/deal.II/lac/la_vector.templates.h b/include/deal.II/lac/la_vector.templates.h index 84ee087754..c823b74777 100644 --- a/include/deal.II/lac/la_vector.templates.h +++ b/include/deal.II/lac/la_vector.templates.h @@ -25,6 +25,51 @@ DEAL_II_NAMESPACE_OPEN namespace LinearAlgebra { + template + void Vector::reinit(const size_type size, + const bool omit_zeroing_entries) + { + ReadWriteVector::reinit(size, omit_zeroing_entries); + } + + + + template + template + void Vector::reinit(const ReadWriteVector &in_vector, + const bool omit_zeroing_entries) + { + ReadWriteVector::reinit(in_vector, omit_zeroing_entries); + } + + + + template + void Vector::reinit(const IndexSet &locally_stored_indices, + const bool omit_zeroing_entries) + { + ReadWriteVector::reinit(locally_stored_indices, omit_zeroing_entries); + } + + + + template + void Vector::reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries) + { + // Check that casting will work. + Assert(dynamic_cast*>(&V)!=NULL, ExcVectorTypeNotCompatible()); + + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast&>(V); + Assert(down_V.size()==this->size(), + ExcMessage("Cannot add two vectors with different numbers of elements")); + + ReadWriteVector::reinit(down_V, omit_zeroing_entries); + } + + + template Vector &Vector::operator= (const Vector &in_vector) { @@ -49,7 +94,7 @@ namespace LinearAlgebra { this->thread_loop_partitioner = in_vector.thread_loop_partitioner; if (this->size() != in_vector.size()) - this->reinit(in_vector, true); + ReadWriteVector::reinit(in_vector, true); dealii::internal::VectorOperations::Vector_copy copier(in_vector.val, this->val); internal::VectorOperations::parallel_for(copier, 0, this->size(), this->thread_loop_partitioner); @@ -446,7 +491,7 @@ namespace LinearAlgebra in.getline(buf,16,'\n'); sz = std::atoi(buf); - this->reinit(sz,true); + ReadWriteVector::reinit(sz,true); char c; // in >> c; diff --git a/include/deal.II/lac/read_write_vector.h b/include/deal.II/lac/read_write_vector.h index 2f38a41ea4..4ae7b29fd3 100644 --- a/include/deal.II/lac/read_write_vector.h +++ b/include/deal.II/lac/read_write_vector.h @@ -31,6 +31,7 @@ #ifdef DEAL_II_WITH_TRILINOS #include +#include #include "Epetra_MultiVector.h" #endif @@ -63,14 +64,6 @@ namespace TrilinosWrappers class Vector; } } - -namespace LinearAlgebra -{ - namespace EpetraWrappers - { - class Vector; - } -} #endif #ifdef DEAL_II_WITH_CUDA @@ -182,8 +175,8 @@ namespace LinearAlgebra * initialized with zero, otherwise the memory will be untouched (and the * user must make sure to fill it with reasonable data before using it). */ - void reinit (const size_type size, - const bool omit_zeroing_entries = false); + virtual void reinit (const size_type size, + const bool omit_zeroing_entries = false); /** * Uses the same IndexSet as the one of the input vector @p in_vector and @@ -206,8 +199,8 @@ namespace LinearAlgebra * user must make sure to fill it with reasonable data before using it). * locally_stored_indices. */ - void reinit (const IndexSet &locally_stored_indices, - const bool omit_zeroing_entries = false); + virtual void reinit (const IndexSet &locally_stored_indices, + const bool omit_zeroing_entries = false); #ifdef DEAL_II_WITH_CXX11 /** diff --git a/include/deal.II/lac/trilinos_epetra_vector.h b/include/deal.II/lac/trilinos_epetra_vector.h index c08d843ca6..1745f94768 100644 --- a/include/deal.II/lac/trilinos_epetra_vector.h +++ b/include/deal.II/lac/trilinos_epetra_vector.h @@ -24,7 +24,8 @@ #ifdef DEAL_II_WITH_MPI #include -#include +#include +#include #include #include #include "mpi.h" @@ -50,7 +51,7 @@ namespace LinearAlgebra * @ingroup Vectors * @author Bruno Turcksin, 2015 */ - class Vector : public VectorSpaceVector + class Vector : public VectorSpaceVector, public Subscriptor { public: /** @@ -83,6 +84,13 @@ namespace LinearAlgebra const MPI_Comm &communicator, const bool omit_zeroing_entries = false); + /** + * Change the dimension to that of the vector V. The elements of V are not + * copied. + */ + virtual void reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries = false); + /** * Copy function. This function takes a Vector and copies all the * elements. The Vector will have the same parallel distribution as @p @@ -90,6 +98,12 @@ namespace LinearAlgebra */ Vector &operator= (const Vector &V); + /** + * Sets all elements of the vector to the scalar @p s. This operation is + * only allowed if @p s is equal to zero. + */ + Vector &operator= (const double s); + /** * Imports all the elements present in the vector's IndexSet from the input * vector @p V. VectorOperation::values @p operation is used to decide if @@ -211,6 +225,11 @@ namespace LinearAlgebra virtual double add_and_dot(const double a, const VectorSpaceVector &V, const VectorSpaceVector &W); + /** + * This function always returns false and is present only for backward + * compatibily. + */ + bool has_ghost_elements() const; /** * Return the global size of the vector, equal to the sum of the number of @@ -308,9 +327,25 @@ namespace LinearAlgebra */ std::shared_ptr epetra_comm_pattern; }; + + + inline + bool Vector::has_ghost_elements() const + { + return false; + } } } + +/** + * Declare dealii::LinearAlgebra::EpetraWrappers::Vector as distributed vector. + */ +template <> +struct is_serial_vector : std_cxx11::false_type +{ +}; + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/include/deal.II/lac/trilinos_sparse_matrix.h b/include/deal.II/lac/trilinos_sparse_matrix.h index 15c53e6f60..36cb54b780 100644 --- a/include/deal.II/lac/trilinos_sparse_matrix.h +++ b/include/deal.II/lac/trilinos_sparse_matrix.h @@ -26,6 +26,7 @@ # include # include # include +# include # include #ifdef DEAL_II_WITH_CXX11 diff --git a/include/deal.II/lac/vector_space_vector.h b/include/deal.II/lac/vector_space_vector.h index ee966cc76d..8a18ff9f65 100644 --- a/include/deal.II/lac/vector_space_vector.h +++ b/include/deal.II/lac/vector_space_vector.h @@ -55,6 +55,18 @@ namespace LinearAlgebra typedef types::global_dof_index size_type; typedef typename numbers::NumberTraits::real_type real_type; + /** + * Change the dimension to that of the vector V. The elements of V are not + * copied. + */ + virtual void reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries = false) = 0; + + /** + * Sets all elements of the vector to the scalar @p s. This operation is + * only allowed if @p s is equal to zero. + */ + virtual VectorSpaceVector &operator= (const Number s) = 0; /** * Multiply the entire vector by a fixed factor. diff --git a/source/lac/cuda_vector.cu b/source/lac/cuda_vector.cu index 64e60a6a8c..4d0921f2f0 100644 --- a/source/lac/cuda_vector.cu +++ b/source/lac/cuda_vector.cu @@ -564,6 +564,15 @@ namespace LinearAlgebra + template + void Vector::reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries) + { + reinit(V.size(), omit_zeroing_entries); + } + + + template void Vector::import(const ReadWriteVector &V, VectorOperation::values operation, @@ -606,6 +615,20 @@ namespace LinearAlgebra + template + Vector &Vector::operator= (const Number s) + { + Assert(s == Number(), ExcMessage("Onlyt 0 can be assigned to a vector.")); + (void)s; + + cudaError_t error_code = cudaMemset(val, 0, n_elements*sizeof(Number)); + AssertCuda(error_code); + + return *this; + } + + + template Vector &Vector::operator*= (const Number factor) { @@ -1029,7 +1052,7 @@ namespace LinearAlgebra void Vector::print(std::ostream &out, const unsigned int precision, const bool scientific, - const bool across) const + const bool ) const { AssertThrow(out, ExcIO()); std::ios::fmtflags old_flags = out.flags(); diff --git a/source/lac/trilinos_epetra_vector.cc b/source/lac/trilinos_epetra_vector.cc index 078c36acf3..e60297315a 100644 --- a/source/lac/trilinos_epetra_vector.cc +++ b/source/lac/trilinos_epetra_vector.cc @@ -23,6 +23,7 @@ #include +#include #include "Epetra_Import.h" #include "Epetra_Map.h" #include "Epetra_MpiComm.h" @@ -43,6 +44,7 @@ namespace LinearAlgebra Vector::Vector(const Vector &V) : + Subscriptor(), vector(new Epetra_FEVector(V.trilinos_vector())) {} @@ -73,6 +75,21 @@ namespace LinearAlgebra + void Vector::reinit(const VectorSpaceVector &V, + const bool omit_zeroing_entries) + { + // Check that casting will work. + Assert(dynamic_cast(&V)!=NULL, 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); + } + + + Vector &Vector::operator= (const Vector &V) { // Distinguish three cases: @@ -100,6 +117,16 @@ namespace LinearAlgebra + Vector &Vector::operator= (const double s) + { + Assert(s==0., ExcMessage("Only 0 can be assigned to a vector.")); + vector->PutScalar(s); + + return *this; + } + + + void Vector::import(const ReadWriteVector &V, VectorOperation::values operation, std::shared_ptr communication_pattern) diff --git a/tests/cuda/cuda_vector_01.cc b/tests/cuda/cuda_vector_01.cc index f3c8b53685..2632adc84a 100644 --- a/tests/cuda/cuda_vector_01.cc +++ b/tests/cuda/cuda_vector_01.cc @@ -29,10 +29,13 @@ void test() LinearAlgebra::CUDAWrappers::Vector a; LinearAlgebra::CUDAWrappers::Vector b(size); LinearAlgebra::CUDAWrappers::Vector c(b); + LinearAlgebra::CUDAWrappers::Vector d; + d.reinit(c); AssertThrow(a.size()==0, ExcMessage("Vector has the wrong size.")); AssertThrow(b.size()==size, ExcMessage("Vector has the wrong size.")); AssertThrow(c.size()==size, ExcMessage("Vector has the wrong size.")); + AssertThrow(d.size()==size, ExcMessage("Vector has the wrong size.")); a.reinit(size); AssertThrow(a.size()==size, ExcMessage("Vector has the wrong size.")); @@ -95,6 +98,11 @@ void test() c.import(read_write_1, VectorOperation::insert); const double val = b*c; AssertThrow(val==328350., ExcMessage("Problem in operator *.")); + + b = 0.; + read_write_3.import(b, VectorOperation::insert); + for (unsigned int i=0; i