--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/vector_space_vector.h>
+
+#include <cstdio>
+#include <iostream>
+#include <cstring>
+#include <vector>
+
+
+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 <typename Number>
+ class Vector : public ReadWriteVector<Number>, public VectorSpaceVector<Number>
+ {
+ 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<Number> &V);
+
+ /**
+ * Constructor. Set dimension to @p n and initialize all elements with
+ * zero.
+ *
+ * The constructor is made explicit to avoid accident like this:
+ * <tt>v=0;</tt>. Presumably, the user wants to set every element of the
+ * vector to zero, but instead, what happens is this call:
+ * <tt>v=Vector@<Number@>(0);</tt>, i.e. the vector is replaced by one of
+ * length zero.
+ */
+ explicit Vector(const typename ReadWriteVector<Number>::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 <typename InputIterator>
+ Vector(const InputIterator first, const InputIterator last);
+
+ /**
+ * Destructor, deallocates memory.
+ */
+ ~Vector();
+
+ /**
+ * Multiply the entire vector by a fixed factor.
+ */
+ VectorSpaceVector<Number> &operator*= (const Number factor);
+
+ /**
+ * Divide the entire vector by a fixed factor.
+ */
+ VectorSpaceVector<Number> &operator/= (const Number factor);
+
+ /**
+ * Add the vector @p V to the present one.
+ */
+ VectorSpaceVector<Number> &operator+= (const VectorSpaceVector<Number> &V);
+
+ /**
+ * Substract the vector @p V from the present one.
+ */
+ VectorSpaceVector<Number> &operator-= (const VectorSpaceVector<Number> &V);
+
+ /**
+ * Return the scalar product of two vectors.
+ */
+ Number operator* (const VectorSpaceVector<Number> &V);
+
+ /**
+ * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
+ */
+ void add(const Number a, const VectorSpaceVector<Number> &V);
+
+ /**
+ * Multiple addition of a multiple of a vector, i.e. <tt>*this += a*V+b*W</tt>.
+ */
+ void add(const Number a, const VectorSpaceVector<Number> &V,
+ const Number b, const VectorSpaceVector<Number> &W);
+
+ /**
+ * Scaling and simple vector addition, i.e. <tt>*this = s*(*this)+V</tt>.
+ */
+ void sadd(const Number s, const VectorSpaceVector<Number> &V);
+
+ /**
+ * Scaling and simple addition of a multiple of a vector, i.e. <tt>*this =
+ * s*(*this)+a*V<tt>.
+ */
+ void sadd(const Number s, const Number a,
+ const VectorSpaceVector<Number> &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<Number> &scaling_factors);
+
+ /**
+ * Assignement <tt>*this = a*V</tt>.
+ */
+ void equ(const Number a, const VectorSpaceVector<Number> &V);
+
+ /**
+ * Returns the l<sub>1</sub> norm of the vector (i.e., the sum of the
+ * absolute values of all entries).
+ */
+ typename VectorSpaceVector<Number>::real_type l1_norm();
+
+ /**
+ * Returns the l<sub>2</sub> norm of the vector (i.e., the square root of
+ * the sum of the square of all entries among all processors).
+ */
+ typename VectorSpaceVector<Number>::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<Number>::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<Number> &V,
+ const VectorSpaceVector<Number> &W);
+
+ /**
+ * Returns the global size of the vector, equal to the sum of the number
+ * of locally owned indices among all processors.
+ */
+ typename ReadWriteVector<Number>::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<Number>::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<Number>::real_type l2_norm_recursive(unsigned int i,
+ unsigned int j);
+ };
+
+ /*@}*/
+ /*----------------------- Inline functions ----------------------------------*/
+
+ template <typename Number>
+ inline
+ Vector<Number>::Vector() {}
+
+
+
+ template <typename Number>
+ inline
+ Vector<Number>::Vector(const Vector<Number> &V)
+ :
+ ReadWriteVector<Number>(V)
+ {}
+
+
+
+ template <typename Number>
+ inline
+ Vector<Number>::Vector(const typename ReadWriteVector<Number>::size_type n)
+ :
+ ReadWriteVector<Number>(n)
+ {}
+
+
+
+ template <typename Number>
+ template <typename InputIterator>
+ inline
+ Vector<Number>::Vector(const InputIterator first, const InputIterator last)
+ {
+ this->reinit(complete_index_set(std::distance (first, last)), true);
+ std::copy(first, last, this->begin());
+ }
+
+
+
+ template <typename Number>
+ inline
+ Vector<Number>::~Vector() {}
+
+
+
+ template <typename Number>
+ inline
+ typename ReadWriteVector<Number>::size_type Vector<Number>::size() const
+ {
+ return ReadWriteVector<Number>::size();
+ }
+
+
+
+ template <typename Number>
+ inline
+ const dealii::IndexSet &Vector<Number>::locally_owned_elements() const
+ {
+ return ReadWriteVector<Number>::get_stored_elements();
+ }
+
+
+
+ template <typename Number>
+ inline
+ void Vector<Number>::print(std::ostream &out,
+ const unsigned int precision,
+ const bool scientific,
+ const bool) const
+ {
+ ReadWriteVector<Number>::print(out, precision, scientific);
+ }
+
+
+
+ template <typename Number>
+ inline
+ std::size_t Vector<Number>::memory_consumption() const
+ {
+ return ReadWriteVector<Number>::memory_consumption();
+ }
+} // end of namespace LinearAlgebra
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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/lac/la_vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+ template <typename Number>
+ VectorSpaceVector<Number> &Vector<Number>::operator*= (const Number factor)
+ {
+ AssertIsFinite(factor);
+ for (unsigned int i=0; i<this->size(); ++i)
+ this->val[i] *= factor;
+
+ return *this;
+ }
+
+
+
+ template <typename Number>
+ VectorSpaceVector<Number> &Vector<Number>::operator/= (const Number factor)
+ {
+ AssertIsFinite(factor);
+ Assert(factor!=Number(0.), ExcZero());
+ this->operator *= (Number(1.)/factor);
+
+ return *this;
+ }
+
+
+
+ template <typename Number>
+ VectorSpaceVector<Number> &Vector<Number>::operator+= (const VectorSpaceVector<Number> &V)
+ {
+ // Downcast V. If fails, throws an exception.
+ const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+ Assert(down_V.size()==this->size(),
+ ExcMessage("Cannot add two vectors with different numbers of elements"));
+ for (unsigned int i=0; i<this->size(); ++i)
+ {
+ AssertIsFinite(down_V[i]);
+ this->val[i] += down_V[i];
+ }
+
+ return *this;
+ }
+
+
+
+ template <typename Number>
+ VectorSpaceVector<Number> &Vector<Number>::operator-= (const VectorSpaceVector<Number> &V)
+ {
+ // Downcast V. If fails, throws an exception.
+ const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+ Assert(down_V.size()==this->size(),
+ ExcMessage("Cannot subtract two vectors with different numbers of elements"));
+ for (unsigned int i=0; i<this->size(); ++i)
+ {
+ AssertIsFinite(down_V[i]);
+ this->val[i] -= down_V[i];
+ }
+
+ return *this;
+ }
+
+
+
+ template <typename Number>
+ Number Vector<Number>::operator* (const VectorSpaceVector<Number> &V)
+ {
+ // Downcast V. If fails, throws an exception.
+ const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(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; i<this->size(); ++i)
+ {
+ AssertIsFinite(down_V[i]);
+ value += this->val[i]*down_V[i];
+ }
+
+ return value;
+ }
+
+
+
+ template <typename Number>
+ void Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &V)
+ {
+ // Downcast V. If fails, throws an exception.
+ const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+ AssertIsFinite(a);
+ Assert(down_V.size()==this->size(),
+ ExcMessage("Cannot add two vectors with different numbers of elements"));
+ for (unsigned int i=0; i<this->size(); ++i)
+ {
+ AssertIsFinite(down_V[i]);
+ this->val[i] += a*down_V[i];
+ }
+ }
+
+
+
+ template <typename Number>
+ void Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &V,
+ const Number b, const VectorSpaceVector<Number> &W)
+ {
+ // Downcast V. If fails, throws an exception.
+ const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+ // Downcast W. If fails, throws an exception.
+ const Vector<Number> &down_W = dynamic_cast<const Vector<Number>&>(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; i<this->size(); ++i)
+ {
+ AssertIsFinite(down_V[i]);
+ AssertIsFinite(down_W[i]);
+ this->val[i] += a*down_V[i]+b*down_W[i];
+ }
+ }
+
+
+
+ template <typename Number>
+ void Vector<Number>::sadd(const Number s, const VectorSpaceVector<Number> &V)
+ {
+ *this *= s;
+ *this += V;
+ }
+
+
+
+ template <typename Number>
+ void Vector<Number>::sadd(const Number s, const Number a,
+ const VectorSpaceVector<Number> &V)
+ {
+ *this *= s;
+ // Downcast V. It fails, throws an exception.
+ const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+ Vector<Number> tmp(down_V);
+ tmp *= a;
+ *this += tmp;
+ }
+
+
+
+ template <typename Number>
+ void Vector<Number>::scale(const VectorSpaceVector<Number> &scaling_factors)
+ {
+ // Downcast scaling_factors. If fails, throws an exception.
+ const Vector<Number> &down_scaling_factors =
+ dynamic_cast<const Vector<Number>&>(scaling_factors);
+ Assert(down_scaling_factors.size()==this->size(),
+ ExcMessage("Cannot add two vectors with different numbers of elements"));
+ for (unsigned int i=0; i<this->size(); ++i)
+ {
+ AssertIsFinite(down_scaling_factors[i]);
+ this->val[i] *= down_scaling_factors[i];
+ }
+ }
+
+
+
+ template <typename Number>
+ void Vector<Number>::equ(const Number a, const VectorSpaceVector<Number> &V)
+ {
+ // Downcast V. If fails, throws an exception.
+ const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+ *this = down_V;
+ *this *= a;
+ }
+
+
+
+ template <typename Number>
+ typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm()
+ {
+ Assert (this->size(), ExcEmptyObject());
+
+ typename ReadWriteVector<Number>::real_type norm = 0.;
+ norm = l1_norm_recursive(0,this->size()-1);
+
+ return norm;
+ }
+
+
+
+
+ template <typename Number>
+ typename VectorSpaceVector<Number>::real_type Vector<Number>::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<Number>::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<real_type>::min())
+ return std::sqrt(norm_square);
+ else
+ {
+ real_type scale = 0.;
+ real_type sum = 1.;
+ for (unsigned int i=0; i<this->size(); ++i)
+ {
+ if (this->val[i] != Number())
+ {
+ const real_type abs_x =
+ numbers::NumberTraits<Number>::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 Number>
+ typename VectorSpaceVector<Number>::real_type Vector<Number>::linfty_norm()
+ {
+ typename ReadWriteVector<Number>::real_type norm = 0.;
+ for (unsigned int i=0; i<this->size(); ++i)
+ if (std::abs(this->val[i])>norm)
+ norm = std::abs(this->val[i]);
+
+ return norm;
+ }
+
+
+
+ template <typename Number>
+ Number Vector<Number>::add_and_dot(const Number a,
+ const VectorSpaceVector<Number> &V,
+ const VectorSpaceVector<Number> &W)
+ {
+ this->add(a,V);
+
+ return *this*W;
+ }
+
+
+
+ template <typename Number>
+ typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm_recursive(unsigned int i,
+ unsigned int j)
+ {
+ Assert(j>=i, ExcInternalError());
+ typename ReadWriteVector<Number>::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 Number>
+ typename VectorSpaceVector<Number>::real_type Vector<Number>::l2_norm_recursive(unsigned int i,
+ unsigned int j)
+ {
+ Assert(j>=i, ExcInternalError());
+
+ typename ReadWriteVector<Number>::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
#include <deal.II/base/mpi.h>
#include <deal.II/base/numbers.h>
+
DEAL_II_NAMESPACE_OPEN
+class IndexSet;
+
namespace LinearAlgebra
{
template <typename Number>
class ReadWriteVector;
- class IndexSet;
/*! @addtogroup Vectors
*@{
typedef types::global_dof_index size_type;
typedef typename numbers::NumberTraits<Number>::real_type real_type;
+
/**
* Multiply the entire vector by a fixed factor.
*/
virtual void sadd(const Number s, const VectorSpaceVector<Number> &V) = 0;
/**
- * Scaling and simple additiion of a multiple of a vector, i.e. <tt>*this =
+ * Scaling and simple addition of a multiple of a vector, i.e. <tt>*this =
* s*(*this)+a*V</tt>.
*/
virtual void sadd(const Number s, const Number a,
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
* 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.
*/
virtual std::size_t memory_consumption() const = 0;
};
-
/*@}*/
}
constraint_matrix.cc
full_matrix.cc
lapack_full_matrix.cc
+ la_vector.cc
matrix_lib.cc
matrix_out.cc
parallel_vector.cc
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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/lac/la_vector.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+#include "la_vector.inst"
+}
+
+DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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<SCALAR>;
+ }
+
+
+
+for (SCALAR : COMPLEX_SCALARS)
+ {
+ template class Vector<SCALAR>;
+ }
+
\}
}
+
+for (SCALAR : COMPLEX_SCALARS)
+{
+ namespace LinearAlgebra
+ \{
+ template class ReadWriteVector<SCALAR>;
+ \}
+}
+
+for (S1, S2 : COMPLEX_SCALARS)
+{
+ namespace LinearAlgebra
+ \{
+ template void ReadWriteVector<S1>::reinit<S2> (const ReadWriteVector<S2>&,
+ const bool);
+ \}
+}
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/logstream.h>
+#include <deal.II/lac/la_vector.h>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+
+
+
+
+template <typename number>
+void check_norms ()
+{
+ for (unsigned int test=0; test<20; ++test)
+ {
+ const unsigned int size = Testing::rand() % 100000;
+ LinearAlgebra::Vector<number> vec (size);
+ for (unsigned int i=0; i<size; ++i)
+ vec(i) = static_cast<number>(Testing::rand())/static_cast<number>(RAND_MAX);
+ const typename LinearAlgebra::ReadWriteVector<number>::real_type norm = vec.l2_norm();
+ for (unsigned int i=0; i<30; ++i)
+ AssertThrow (vec.l2_norm() == norm, ExcInternalError());
+
+ LinearAlgebra::Vector<number> 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<float>();
+ check_norms<double>();
+ check_norms<long double>();
+ check_norms<std::complex<double> >();
+ deallog << "OK" << std::endl;
+}
+
+
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/logstream.h>
+#include <deal.II/lac/la_vector.h>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+
+
+
+
+template <typename number>
+void check ()
+{
+ for (unsigned int test=0; test<5; ++test)
+ {
+ const unsigned int size = 17 + test*1101;
+ LinearAlgebra::Vector<number> v1 (size), v2(size), v3(size), check(size);
+ for (unsigned int i=0; i<size; ++i)
+ {
+ v1(i) = 0.1 + 0.005 * i;
+ v2(i) = -5.2 + 0.18 * i;
+ v3(i) = 3.14159 + 2.7183/(1.+i);
+ }
+ check = v1;
+ const number factor = 0.01432;
+
+ v1.add(factor, v2);
+ const number prod = v1 * v3;
+ const number prod_check = check.add_and_dot(factor, v2, v3);
+ if (test == 0 && types_are_equal<number,double>::value)
+ {
+ deallog << "Vector add reference: ";
+ for (unsigned int i=0; i<size; ++i)
+ deallog << v1[i] << " ";
+ deallog << std::endl;
+ deallog << "Vector check reference: ";
+ for (unsigned int i=0; i<size; ++i)
+ deallog << check[i] << " ";
+ deallog << std::endl;
+ }
+
+ deallog << "Add and dot should be " << prod/static_cast<number>(size)
+ << ", is " << prod_check/static_cast<number>(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<float>();
+ check<double>();
+ check<long double>();
+ check<std::complex<double> >();
+ deallog << "OK" << std::endl;
+}
+
+
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/logstream.h>
+#include <deal.II/lac/la_vector.h>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <limits>
+
+
+
+void check_large_numbers()
+{
+ LinearAlgebra::Vector<float> 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<v.size(); ++i)
+ v[i] = (float) 0.;
+ v(5) = 1e-30;
+ v(9) = 1e-32;
+ const double correct2 = std::sqrt(1e-60 + 1e-64);
+ AssertThrow (std::abs(v.l2_norm() - correct2) < 1e-6*correct2, ExcInternalError());
+
+ LinearAlgebra::Vector<double> 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<w.size(); ++i)
+ w[i] = (float) 0.;
+ w(1) = 1e-302;
+ w(2) = 1e-303;
+ w(3) = 2e-303;
+ w(4) = 3e-303;
+ w(5) = -1e-303;
+ const double correct4 = std::sqrt(100. + 1. + 4. + 9 + 1.) * 1e-303;
+ AssertThrow (std::abs(w.l2_norm() - correct4) <= 1e-13*correct4, ExcInternalError());
+
+ deallog << "OK" << 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_large_numbers();
+}
+
+
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// checks that the l1, l2, lp norm is computed correctly for deal.II vectors
+// (check the summation algorithm), including an accuracy test (should not
+// lose more than 1 decimal also for 200000 vector entries)
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/la_vector.h>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+
+
+
+
+template <typename number>
+void check_norms ()
+{
+ const number acc = 1e1*std::numeric_limits<number>::epsilon();
+ unsigned int skip = 73;
+ for (unsigned int size=1; size<200000; size+=skip)
+ {
+ // test correctness
+ if (size > 10000)
+ skip += 17;
+ LinearAlgebra::Vector<number> vec(size);
+ for (unsigned int i=0; i<size; ++i)
+ vec(i) = i+1;
+
+ const number l1_norm = vec.l1_norm();
+ AssertThrow (std::abs(l1_norm-0.5*size*(size+1)) < acc*0.5*size*(size+1),
+ ExcInternalError());
+
+ // test accuracy of summation
+ const long double value = 3.14159265358979323846;
+ for (unsigned int i=0; i<size; ++i)
+ vec[i] = (number)value;
+ const number l1_norma = vec.l1_norm();
+ std::cout<<l1_norma<<std::endl;
+ AssertThrow (std::abs(l1_norma-value*size) < acc*size*value,
+ ExcInternalError());
+ const number l2_norma = vec.l2_norm();
+ AssertThrow (std::abs(l2_norma-value*std::sqrt((number)size)) < acc*std::sqrt(size)*value,
+ ExcInternalError());
+ }
+}
+
+
+template <typename number>
+void check_complex_norms ()
+{
+ const number acc = 1e1*std::numeric_limits<number>::epsilon();
+ unsigned int skip = 73;
+ for (unsigned int size=1; size<100000; size+=skip)
+ {
+ // test correctness
+ if (size > 10000)
+ skip += 17;
+ LinearAlgebra::Vector<std::complex<number> > vec(size);
+ long double sum = 0.;
+ for (unsigned int i=0; i<size; ++i)
+ {
+ vec(i) = std::complex<number>(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<long double> value (3.14159265358979323846, 0.1);
+ for (unsigned int i=0; i<size; ++i)
+ vec[i] = std::complex<number>(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<float>();
+ check_norms<double>();
+ check_norms<long double>();
+ check_complex_norms<double>();
+ check_complex_norms<float>();
+ deallog << "OK" << std::endl;
+}
+
+
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/logstream.h>
+#include <deal.II/lac/la_vector.h>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+
+
+
+
+const unsigned int N=10;
+unsigned int check_point = 0;
+
+
+
+
+template <typename number>
+void print (const LinearAlgebra::Vector<number> &v)
+{
+ for (unsigned int i=0; i<v.size(); ++i)
+ deallog << v(i) << '\t';
+ deallog << std::endl;
+}
+
+
+
+template <typename number1, typename number2>
+void check_vectors (LinearAlgebra::Vector<number1> &d1, LinearAlgebra::Vector<number2> &d2)
+{
+ deallog << "Fill & Swap" << std::endl;
+ LinearAlgebra::Vector<number1> d3(d1.size());
+ print (d3);
+
+ for (unsigned int i=0; i<N; ++i)
+ {
+ d1(i) = 2. * i;
+ d2(i) = .5 * d1(i)*d1(i);
+ d3(i) = 2. - .5 * i;
+ };
+
+ print (d1);
+ print (d2);
+ print (d3);
+
+ swap (d1, d2);
+ print (d1);
+
+ d1 = d2;
+ print (d1);
+
+ for (unsigned int i=0; i<N; ++i)
+ d1[i] = 2.5;
+ print (d1);
+
+ // initialize with iterators
+ number1 array[] = { 0.0, 1.1, 2.2, 3.3 };
+ LinearAlgebra::Vector<number1> 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; i<N; ++i)
+ sum += 4.*i-i*i;
+ deallog << d3 *d2 << '\t' << sum << std::endl;
+
+ sum = 0.;
+ for (unsigned int i=0; i<N; ++i)
+ sum += 4.*i*i;
+ sum = std::sqrt(sum);
+ deallog << d2.l2_norm() << '\t' << sum << std::endl;
+
+ sum = 0.;
+ for (unsigned int i=0; i<N; ++i)
+ sum += std::fabs(2.-.5*i);
+ deallog << d3.l1_norm() << '\t' << sum << std::endl;
+
+ sum = 0.;
+ for (unsigned int i=0; i<N; ++i)
+ {
+ double t = std::fabs(2.-.5*i);
+ if (t>sum) 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<double> d1(N), d2(N);
+ LinearAlgebra::Vector<float> f1(N), f2(N);
+ LinearAlgebra::Vector<long double> 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);
+}
+
+
--- /dev/null
+
+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