From b94eefc4f958f897b1b2d2ca5eb1cf0e45542759 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Wed, 28 Oct 2015 15:49:26 -0500 Subject: [PATCH] Add exceptions when we cannot downcast a VectorSpaceVector object to Vector. --- include/deal.II/lac/la_vector.h | 89 +++++++++++++---------- include/deal.II/lac/la_vector.templates.h | 52 ++++++++++--- include/deal.II/lac/vector_space_vector.h | 16 ++-- 3 files changed, 101 insertions(+), 56 deletions(-) diff --git a/include/deal.II/lac/la_vector.h b/include/deal.II/lac/la_vector.h index 20d972f2a5..0e0ea0acb7 100644 --- a/include/deal.II/lac/la_vector.h +++ b/include/deal.II/lac/la_vector.h @@ -51,6 +51,8 @@ namespace LinearAlgebra class Vector : public ReadWriteVector, public VectorSpaceVector { public: + typedef types::global_dof_index size_type; + /** * Constructor. Create a vector of dimension zero. */ @@ -72,11 +74,11 @@ namespace LinearAlgebra * v=Vector@(0);, i.e. the vector is replaced by one of * length zero. */ - explicit Vector(const typename ReadWriteVector::size_type n); + explicit Vector(const size_type n); /** * Initialize the vector with a given range of values pointed to by the - * iterators. This function is there in analogy to the @p std::vector class. + * iterators. This function exists in analogy to the @p std::vector class. */ template Vector(const InputIterator first, const InputIterator last); @@ -84,88 +86,88 @@ namespace LinearAlgebra /** * Destructor, deallocates memory. */ - ~Vector(); + virtual ~Vector(); /** * Multiply the entire vector by a fixed factor. */ - VectorSpaceVector &operator*= (const Number factor); + virtual VectorSpaceVector &operator*= (const Number factor); /** * Divide the entire vector by a fixed factor. */ - VectorSpaceVector &operator/= (const Number factor); + virtual VectorSpaceVector &operator/= (const Number factor); /** * Add the vector @p V to the present one. */ - VectorSpaceVector &operator+= (const VectorSpaceVector &V); + virtual VectorSpaceVector &operator+= (const VectorSpaceVector &V); /** * Substract the vector @p V from the present one. */ - VectorSpaceVector &operator-= (const VectorSpaceVector &V); + virtual VectorSpaceVector &operator-= (const VectorSpaceVector &V); /** * Return the scalar product of two vectors. */ - Number operator* (const VectorSpaceVector &V); + virtual Number operator* (const VectorSpaceVector &V); /** * Simple addition of a multiple of a vector, i.e. *this += a*V. */ - void add(const Number a, const VectorSpaceVector &V); + virtual void add(const Number a, const VectorSpaceVector &V); /** * Multiple addition of a multiple of a vector, i.e. *this += a*V+b*W. */ - void add(const Number a, const VectorSpaceVector &V, - const Number b, const VectorSpaceVector &W); + virtual void add(const Number a, const VectorSpaceVector &V, + const Number b, const VectorSpaceVector &W); /** * Scaling and simple vector addition, i.e. *this = s*(*this)+V. */ - void sadd(const Number s, const VectorSpaceVector &V); + virtual void sadd(const Number s, const VectorSpaceVector &V); /** * Scaling and simple addition of a multiple of a vector, i.e. *this = * s*(*this)+a*V. */ - void sadd(const Number s, const Number a, - const VectorSpaceVector &V); + virtual void sadd(const Number s, const Number a, + const VectorSpaceVector &V); /** * Scale each element of this vector by the corresponding element in the * argument. This function is mostly meant to simulate multiplication (and * immediate re-assignement) by a diagonal scaling matrix. */ - void scale(const VectorSpaceVector &scaling_factors); + virtual void scale(const VectorSpaceVector &scaling_factors); /** * Assignement *this = a*V. */ - void equ(const Number a, const VectorSpaceVector &V); + virtual void equ(const Number a, const VectorSpaceVector &V); /** - * Returns the l1 norm of the vector (i.e., the sum of the + * Return the l1 norm of the vector (i.e., the sum of the * absolute values of all entries). */ - typename VectorSpaceVector::real_type l1_norm(); + virtual typename VectorSpaceVector::real_type l1_norm(); /** - * Returns the l2 norm of the vector (i.e., the square root of + * Return the l2 norm of the vector (i.e., the square root of * the sum of the square of all entries among all processors). */ - typename VectorSpaceVector::real_type l2_norm(); + virtual typename VectorSpaceVector::real_type l2_norm(); /** - * Returns the maximum norm of the vector (i.e., the maximum absolute + * Return the maximum norm of the vector (i.e., the maximum absolute * value among all entries and among all processors). */ - typename VectorSpaceVector::real_type linfty_norm(); + virtual typename VectorSpaceVector::real_type linfty_norm(); /** - * Performs a combined operation of a vector addition and a subsequent + * Perform a combined operation of a vector addition and a subsequent * inner product, returning the value of the inner product. In other * words, the result of this function is the same as if the user called * @code @@ -173,18 +175,18 @@ namespace LinearAlgebra * return_value = *this * W; * @endcode */ - Number add_and_dot(const Number a, - const VectorSpaceVector &V, - const VectorSpaceVector &W); + virtual Number add_and_dot(const Number a, + const VectorSpaceVector &V, + const VectorSpaceVector &W); /** - * Returns the global size of the vector, equal to the sum of the number + * Return the global size of the vector, equal to the sum of the number * of locally owned indices among all processors. */ - typename ReadWriteVector::size_type size() const; + size_type size() const; /** - * Returns an index set that describes which elements of this vector are + * Return an index set that describes which elements of this vector are * owned by the current processor. As a consequence, the index sets returned * on different procesors if this is a distributed vector will form disjoint * sets that add up to the complete index set. Obviously, if a vector is @@ -198,16 +200,23 @@ namespace LinearAlgebra /** * Prints the vector to the output stream @p out. */ - void print(std::ostream &out, - const unsigned int precision=3, - const bool scientific=true, - const bool across=true) const; + virtual void print(std::ostream &out, + const unsigned int precision=3, + const bool scientific=true, + const bool across=true) const; /** * Returns the memory consumption of this class in bytes. */ std::size_t memory_consumption() const; + /** + * Attempt to perform an operation between two incompatible vector types. + * + * @ingroup Exceptions + */ + DeclException0(ExcVectorTypeNotCompatible); + private: /** * Compute the L1 norm in a recursive way by dividing the vector on smaller @@ -217,11 +226,13 @@ namespace LinearAlgebra unsigned int j); /** - * Compute the L2 norm in a recursive way by dividing the vector on smaller - * and smaller intervals. This reduces the numerical error on large vector. + * Compute the squared L2 norm in a recursive way by dividing the vector on + * smaller and smaller intervals. This reduces the numerical error on large + * vector. */ - typename VectorSpaceVector::real_type l2_norm_recursive(unsigned int i, - unsigned int j); + typename VectorSpaceVector::real_type l2_norm_squared_recursive( + unsigned int i, + unsigned int j); }; /*@}*/ @@ -244,7 +255,7 @@ namespace LinearAlgebra template inline - Vector::Vector(const typename ReadWriteVector::size_type n) + Vector::Vector(const size_type n) : ReadWriteVector(n) {} @@ -270,7 +281,7 @@ namespace LinearAlgebra template inline - typename ReadWriteVector::size_type Vector::size() const + typename Vector::size_type Vector::size() const { return ReadWriteVector::size(); } diff --git a/include/deal.II/lac/la_vector.templates.h b/include/deal.II/lac/la_vector.templates.h index 9057e54249..f436504d48 100644 --- a/include/deal.II/lac/la_vector.templates.h +++ b/include/deal.II/lac/la_vector.templates.h @@ -49,6 +49,9 @@ namespace LinearAlgebra template VectorSpaceVector &Vector::operator+= (const VectorSpaceVector &V) { + // 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(), @@ -67,6 +70,9 @@ namespace LinearAlgebra template VectorSpaceVector &Vector::operator-= (const VectorSpaceVector &V) { + // 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(), @@ -85,6 +91,10 @@ namespace LinearAlgebra template Number Vector::operator* (const VectorSpaceVector &V) { + // 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(), @@ -105,6 +115,10 @@ namespace LinearAlgebra template void Vector::add(const Number a, const VectorSpaceVector &V) { + // 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); AssertIsFinite(a); @@ -123,6 +137,13 @@ namespace LinearAlgebra void Vector::add(const Number a, const VectorSpaceVector &V, const Number b, const VectorSpaceVector &W) { + // Check that casting will work. + Assert(dynamic_cast*>(&V)!=NULL, + ExcVectorTypeNotCompatible()); + // Check that casting will work. + Assert(dynamic_cast*>(&W)!=NULL, + ExcVectorTypeNotCompatible()); + // Downcast V. If fails, throws an exception. const Vector &down_V = dynamic_cast&>(V); // Downcast W. If fails, throws an exception. @@ -156,6 +177,10 @@ namespace LinearAlgebra void Vector::sadd(const Number s, const Number a, const VectorSpaceVector &V) { + // Check that casting will work. + Assert(dynamic_cast*>(&V)!=NULL, + ExcVectorTypeNotCompatible()); + *this *= s; // Downcast V. It fails, throws an exception. const Vector &down_V = dynamic_cast&>(V); @@ -169,6 +194,10 @@ namespace LinearAlgebra template void Vector::scale(const VectorSpaceVector &scaling_factors) { + // Check that casting will work. + Assert(dynamic_cast*>(&scaling_factors)!=NULL, + ExcVectorTypeNotCompatible()); + // Downcast scaling_factors. If fails, throws an exception. const Vector &down_scaling_factors = dynamic_cast&>(scaling_factors); @@ -186,6 +215,10 @@ namespace LinearAlgebra template void Vector::equ(const Number a, const VectorSpaceVector &V) { + // 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); *this = down_V; @@ -220,7 +253,7 @@ namespace LinearAlgebra // precision) using the BLAS approach with a weight, see e.g. dnrm2.f. typedef typename ReadWriteVector::real_type real_type; real_type norm_square = 0.; - norm_square = l2_norm_recursive(0,this->size()-1); + norm_square = l2_norm_squared_recursive(0,this->size()-1); if (numbers::is_finite(norm_square) && norm_square>=std::numeric_limits::min()) return std::sqrt(norm_square); @@ -255,8 +288,7 @@ namespace LinearAlgebra { typename ReadWriteVector::real_type norm = 0.; for (unsigned int i=0; isize(); ++i) - if (std::abs(this->val[i])>norm) - norm = std::abs(this->val[i]); + norm = std::max(std::abs(this->val[i]),norm); return norm; } @@ -276,8 +308,9 @@ namespace LinearAlgebra template - typename VectorSpaceVector::real_type Vector::l1_norm_recursive(unsigned int i, - unsigned int j) + typename VectorSpaceVector::real_type Vector::l1_norm_recursive( + unsigned int i, + unsigned int j) { Assert(j>=i, ExcInternalError()); typename ReadWriteVector::real_type norm = 0.; @@ -296,16 +329,17 @@ namespace LinearAlgebra template - typename VectorSpaceVector::real_type Vector::l2_norm_recursive(unsigned int i, - unsigned int j) + typename VectorSpaceVector::real_type Vector::l2_norm_squared_recursive( + unsigned int i, + unsigned int j) { Assert(j>=i, ExcInternalError()); typename ReadWriteVector::real_type norm = 0.; if ((j-i)!=0) { - norm += l2_norm_recursive(i,(i+j)/2); - norm += l2_norm_recursive((i+j)/2+1,j); + norm += l2_norm_squared_recursive(i,(i+j)/2); + norm += l2_norm_squared_recursive((i+j)/2+1,j); } else norm += std::pow(std::abs(this->val[i]),2); diff --git a/include/deal.II/lac/vector_space_vector.h b/include/deal.II/lac/vector_space_vector.h index 60e47ee771..31ba1b1c72 100644 --- a/include/deal.II/lac/vector_space_vector.h +++ b/include/deal.II/lac/vector_space_vector.h @@ -113,25 +113,25 @@ namespace LinearAlgebra virtual void equ(const Number a, const VectorSpaceVector &V) = 0; /** - * Returns the l1 norm of the vector (i.e., the sum of the + * Return the l1 norm of the vector (i.e., the sum of the * absolute values of all entries among all processors). */ virtual real_type l1_norm() = 0; /** - * Returns the l2 norm of the vector (i.e., the square root of + * Return the l2 norm of the vector (i.e., the square root of * the sum of the square of all entries among all processors). */ virtual real_type l2_norm() = 0; /** - * Returns the maximum norm of the vector (i.e., the maximum absolute value + * Return the maximum norm of the vector (i.e., the maximum absolute value * among all entries and among all processors). */ virtual real_type linfty_norm() = 0; /** - * Performs a combined operation of a vector addition and a subsequent + * Perform a combined operation of a vector addition and a subsequent * inner product, returning the value of the inner product. In other * words, the result of this function is the same as if the user called * @code @@ -151,13 +151,13 @@ namespace LinearAlgebra const VectorSpaceVector &W) = 0; /** - * Returns the global size of the vector, equal to the sum of the number of + * Return the global size of the vector, equal to the sum of the number of * locally owned indices among all processors. */ virtual size_type size() const = 0; /** - * Returns an index set that describes which elements of this vector are + * Return an index set that describes which elements of this vector are * owned by the current processor. As a consequence, the index sets returned * on different procesors if this is a distributed vector will form disjoint * sets that add up to the complete index set. Obviously, if a vector is @@ -169,7 +169,7 @@ namespace LinearAlgebra virtual const dealii::IndexSet &locally_owned_elements() const = 0; /** - * Prints the vector to the output stream @p out. + * Print the vector to the output stream @p out. */ virtual void print(std::ostream &out, const unsigned int precision=3, @@ -177,7 +177,7 @@ namespace LinearAlgebra const bool across=true) const = 0; /** - * Returns the memory consumption of this class in bytes. + * Return the memory consumption of this class in bytes. */ virtual std::size_t memory_consumption() const = 0; }; -- 2.39.5