From: Bruno Turcksin Date: Tue, 27 Oct 2015 14:19:04 +0000 (-0500) Subject: Add LinearAlgebra::Vector derived from LinearAlgebra::ReadWriteVector and LinearAlgeb... X-Git-Tag: v8.4.0-rc2~215^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b5229bec209bbe9dd315c50e98ad34519413f968;p=dealii.git Add LinearAlgebra::Vector derived from LinearAlgebra::ReadWriteVector and LinearAlgebra::VectorSpaceVector. --- diff --git a/include/deal.II/lac/la_vector.h b/include/deal.II/lac/la_vector.h new file mode 100644 index 0000000000..20d972f2a5 --- /dev/null +++ b/include/deal.II/lac/la_vector.h @@ -0,0 +1,311 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii__la_vector_h +#define dealii__la_vector_h + + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + /*! @addtogroup Vectors + *@{ + */ + + /** + * Numerical vector of data. This class derives from both + * ::dealii::LinearAlgebra::ReadWriteVector and + * ::dealii::LinearAlgebra::VectorSpaceVector. As opposed to the array of the + * C++ standard library, this class implements an element of a vector space + * suitable for numerical computations. + * + * @authod Bruno Turcksin, 2015. + */ + template + class Vector : public ReadWriteVector, public VectorSpaceVector + { + public: + /** + * Constructor. Create a vector of dimension zero. + */ + Vector(); + + /** + * Copy constructor. Sets the dimension to that of the given vector and + * copies all elements. + */ + Vector(const Vector &V); + + /** + * Constructor. Set dimension to @p n and initialize all elements with + * zero. + * + * The constructor is made explicit to avoid accident like this: + * v=0;. Presumably, the user wants to set every element of the + * vector to zero, but instead, what happens is this call: + * v=Vector@(0);, i.e. the vector is replaced by one of + * length zero. + */ + explicit Vector(const typename ReadWriteVector::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. + */ + template + Vector(const InputIterator first, const InputIterator last); + + /** + * Destructor, deallocates memory. + */ + ~Vector(); + + /** + * Multiply the entire vector by a fixed factor. + */ + VectorSpaceVector &operator*= (const Number factor); + + /** + * Divide the entire vector by a fixed factor. + */ + VectorSpaceVector &operator/= (const Number factor); + + /** + * Add the vector @p V to the present one. + */ + VectorSpaceVector &operator+= (const VectorSpaceVector &V); + + /** + * Substract the vector @p V from the present one. + */ + VectorSpaceVector &operator-= (const VectorSpaceVector &V); + + /** + * Return the scalar product of two vectors. + */ + 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); + + /** + * 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); + + /** + * Scaling and simple vector addition, i.e. *this = s*(*this)+V. + */ + 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); + + /** + * 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); + + /** + * Assignement *this = a*V. + */ + void equ(const Number a, const VectorSpaceVector &V); + + /** + * Returns the l1 norm of the vector (i.e., the sum of the + * absolute values of all entries). + */ + typename VectorSpaceVector::real_type l1_norm(); + + /** + * Returns 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(); + + /** + * Returns 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(); + + /** + * Performs 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 + * this->add(a, V); + * return_value = *this * W; + * @endcode + */ + 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 + * of locally owned indices among all processors. + */ + typename ReadWriteVector::size_type size() const; + + /** + * Returns 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 + * created on only one processor, then the result would satisfy + * @code + * vec.locally_owned_elements() == complete_index_set(vec.size()) + * @endcode + */ + const dealii::IndexSet &locally_owned_elements() const; + + /** + * 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; + + /** + * Returns the memory consumption of this class in bytes. + */ + std::size_t memory_consumption() const; + + private: + /** + * Compute the L1 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 l1_norm_recursive(unsigned int i, + 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. + */ + typename VectorSpaceVector::real_type l2_norm_recursive(unsigned int i, + unsigned int j); + }; + + /*@}*/ + /*----------------------- Inline functions ----------------------------------*/ + + template + inline + Vector::Vector() {} + + + + template + inline + Vector::Vector(const Vector &V) + : + ReadWriteVector(V) + {} + + + + template + inline + Vector::Vector(const typename ReadWriteVector::size_type n) + : + ReadWriteVector(n) + {} + + + + template + template + inline + Vector::Vector(const InputIterator first, const InputIterator last) + { + this->reinit(complete_index_set(std::distance (first, last)), true); + std::copy(first, last, this->begin()); + } + + + + template + inline + Vector::~Vector() {} + + + + template + inline + typename ReadWriteVector::size_type Vector::size() const + { + return ReadWriteVector::size(); + } + + + + template + inline + const dealii::IndexSet &Vector::locally_owned_elements() const + { + return ReadWriteVector::get_stored_elements(); + } + + + + template + inline + void Vector::print(std::ostream &out, + const unsigned int precision, + const bool scientific, + const bool) const + { + ReadWriteVector::print(out, precision, scientific); + } + + + + template + inline + std::size_t Vector::memory_consumption() const + { + return ReadWriteVector::memory_consumption(); + } +} // end of namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/lac/la_vector.templates.h b/include/deal.II/lac/la_vector.templates.h new file mode 100644 index 0000000000..9057e54249 --- /dev/null +++ b/include/deal.II/lac/la_vector.templates.h @@ -0,0 +1,319 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii__la_vector_templates_h +#define dealii__la_vector_templates_h + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + template + VectorSpaceVector &Vector::operator*= (const Number factor) + { + AssertIsFinite(factor); + for (unsigned int i=0; isize(); ++i) + this->val[i] *= factor; + + return *this; + } + + + + template + VectorSpaceVector &Vector::operator/= (const Number factor) + { + AssertIsFinite(factor); + Assert(factor!=Number(0.), ExcZero()); + this->operator *= (Number(1.)/factor); + + return *this; + } + + + + template + VectorSpaceVector &Vector::operator+= (const VectorSpaceVector &V) + { + // 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")); + for (unsigned int i=0; isize(); ++i) + { + AssertIsFinite(down_V[i]); + this->val[i] += down_V[i]; + } + + return *this; + } + + + + template + VectorSpaceVector &Vector::operator-= (const VectorSpaceVector &V) + { + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast&>(V); + Assert(down_V.size()==this->size(), + ExcMessage("Cannot subtract two vectors with different numbers of elements")); + for (unsigned int i=0; isize(); ++i) + { + AssertIsFinite(down_V[i]); + this->val[i] -= down_V[i]; + } + + return *this; + } + + + + template + Number Vector::operator* (const VectorSpaceVector &V) + { + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast&>(V); + Assert(down_V.size()==this->size(), + ExcMessage("Cannot compute the scalar product " + "of two vectors with different numbers of elements")); + Number value = 0.; + for (unsigned int i=0; isize(); ++i) + { + AssertIsFinite(down_V[i]); + value += this->val[i]*down_V[i]; + } + + return value; + } + + + + template + void Vector::add(const Number a, const VectorSpaceVector &V) + { + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast&>(V); + AssertIsFinite(a); + Assert(down_V.size()==this->size(), + ExcMessage("Cannot add two vectors with different numbers of elements")); + for (unsigned int i=0; isize(); ++i) + { + AssertIsFinite(down_V[i]); + this->val[i] += a*down_V[i]; + } + } + + + + template + void Vector::add(const Number a, const VectorSpaceVector &V, + const Number b, const VectorSpaceVector &W) + { + // 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); + AssertIsFinite(a); + Assert(down_V.size()==this->size(), + ExcMessage("Cannot add two vectors with different numbers of elements")); + AssertIsFinite(b); + Assert(down_W.size()==this->size(), + ExcMessage("Cannot add two vectors with different numbers of elements")); + for (unsigned int i=0; isize(); ++i) + { + AssertIsFinite(down_V[i]); + AssertIsFinite(down_W[i]); + this->val[i] += a*down_V[i]+b*down_W[i]; + } + } + + + + template + void Vector::sadd(const Number s, const VectorSpaceVector &V) + { + *this *= s; + *this += V; + } + + + + template + void Vector::sadd(const Number s, const Number a, + const VectorSpaceVector &V) + { + *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) + { + // Downcast scaling_factors. If fails, throws an exception. + const Vector &down_scaling_factors = + dynamic_cast&>(scaling_factors); + Assert(down_scaling_factors.size()==this->size(), + ExcMessage("Cannot add two vectors with different numbers of elements")); + for (unsigned int i=0; isize(); ++i) + { + AssertIsFinite(down_scaling_factors[i]); + this->val[i] *= down_scaling_factors[i]; + } + } + + + + template + void Vector::equ(const Number a, const VectorSpaceVector &V) + { + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast&>(V); + *this = down_V; + *this *= a; + } + + + + template + typename VectorSpaceVector::real_type Vector::l1_norm() + { + Assert (this->size(), ExcEmptyObject()); + + typename ReadWriteVector::real_type norm = 0.; + norm = l1_norm_recursive(0,this->size()-1); + + return norm; + } + + + + + template + typename VectorSpaceVector::real_type Vector::l2_norm() + { + Assert (this->size(), ExcEmptyObject()); + + // if l2_norm()^2 is finite and non-zero, the answer is computed as + // std::sqrt(norm_sqr()). If norm_sqrt() is infinite or zero, the l2 norm + // might still be finite. In that case, recompute ig (this is a rare case, + // so working on the vector twice is uncritical and paid off by the extended + // 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); + if (numbers::is_finite(norm_square) && + norm_square>=std::numeric_limits::min()) + return std::sqrt(norm_square); + else + { + real_type scale = 0.; + real_type sum = 1.; + for (unsigned int i=0; isize(); ++i) + { + if (this->val[i] != Number()) + { + const real_type abs_x = + numbers::NumberTraits::abs(this->val[i]); + if (scale < abs_x) + { + sum = 1. + sum * (scale/abs_x) * (scale/abs_x); + scale = abs_x; + } + else + sum += (abs_x/scale) * (abs_x/scale); + } + } + AssertIsFinite(scale*std::sqrt(sum)); + return scale*std::sqrt(sum); + } + } + + + + template + typename VectorSpaceVector::real_type Vector::linfty_norm() + { + 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]); + + return norm; + } + + + + template + Number Vector::add_and_dot(const Number a, + const VectorSpaceVector &V, + const VectorSpaceVector &W) + { + this->add(a,V); + + return *this*W; + } + + + + template + typename VectorSpaceVector::real_type Vector::l1_norm_recursive(unsigned int i, + unsigned int j) + { + Assert(j>=i, ExcInternalError()); + typename ReadWriteVector::real_type norm = 0.; + + if ((j-i)!=0) + { + norm += l1_norm_recursive(i,(i+j)/2); + norm += l1_norm_recursive((i+j)/2+1,j); + } + else + norm += std::abs(this->val[i]); + + return norm; + } + + + + template + typename VectorSpaceVector::real_type Vector::l2_norm_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); + } + else + norm += std::pow(std::abs(this->val[i]),2); + + return norm; + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/lac/vector_space_vector.h b/include/deal.II/lac/vector_space_vector.h index 47184f60e0..60e47ee771 100644 --- a/include/deal.II/lac/vector_space_vector.h +++ b/include/deal.II/lac/vector_space_vector.h @@ -20,14 +20,16 @@ #include #include + DEAL_II_NAMESPACE_OPEN +class IndexSet; + namespace LinearAlgebra { template class ReadWriteVector; - class IndexSet; /*! @addtogroup Vectors *@{ @@ -49,6 +51,7 @@ namespace LinearAlgebra typedef types::global_dof_index size_type; typedef typename numbers::NumberTraits::real_type real_type; + /** * Multiply the entire vector by a fixed factor. */ @@ -91,7 +94,7 @@ namespace LinearAlgebra virtual void sadd(const Number s, const VectorSpaceVector &V) = 0; /** - * Scaling and simple additiion of a multiple of a vector, i.e. *this = + * Scaling and simple addition of a multiple of a vector, i.e. *this = * s*(*this)+a*V. */ virtual void sadd(const Number s, const Number a, @@ -154,7 +157,7 @@ namespace LinearAlgebra virtual size_type size() const = 0; /** - * Return an index set that describes which elements of this vector are + * Returns 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 @@ -163,7 +166,7 @@ namespace LinearAlgebra * vec.locally_owned_elements() == complete_index_set(vec.size()) * @endcode */ - virtual IndexSet locally_owned_elements() const = 0; + virtual const dealii::IndexSet &locally_owned_elements() const = 0; /** * Prints the vector to the output stream @p out. @@ -178,7 +181,6 @@ namespace LinearAlgebra */ virtual std::size_t memory_consumption() const = 0; }; - /*@}*/ } diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index b4de9e64e3..648edcfc28 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -27,6 +27,7 @@ SET(_src constraint_matrix.cc full_matrix.cc lapack_full_matrix.cc + la_vector.cc matrix_lib.cc matrix_out.cc parallel_vector.cc @@ -60,6 +61,7 @@ SET(_inst constraint_matrix.inst.in full_matrix.inst.in lapack_full_matrix.inst.in + la_vector.inst.in parallel_vector.inst.in precondition_block.inst.in relaxation_block.inst.in diff --git a/source/lac/la_vector.cc b/source/lac/la_vector.cc new file mode 100644 index 0000000000..ea6500ae7f --- /dev/null +++ b/source/lac/la_vector.cc @@ -0,0 +1,25 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ +#include "la_vector.inst" +} + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/lac/la_vector.inst.in b/source/lac/la_vector.inst.in new file mode 100644 index 0000000000..5402820eb2 --- /dev/null +++ b/source/lac/la_vector.inst.in @@ -0,0 +1,29 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1999 - 2014 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. +// +// --------------------------------------------------------------------- + + + +for (SCALAR : REAL_SCALARS) + { + template class Vector; + } + + + +for (SCALAR : COMPLEX_SCALARS) + { + template class Vector; + } + diff --git a/source/lac/read_write_vector.inst.in b/source/lac/read_write_vector.inst.in index eceace6fd4..dd83ecdbca 100644 --- a/source/lac/read_write_vector.inst.in +++ b/source/lac/read_write_vector.inst.in @@ -32,3 +32,21 @@ for (S1, S2 : REAL_SCALARS) \} } + +for (SCALAR : COMPLEX_SCALARS) +{ + namespace LinearAlgebra + \{ + template class ReadWriteVector; + \} +} + +for (S1, S2 : COMPLEX_SCALARS) +{ + namespace LinearAlgebra + \{ + template void ReadWriteVector::reinit (const ReadWriteVector&, + const bool); + \} +} + diff --git a/tests/lac/la_vector_accumulation_01.cc b/tests/lac/la_vector_accumulation_01.cc new file mode 100644 index 0000000000..a4c00f806d --- /dev/null +++ b/tests/lac/la_vector_accumulation_01.cc @@ -0,0 +1,67 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2012 - 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. +// +// --------------------------------------------------------------------- + + +// check that the l2 norm is exactly the same for many runs on random vector +// size with random vector entries. this is to ensure that the result is +// reproducible also when parallel evaluation is done + +#include "../tests.h" +#include +#include +#include +#include +#include + + + + +template +void check_norms () +{ + for (unsigned int test=0; test<20; ++test) + { + const unsigned int size = Testing::rand() % 100000; + LinearAlgebra::Vector vec (size); + for (unsigned int i=0; i(Testing::rand())/static_cast(RAND_MAX); + const typename LinearAlgebra::ReadWriteVector::real_type norm = vec.l2_norm(); + for (unsigned int i=0; i<30; ++i) + AssertThrow (vec.l2_norm() == norm, ExcInternalError()); + + LinearAlgebra::Vector vec2 (vec); + for (unsigned int i=0; i<10; ++i) + AssertThrow (vec2.l2_norm() == norm, ExcInternalError()); + } +} + + +int main() +{ + std::ofstream logfile("output"); + deallog << std::fixed; + deallog << std::setprecision(2); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + check_norms(); + check_norms(); + check_norms(); + check_norms >(); + deallog << "OK" << std::endl; +} + + diff --git a/tests/lac/la_vector_accumulation_01.output b/tests/lac/la_vector_accumulation_01.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/lac/la_vector_accumulation_01.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/lac/la_vector_add_and_dot.cc b/tests/lac/la_vector_add_and_dot.cc new file mode 100644 index 0000000000..4441163888 --- /dev/null +++ b/tests/lac/la_vector_add_and_dot.cc @@ -0,0 +1,83 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2012 - 2014 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 that LinearAlgebra::Vector::add_and_dot works correctly + +#include "../tests.h" +#include +#include +#include +#include +#include + + + + +template +void check () +{ + for (unsigned int test=0; test<5; ++test) + { + const unsigned int size = 17 + test*1101; + LinearAlgebra::Vector v1 (size), v2(size), v3(size), check(size); + for (unsigned int i=0; i::value) + { + deallog << "Vector add reference: "; + for (unsigned int i=0; i(size) + << ", is " << prod_check/static_cast(size) + << std::endl; + } +} + + +int main() +{ + std::ofstream logfile("output"); + deallog << std::fixed; + deallog << std::setprecision(2); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + check(); + check(); + check(); + check >(); + deallog << "OK" << std::endl; +} + + diff --git a/tests/lac/la_vector_add_and_dot.output b/tests/lac/la_vector_add_and_dot.output new file mode 100644 index 0000000000..329e9f8d54 --- /dev/null +++ b/tests/lac/la_vector_add_and_dot.output @@ -0,0 +1,24 @@ + +DEAL::Add and dot should be 0.30, is 0.30 +DEAL::Add and dot should be 13.40, is 13.40 +DEAL::Add and dot should be 26.50, is 26.50 +DEAL::Add and dot should be 39.61, is 39.61 +DEAL::Add and dot should be 52.71, is 52.71 +DEAL::Vector add reference: 0.03 0.03 0.04 0.05 0.06 0.06 0.07 0.08 0.09 0.09 0.10 0.11 0.12 0.12 0.13 0.14 0.15 +DEAL::Vector check reference: 0.03 0.03 0.04 0.05 0.06 0.06 0.07 0.08 0.09 0.09 0.10 0.11 0.12 0.12 0.13 0.14 0.15 +DEAL::Add and dot should be 0.30, is 0.30 +DEAL::Add and dot should be 13.40, is 13.40 +DEAL::Add and dot should be 26.50, is 26.50 +DEAL::Add and dot should be 39.61, is 39.61 +DEAL::Add and dot should be 52.71, is 52.71 +DEAL::Add and dot should be 0.30, is 0.30 +DEAL::Add and dot should be 13.40, is 13.40 +DEAL::Add and dot should be 26.50, is 26.50 +DEAL::Add and dot should be 39.61, is 39.61 +DEAL::Add and dot should be 52.71, is 52.71 +DEAL::Add and dot should be (0.30,0.00), is (0.30,0.00) +DEAL::Add and dot should be (13.40,0.00), is (13.40,0.00) +DEAL::Add and dot should be (26.50,0.00), is (26.50,0.00) +DEAL::Add and dot should be (39.61,0.00), is (39.61,0.00) +DEAL::Add and dot should be (52.71,0.00), is (52.71,0.00) +DEAL::OK diff --git a/tests/lac/la_vector_large_numbers.cc b/tests/lac/la_vector_large_numbers.cc new file mode 100644 index 0000000000..884e053ea4 --- /dev/null +++ b/tests/lac/la_vector_large_numbers.cc @@ -0,0 +1,81 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2012 - 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. +// +// --------------------------------------------------------------------- + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include + + + +void check_large_numbers() +{ + LinearAlgebra::Vector v(10); + v(0) = 1e13; + v(1) = 1e19; + v(2) = 1e21; + v(4) = 1.; + v(7) = 1e20; + + const double correct = std::sqrt(1e26 + 1e38 + 1e42 + 1e40 + 1.); + AssertThrow (std::abs(v.l2_norm() - correct) < 1e-6*correct, ExcInternalError()); + + for (unsigned int i=0; i w(7); + w(0) = 1e232; + w(6) = 1e231; + w(3) = 1e200; + const double correct3 = std::sqrt(100. + 1.) * 1e231; + AssertThrow (std::abs(w.l2_norm() - correct3) < 1e-13*correct3, ExcInternalError()); + + for (unsigned int i=0; i +#include +#include +#include +#include + + + + +template +void check_norms () +{ + const number acc = 1e1*std::numeric_limits::epsilon(); + unsigned int skip = 73; + for (unsigned int size=1; size<200000; size+=skip) + { + // test correctness + if (size > 10000) + skip += 17; + LinearAlgebra::Vector vec(size); + for (unsigned int i=0; i +void check_complex_norms () +{ + const number acc = 1e1*std::numeric_limits::epsilon(); + unsigned int skip = 73; + for (unsigned int size=1; size<100000; size+=skip) + { + // test correctness + if (size > 10000) + skip += 17; + LinearAlgebra::Vector > vec(size); + long double sum = 0.; + for (unsigned int i=0; i(i+1,i+2); + sum += std::sqrt((long double)(i+1)*(i+1) + (long double)(i+2)*(i+2)); + } + + const number l1_norm = vec.l1_norm(); + AssertThrow (std::abs(l1_norm-sum) < acc*sum, + ExcInternalError()); + + // test accuracy of summation + const std::complex value (3.14159265358979323846, 0.1); + for (unsigned int i=0; i(value); + const number l1_norma = vec.l1_norm(); + AssertThrow (std::abs(l1_norma-std::abs(value)*size) < + acc*size*std::abs(value), + ExcInternalError()); + const number l2_norma = vec.l2_norm(); + AssertThrow (std::abs(l2_norma-std::abs(value)*std::sqrt((number)size)) < + acc*std::sqrt((number)size)*std::abs(value), + ExcInternalError()); + } +} + + +int main() +{ + std::ofstream logfile("output"); + deallog << std::fixed; + deallog << std::setprecision(2); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + check_norms(); + check_norms(); + check_norms(); + check_complex_norms(); + check_complex_norms(); + deallog << "OK" << std::endl; +} + + diff --git a/tests/lac/la_vector_norms.output b/tests/lac/la_vector_norms.output new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/lac/la_vector_norms.output @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/lac/la_vector_vector.cc b/tests/lac/la_vector_vector.cc new file mode 100644 index 0000000000..1dbc2c4be9 --- /dev/null +++ b/tests/lac/la_vector_vector.cc @@ -0,0 +1,159 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 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. +// +// --------------------------------------------------------------------- + + + +#include "../tests.h" +#include +#include +#include +#include +#include + + + + +const unsigned int N=10; +unsigned int check_point = 0; + + + + +template +void print (const LinearAlgebra::Vector &v) +{ + for (unsigned int i=0; i +void check_vectors (LinearAlgebra::Vector &d1, LinearAlgebra::Vector &d2) +{ + deallog << "Fill & Swap" << std::endl; + LinearAlgebra::Vector d3(d1.size()); + print (d3); + + for (unsigned int i=0; i d4 (&array[0], &array[4]); + print (d4); + + deallog << "Extract number" << std::endl; + // Each line should contain two equal numbers + double sum = 0.; + for (unsigned int i=0; isum) sum = t; + } + deallog << d3.linfty_norm() << '\t' << sum << std::endl; + + deallog << "add & sub" << std::endl; + + d1 += d2; + print (d1); + + d2 -= d1; + print (d2); + + d1.add (2, d3); + print (d1); + + d1.add (2., d2, .5, d3); + print (d1); + + deallog << "sadd & scale" << std::endl; + + d2.sadd (2., d1); + print (d2); + + d2.sadd (2., .5, d1); + print (d2); + + d1 *= 4.; + print (d1); + + deallog << "equ" << std::endl; + + d2.equ (.25, d1); + print (d2); +} + + +int main() +{ + std::ofstream logfile("output"); + deallog << std::fixed; + deallog << std::setprecision(2); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + LinearAlgebra::Vector d1(N), d2(N); + LinearAlgebra::Vector f1(N), f2(N); + LinearAlgebra::Vector l1(N), l2(N); + + // cross-tests with double/float + // vectors at the same time are not + // supported at present, + // unfortunately, as many functions + // don't accept other data types as + // arguments + check_vectors (d1, d2); + check_vectors (f1, f2); + check_vectors (l1, l2); +} + + diff --git a/tests/lac/la_vector_vector.output b/tests/lac/la_vector_vector.output new file mode 100644 index 0000000000..46bd0872f1 --- /dev/null +++ b/tests/lac/la_vector_vector.output @@ -0,0 +1,76 @@ + +DEAL::Fill & Swap +DEAL::0 0 0 0 0 0 0 0 0 0 +DEAL::0 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 +DEAL::0 2.00 8.00 18.00 32.00 50.00 72.00 98.00 128.00 162.00 +DEAL::2.00 1.50 1.00 0.50 0 -0.50 -1.00 -1.50 -2.00 -2.50 +DEAL::0 2.00 8.00 18.00 32.00 50.00 72.00 98.00 128.00 162.00 +DEAL::0 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 +DEAL::2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 +DEAL::0 1.10 2.20 3.30 +DEAL::Extract number +DEAL::-105.00 -105.00 +DEAL::33.76 33.76 +DEAL::12.50 12.50 +DEAL::2.50 2.50 +DEAL::add & sub +DEAL::2.50 4.50 6.50 8.50 10.50 12.50 14.50 16.50 18.50 20.50 +DEAL::-2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 +DEAL::6.50 7.50 8.50 9.50 10.50 11.50 12.50 13.50 14.50 15.50 +DEAL::2.50 3.25 4.00 4.75 5.50 6.25 7.00 7.75 8.50 9.25 +DEAL::sadd & scale +DEAL::-2.50 -1.75 -1.00 -0.25 0.50 1.25 2.00 2.75 3.50 4.25 +DEAL::-3.75 -1.88 0 1.88 3.75 5.62 7.50 9.38 11.25 13.12 +DEAL::10.00 13.00 16.00 19.00 22.00 25.00 28.00 31.00 34.00 37.00 +DEAL::equ +DEAL::2.50 3.25 4.00 4.75 5.50 6.25 7.00 7.75 8.50 9.25 +DEAL::Fill & Swap +DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 +DEAL::0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 +DEAL::0.00 2.00 8.00 18.00 32.00 50.00 72.00 98.00 128.00 162.00 +DEAL::2.00 1.50 1.00 0.50 0.00 -0.50 -1.00 -1.50 -2.00 -2.50 +DEAL::0.00 2.00 8.00 18.00 32.00 50.00 72.00 98.00 128.00 162.00 +DEAL::0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 +DEAL::2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 +DEAL::0.00 1.10 2.20 3.30 +DEAL::Extract number +DEAL::-105.00 -105.00 +DEAL::33.76 33.76 +DEAL::12.50 12.50 +DEAL::2.50 2.50 +DEAL::add & sub +DEAL::2.50 4.50 6.50 8.50 10.50 12.50 14.50 16.50 18.50 20.50 +DEAL::-2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 +DEAL::6.50 7.50 8.50 9.50 10.50 11.50 12.50 13.50 14.50 15.50 +DEAL::2.50 3.25 4.00 4.75 5.50 6.25 7.00 7.75 8.50 9.25 +DEAL::sadd & scale +DEAL::-2.50 -1.75 -1.00 -0.25 0.50 1.25 2.00 2.75 3.50 4.25 +DEAL::-3.75 -1.88 0.00 1.88 3.75 5.62 7.50 9.38 11.25 13.12 +DEAL::10.00 13.00 16.00 19.00 22.00 25.00 28.00 31.00 34.00 37.00 +DEAL::equ +DEAL::2.50 3.25 4.00 4.75 5.50 6.25 7.00 7.75 8.50 9.25 +DEAL::Fill & Swap +DEAL::0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 +DEAL::0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 +DEAL::0.00 2.00 8.00 18.00 32.00 50.00 72.00 98.00 128.00 162.00 +DEAL::2.00 1.50 1.00 0.50 0.00 -0.50 -1.00 -1.50 -2.00 -2.50 +DEAL::0.00 2.00 8.00 18.00 32.00 50.00 72.00 98.00 128.00 162.00 +DEAL::0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 +DEAL::2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 +DEAL::0.00 1.10 2.20 3.30 +DEAL::Extract number +DEAL::-105.00 -105.00 +DEAL::33.76 33.76 +DEAL::12.50 12.50 +DEAL::2.50 2.50 +DEAL::add & sub +DEAL::2.50 4.50 6.50 8.50 10.50 12.50 14.50 16.50 18.50 20.50 +DEAL::-2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 -2.50 +DEAL::6.50 7.50 8.50 9.50 10.50 11.50 12.50 13.50 14.50 15.50 +DEAL::2.50 3.25 4.00 4.75 5.50 6.25 7.00 7.75 8.50 9.25 +DEAL::sadd & scale +DEAL::-2.50 -1.75 -1.00 -0.25 0.50 1.25 2.00 2.75 3.50 4.25 +DEAL::-3.75 -1.88 0.00 1.88 3.75 5.62 7.50 9.38 11.25 13.12 +DEAL::10.00 13.00 16.00 19.00 22.00 25.00 28.00 31.00 34.00 37.00 +DEAL::equ +DEAL::2.50 3.25 4.00 4.75 5.50 6.25 7.00 7.75 8.50 9.25