+++ /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 <boost/serialization/array.hpp>
-#include <boost/serialization/split_member.hpp>
-
-#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.
- *
- * @author Bruno Turcksin, 2015.
- */
- template <typename Number>
- class Vector : public ReadWriteVector<Number>, public VectorSpaceVector<Number>
- {
- public:
- typedef types::global_dof_index size_type;
-
- /**
- * 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 size_type n);
-
- /**
- * Initialize the vector with a given range of values pointed to by the
- * iterators. This function exists in analogy to the @p std::vector class.
- */
- template <typename InputIterator>
- Vector(const InputIterator first, const InputIterator last);
-
- /**
- * Destructor, deallocates memory.
- */
- virtual ~Vector();
-
- /**
- * Copies the data of the input vector @p in_vector.
- */
- Vector<Number> &operator= (const Vector<Number> &in_vector);
-
- /**
- * Copies the data of the input vector @p in_vector.
- */
- template <typename Number2>
- Vector<Number> &operator= (const Vector<Number2> &in_vector);
-
- /**
- * Sets all elements of the vector to the scalar @p s. This operation is
- * only allowed if @p s is equal to zero.
- */
- Vector<Number> &operator= (const Number s);
-
- /**
- * Multiply the entire vector by a fixed factor.
- */
- virtual VectorSpaceVector<Number> &operator*= (const Number factor);
-
- /**
- * Divide the entire vector by a fixed factor.
- */
- virtual VectorSpaceVector<Number> &operator/= (const Number factor);
-
- /**
- * Add the vector @p V to the present one.
- */
- virtual VectorSpaceVector<Number> &operator+= (const VectorSpaceVector<Number> &V);
-
- /**
- * Substract the vector @p V from the present one.
- */
- virtual VectorSpaceVector<Number> &operator-= (const VectorSpaceVector<Number> &V);
-
- /**
- * Return the scalar product of two vectors.
- */
- virtual Number operator* (const VectorSpaceVector<Number> &V);
-
- /**
- * Add @p a to all components. Note that @p a is a scalar not a vector.
- */
- virtual void add(const Number a);
-
- /**
- * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
- */
- virtual 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>.
- */
- virtual void add(const Number a, const VectorSpaceVector<Number> &V,
- const Number b, const VectorSpaceVector<Number> &W);
-
- /**
- * 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,
- 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.
- */
- virtual void scale(const VectorSpaceVector<Number> &scaling_factors);
-
- /**
- * Assignement <tt>*this = a*V</tt>.
- */
- virtual void equ(const Number a, const VectorSpaceVector<Number> &V);
-
- /**
- * Return the l<sub>1</sub> norm of the vector (i.e., the sum of the
- * absolute values of all entries).
- */
- virtual typename VectorSpaceVector<Number>::real_type l1_norm();
-
- /**
- * Return 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).
- */
- virtual typename VectorSpaceVector<Number>::real_type l2_norm();
-
- /**
- * Return the maximum norm of the vector (i.e., the maximum absolute value
- * among all entries and among all processors).
- */
- virtual typename VectorSpaceVector<Number>::real_type linfty_norm();
-
- /**
- * 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
- * this->add(a, V);
- * return_value = *this * W;
- * @endcode
- */
- virtual Number add_and_dot(const Number a,
- const VectorSpaceVector<Number> &V,
- const VectorSpaceVector<Number> &W);
-
- /**
- * Return the global size of the vector, equal to the sum of the number of
- * locally owned indices among all processors.
- */
- size_type size() const;
-
- /**
- * 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 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.
- */
- virtual void print(std::ostream &out,
- const unsigned int precision=3,
- const bool scientific=true,
- const bool across=true) const;
- /**
- * Write the vector en bloc to a file. This is done in a binary mode, so
- * the output is neither readable by humans nor (probably) by other
- * computers using a different operating system or number format.
- */
- void block_write (std::ostream &out) const;
-
- /**
- * Read a vector en block from a file. This is done using the inverse
- * operations to the above function, so it is reasonably fast because the
- * bitstream is not interpreted.
- *
- * The vector is resized if necessary.
- *
- * A primitive form of error checking is performed which will recognize
- * the bluntest attempts to interpret some data as a vector stored bitwise
- * to a file, but not more.
- */
- void block_read (std::istream &in);
-
- /**
- * 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 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 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<Number>::real_type l2_norm_squared_recursive(
- unsigned int i,
- unsigned int j);
-
- /**
- * Serialize the data of this object using boost. This function is
- * necessary to use boost::archive::text_iarchive and
- * boost::archive::text_oarchive.
- */
- template <typename Archive>
- void serialize(Archive &ar, const unsigned int version);
-
- friend class boost::serialization::access;
- };
-
- /*@}*/
- /*----------------------- 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 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
- Vector<Number> &Vector<Number>::operator= (const Vector<Number> &in_vector)
- {
- this->reinit(in_vector.size(),true);
- std::copy(in_vector.begin(), in_vector.end(), this->begin());
-
- return *this;
- }
-
-
-
- template <typename Number>
- template <typename Number2>
- inline
- Vector<Number> &Vector<Number>::operator= (const Vector<Number2> &in_vector)
- {
- this->reinit(in_vector.size(in_vector.size()),true);
- std::copy(in_vector.begin(), in_vector.end(), this->begin());
-
- return *this;
- }
-
-
-
- template <typename Number>
- inline
- Vector<Number> &Vector<Number>::operator= (const Number s)
- {
- Assert(s==static_cast<Number>(0), ExcMessage("Only 0 can be assigned to a vector."));
- (void) s;
-
- std::fill(this->begin(),this->end(),Number());
-
- return *this;
- }
-
-
-
- template <typename Number>
- inline
- void Vector<Number>::add(const Number a)
- {
- AssertIsFinite(a);
- for (unsigned int i=0; i<this->size(); ++i)
- this->val[i] += a;
- }
-
-
-
- template <typename Number>
- inline
- typename Vector<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>
- template <typename Archive>
- inline
- void Vector<Number>::serialize(Archive &ar, const unsigned int)
- {
- unsigned int current_size = this->size();
- ar &static_cast<Subscriptor &>(*this);
- ar &this->stored_elements;
- // If necessary, resize the vector during a read operation
- if (this->size()!=current_size)
- this->reinit(this->size());
- ar &boost::serialization::make_array(this->val,this->size());
- }
-
-
-
- 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>
-#include <iostream>
-#include <iomanip>
-
-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)
- {
- // Check that casting will work.
- Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL, ExcVectorTypeNotCompatible());
-
- // 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)
- {
- // Check that casting will work.
- Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL, ExcVectorTypeNotCompatible());
-
- // 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)
- {
- // Check that casting will work.
- Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
- ExcVectorTypeNotCompatible());
-
- // 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)
- {
- // Check that casting will work.
- Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
- ExcVectorTypeNotCompatible());
-
- // 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)
- {
- // Check that casting will work.
- Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
- ExcVectorTypeNotCompatible());
- // Check that casting will work.
- Assert(dynamic_cast<const Vector<Number>*>(&W)!=NULL,
- ExcVectorTypeNotCompatible());
-
- // 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 Number a,
- const VectorSpaceVector<Number> &V)
- {
- // Check that casting will work.
- Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
- ExcVectorTypeNotCompatible());
-
- *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)
- {
- // Check that casting will work.
- Assert(dynamic_cast<const Vector<Number>*>(&scaling_factors)!=NULL,
- ExcVectorTypeNotCompatible());
-
- // 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)
- {
- // Check that casting will work.
- Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
- ExcVectorTypeNotCompatible());
-
- // 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_squared_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)
- norm = std::max(std::abs(this->val[i]),norm);
-
- 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_squared_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_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);
-
- return norm;
- }
-
-
-
- template <typename Number>
- void Vector<Number>::block_write(std::ostream &out) const
- {
- AssertThrow(out, ExcIO());
-
- // Other version of the following
- // out << size() << std::endl << '[';
- // Reason: operator<< seems to use some resources that lead to problems in
- // a multithreaded environment.
- const size_type sz = this->size();
- char buf[16];
-#ifdef DEAL_II_WITH_64BIT_INDICES
- std::sprintf(buf, "%llu", sz);
-#else
- std::sprintf(buf, "%u", sz);
-#endif
- std::strcat(buf, "\n[");
-
- out.write(buf, std::strlen(buf));
- out.write(reinterpret_cast<const char *>(this->begin()),
- reinterpret_cast<const char *>(this->end())
- - reinterpret_cast<const char *>(this->begin()));
-
- // out << ']';
- const char outro = ']';
- out.write(&outro, 1);
-
- AssertThrow(out, ExcIO());
- }
-
-
-
- template <typename Number>
- void Vector<Number>::block_read(std::istream &in)
- {
- AssertThrow(in, ExcIO());
-
- size_type sz;
-
- char buf[16];
-
- in.getline(buf,16,'\n');
- sz = std::atoi(buf);
- this->reinit(sz,true);
-
- char c;
- // in >> c;
- in.read(&c,1);
- AssertThrow(c=='[', ExcIO());
-
- in.read(reinterpret_cast<char *>(this->begin()),
- reinterpret_cast<const char *>(this->end())
- - reinterpret_cast<const char *>(this->begin()));
-
- // in >> c;
- in.read(&c,1);
- AssertThrow(c=']', ExcIO());
- }
-}
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2015 - 2016 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__read_write_vector_h
-#define dealii__read_write_vector_h
-
-#include <deal.II/base/config.h>
-#include <deal.II/base/index_set.h>
-#include <deal.II/base/mpi.h>
-#include <deal.II/base/template_constraints.h>
-#include <deal.II/base/types.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/base/memory_consumption.h>
-#include <deal.II/base/thread_management.h>
-#include <deal.II/lac/vector_view.h>
-
-#include <cstring>
-#include <iomanip>
-
-
-DEAL_II_NAMESPACE_OPEN
-
-#ifdef DEAL_II_WITH_PETSC
-namespace PETScWrappers
-{
- namespace MPI
- {
- class Vector;
- }
-}
-#endif
-
-#ifdef DEAL_II_WITH_TRILINOS
-namespace TrilinosWrappers
-{
- namespace MPI
- {
- class Vector;
- }
-}
-#endif
-
-namespace LinearAlgebra
-{
- /*! @addtogroup Vectors
- *@{
- */
-
- /**
- * ReadWriteVector is intended to represent vectors in ${\mathbb R}^N$ for
- * which it stores all or a subset of elements. The latter case in important
- * in parallel computations, where $N$ may be so large that no processor can
- * actually all elements of a solution vector, but where this is also not
- * necessary: one typically only has to store the values of degrees of
- * freedom that live on cells that are locally owned plus potentially those
- * degrees of freedom that live on ghost cells.
- *
- * This class allows to access individual elements to be read or written.
- * However, it does not allow global operations such as taking the norm.
- * ReadWriteVector can be used to read and write elements in vectors derived
- * from VectorSpaceVector such as TrilinosWrappers::MPI::Vector and
- * PETScWrappers::MPI::Vector.
- *
- * <h3>Storing elements</h3> Most of the time, one will simply read from or
- * write into a vector of the current class using the global numbers of
- * these degrees of freedom. This is done using operator() or operator[]
- * which call global_to_local() to transform the <i>global</i> index into a
- * <i>local</i> one. In such cases, it is clear that one can only access
- * elements of the vector that the current object indeed stores.
- *
- * However, it is also possible to access elements in the order in which
- * they are stored by the current object. In other words, one is not
- * interested in accessing elements with their <i>global</i> indices, but
- * instead using an enumeration that only takes into account the elements
- * that are actually stored. This is facilitated by the local_element()
- * function. To this end, it is necessary to know <i>in which order</i> the
- * current class stores its element. The elements of all the consecutive
- * ranges are stored in ascending order of the first index of each range.
- * The function largest_range_starting_index() of IndexSet can be used to
- * get the first index of the largest range.
- *
- * @author Bruno Turcksin, 2015.
- */
- template <typename Number>
- class ReadWriteVector : public Subscriptor
- {
- public:
- /**
- * Declare standard types used in all containers. These types parallel
- * those in the <tt>C++</tt> standard libraries <tt>vector<...></tt>
- * class.
- */
- typedef Number value_type;
- typedef value_type *pointer;
- typedef const value_type *const_pointer;
- typedef value_type *iterator;
- typedef const value_type *const_iterator;
- typedef value_type &reference;
- typedef const value_type &const_reference;
- typedef types::global_dof_index size_type;
- typedef typename numbers::NumberTraits<Number>::real_type real_type;
-
- /**
- * @name 1: Basic Object-handling
- */
- //@{
- /**
- * Empty constructor.
- */
- ReadWriteVector ();
-
- /**
- * Copy constructor.
- */
- ReadWriteVector (const ReadWriteVector<Number> &in_vector);
-
- /**
- * Constructs a vector given the size, the stored elements have their
- * index in [0,size).
- */
- explicit ReadWriteVector (const size_type size);
-
- /**
- * Constructs a vector whose stored elements indices are given by the
- * IndexSet @p locally_stored_indices.
- */
- explicit ReadWriteVector (const IndexSet &locally_stored_indices);
-
- /**
- * Destructor.
- */
- ~ReadWriteVector ();
-
- /**
- * Sets 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).
- */
- 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 <typename Number2>
- void reinit(const ReadWriteVector<Number2> &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.
- */
- void reinit (const IndexSet &locally_stored_indices,
- const bool omit_zeroing_entries = false);
-
- /**
- * Swap the contents of this vector and the other vector @p v. One could
- * do this operation with a temporary variable and copying over the data
- * elements, but this function is significantly more efficient since it
- * only swaps the pointers to the data of the two vectors and therefore
- * does not need to allocate temporary storage and move data around.
- *
- * This function is analog to the the @p swap function of all C++ standard
- * containers. Also, there is a global function <tt>swap(u,v)</tt> that
- * simply calls <tt>u.swap(v)</tt>, again in analogy to standard
- * functions.
- */
- void swap (ReadWriteVector<Number> &v);
-
- /**
- * Copies the data and the IndexSet of the input vector @p in_vector.
- */
- ReadWriteVector<Number> &
- operator= (const ReadWriteVector<Number> &in_vector);
-
- /**
- * Copies the data and the IndexSet of the input vector @p in_vector.
- */
- template <typename Number2>
- ReadWriteVector<Number> &
- operator= (const ReadWriteVector<Number2> &in_vector);
-
-#ifdef DEAL_II_WITH_PETSC
- /**
- * Imports all the elements present in the vector's IndexSet from the
- * input vector @p petsc_vec.
- */
- ReadWriteVector<Number> &
- operator= (const PETScWrappers::MPI::Vector &petsc_vec);
-#endif
-
-#ifdef DEAL_II_WITH_TRILINOS
- /**
- * Imports all the elements present in the vector's IndexSet from the
- * input vector @p trilinos_vec.
- */
- ReadWriteVector<Number> &
- operator= (const TrilinosWrappers::MPI::Vector &trilinos_vec);
-#endif
-
- /**
- * Sets all elements of the vector to the scalar @p s. This operation is
- * only allowed if @p s is equal to zero.
- */
- ReadWriteVector<Number> &operator = (const Number s);
-
- /**
- * The value returned by this function denotes the dimension of the vector
- * spaces that are modeled by objects of this kind. However, objects of
- * the current class do not actually stores all elements of vectors of
- * this space but may, in fact store only a subset. The number of elements
- * stored is returned by n_elements() and is smaller or equal to the
- * number returned by the current function.
- */
- size_type size() const;
-
- /**
- * This function returns the number of elements stored. It is smaller or
- * equal to the dimension of the vector space that is modeled by an object
- * of this kind. This dimension is return by size().
- */
- size_type n_elements() const;
-
- /**
- * Return the IndexSet that represents the indices of the elements stored.
- */
- const IndexSet &get_stored_elements () const;
-
- /**
- * Make the @p ReadWriteVector class a bit like the <tt>vector<></tt>
- * class of the C++ standard library by returning iterators to the start
- * and end of the <i>locally stored</i> elements of this vector.
- */
- iterator begin ();
-
- /**
- * Returns constant iterator to the start of the locally stored elements
- * of the vector.
- */
- const_iterator begin () const;
-
- /**
- * Returns an iterator pointing to the element past the end of the array
- * of locally stored entries.
- */
- iterator end ();
-
- /**
- * Returns a constant iterator pointing to the element past the end of the
- * array of the locally stored entries.
- */
- const_iterator end () const;
- //@}
-
-
- /**
- * @name 2: Data-Access
- */
- //@{
-
- /**
- * Read access to the data in the position corresponding to @p
- * global_index. An exception is thrown if @p global_index is not stored
- * by the current object.
- */
- Number operator () (const size_type global_index) const;
-
- /**
- * Read and write access to the data in the position corresponding to @p
- * global_index. An exception is thrown if @p global_index is not stored
- * by the current object.
- */
- Number &operator () (const size_type global_index);
-
- /**
- * Read access to the data in the position corresponding to @p
- * global_index. An exception is thrown if @p global_index is not stored
- * by the current object.
- *
- * This function does the same thing as operator().
- */
- Number operator [] (const size_type global_index) const;
-
- /**
- * Read and write access to the data in the position corresponding to @p
- * global_index. An exception is thrown if @p global_index is not stored
- * by the current object.
- *
- * This function does the same thing as operator().
- */
- Number &operator [] (const size_type global_index);
-
- /**
- * Instead of getting individual elements of a vector, this function
- * allows to get a whole set of elements at once. The indices of the
- * elements to be read are stated in the first argument, the corresponding
- * values are returned in the second.
- */
- template <typename Number2>
- void extract_subvector_to (const std::vector<size_type> &indices,
- std::vector<Number2> &values) const;
-
- /**
- * Just as the above, but with pointers. Useful in minimizing copying of
- * data around.
- */
- template <typename ForwardIterator, typename OutputIterator>
- void extract_subvector_to (ForwardIterator indices_begin,
- const ForwardIterator indices_end,
- OutputIterator values_begin) const;
-
- /**
- * Read access to the data field specified by @p local_index. When you
- * access elements in the order in which they are stored, it is necessary
- * that you know in which they are stored. In other words, you need to
- * know the map between the global indices of the elements this class
- * stores, and the local indices into the contiguous array of these global
- * elements. For this, see the general documentation of this class.
- *
- * Performance: Direct array access (fast).
- */
- Number local_element (const size_type local_index) const;
-
- /**
- * Read and write access to the data field specified by @p local_index.
- * When you access elements in the order in which they are stored, it is
- * necessary that you know in which they are stored. In other words, you
- * need to know the map between the global indices of the elements this
- * class stores, and the local indices into the contiguous array of these
- * global elements. For this, see the general documentation of this class.
- *
- * Performance: Direct array access (fast).
- */
- Number &local_element (const size_type local_index);
- //@}
-
-
- /**
- * @name 3: Modification of vectors
- */
- //@{
-
- /**
- * This function adds a whole set of values stored in @p values to the
- * vector components specified by @p indices.
- */
- template <typename Number2>
- void add (const std::vector<size_type> &indices,
- const std::vector<Number2> &values);
-
- /**
- * This function is similar to the previous one but takes a
- * ReadWriteVector of values.
- */
- template <typename Number2>
- void add (const std::vector<size_type> &indices,
- const ReadWriteVector<Number2> &values);
-
- /**
- * Take an address where <tt>n_elements</tt> are stored contiguously and
- * add them into the vector. Handles all cases which are not covered by
- * the other two <tt>add()</tt> functions above.
- */
- template <typename Number2>
- void add (const size_type n_elements,
- const size_type *indices,
- const Number2 *values);
-
- /**
- * Prints the vector to the output stream @p out.
- */
- void print (std::ostream &out,
- const unsigned int precision = 3,
- const bool scientific = true) const;
-
- /**
- * Returns the memory consumption of this class in bytes.
- */
- std::size_t memory_consumption () const;
- //@}
-
- protected:
-
- /**
- * Return the local position of @p global_index.
- */
- unsigned int
- global_to_local (const types::global_dof_index global_index) const
- {
- // the following will throw an exception if the global_index is not
- // in the remaining_elements
- return static_cast<unsigned int>(stored_elements.index_within_set(global_index));
- }
-
- /**
- * A helper function that is used to resize the val array.
- */
- void resize_val (const size_type new_allocated_size);
-
- /**
- * Indices of the elements stored.
- */
- IndexSet stored_elements;
-
- /**
- * Pointer to the array of local elements of this vector.
- */
- // TODO: use AlignedVector here for storage
- Number *val;
- };
-
- /*@}*/
-
-
- /*----------------------- Inline functions ----------------------------------*/
-
-#ifndef DOXYGEN
-
- template <typename Number>
- inline
- ReadWriteVector<Number>::ReadWriteVector ()
- :
- val (NULL)
- {}
-
-
-
- template <typename Number>
- inline
- ReadWriteVector<Number>::ReadWriteVector (const ReadWriteVector<Number> &v)
- :
- Subscriptor(),
- val (NULL)
- {
- this->operator=(v);
- }
-
-
-
- template <typename Number>
- inline
- ReadWriteVector<Number>::ReadWriteVector (const size_type size)
- :
- val (NULL)
- {
- reinit (size, false);
- }
-
-
-
- template <typename Number>
- inline
- ReadWriteVector<Number>::ReadWriteVector (const IndexSet &locally_stored_indices)
- :
- val (NULL)
- {
- reinit (locally_stored_indices);
- }
-
-
-
- template <typename Number>
- inline
- ReadWriteVector<Number>::~ReadWriteVector ()
- {
- resize_val(0);
- }
-
-
-
- template <typename Number>
- inline
- ReadWriteVector<Number> &
- ReadWriteVector<Number>::operator= (const ReadWriteVector<Number> &in_vector)
- {
- resize_val(in_vector.n_elements());
- stored_elements = in_vector.get_stored_elements();
- std::copy(in_vector.begin(),in_vector.end(),begin());
-
- return *this;
- }
-
-
-
- template <typename Number>
- template <typename Number2>
- inline
- ReadWriteVector<Number> &
- ReadWriteVector<Number>::operator= (const ReadWriteVector<Number2> &in_vector)
- {
- resize_val(in_vector.n_elements());
- stored_elements = in_vector.get_stored_elements();
- std::copy(in_vector.begin(),in_vector.end(),begin());
-
- return *this;
- }
-
-
-
- template <typename Number>
- inline
- typename ReadWriteVector<Number>::size_type
- ReadWriteVector<Number>::size() const
- {
- return stored_elements.size();
- }
-
-
-
- template <typename Number>
- inline
- typename ReadWriteVector<Number>::size_type
- ReadWriteVector<Number>::n_elements() const
- {
- return stored_elements.n_elements();
- }
-
-
-
- template <typename Number>
- inline
- const IndexSet &
- ReadWriteVector<Number>::get_stored_elements () const
- {
- return stored_elements;
- }
-
-
-
- template <typename Number>
- inline
- typename ReadWriteVector<Number>::iterator
- ReadWriteVector<Number>::begin ()
- {
- return &val[0];
- }
-
-
-
- template <typename Number>
- inline
- typename ReadWriteVector<Number>::const_iterator
- ReadWriteVector<Number>::begin () const
- {
- return &val[0];
- }
-
-
-
- template <typename Number>
- inline
- typename ReadWriteVector<Number>::iterator
- ReadWriteVector<Number>::end ()
- {
- return &val[this->n_elements()];
- }
-
-
-
- template <typename Number>
- inline
- typename ReadWriteVector<Number>::const_iterator
- ReadWriteVector<Number>::end () const
- {
- return &val[this->n_elements()];
- }
-
-
-
- template <typename Number>
- inline
- Number
- ReadWriteVector<Number>::operator() (const size_type global_index) const
- {
- return val[global_to_local(global_index)];
- }
-
-
-
- template <typename Number>
- inline
- Number &
- ReadWriteVector<Number>::operator() (const size_type global_index)
- {
- return val[global_to_local (global_index)];
- }
-
-
-
- template <typename Number>
- inline
- Number
- ReadWriteVector<Number>::operator[] (const size_type global_index) const
- {
- return operator()(global_index);
- }
-
-
-
- template <typename Number>
- inline
- Number &
- ReadWriteVector<Number>::operator[] (const size_type global_index)
- {
- return operator()(global_index);
- }
-
-
-
- template <typename Number>
- template <typename Number2>
- inline
- void ReadWriteVector<Number>::extract_subvector_to (const std::vector<size_type> &indices,
- std::vector<Number2> &values) const
- {
- for (size_type i = 0; i < indices.size(); ++i)
- values[i] = operator()(indices[i]);
- }
-
-
-
- template <typename Number>
- template <typename ForwardIterator, typename OutputIterator>
- inline
- void ReadWriteVector<Number>::extract_subvector_to (ForwardIterator indices_begin,
- const ForwardIterator indices_end,
- OutputIterator values_begin) const
- {
- while (indices_begin != indices_end)
- {
- *values_begin = operator()(*indices_begin);
- indices_begin++;
- values_begin++;
- }
- }
-
-
-
- template <typename Number>
- inline
- Number
- ReadWriteVector<Number>::local_element (const size_type local_index) const
- {
- AssertIndexRange (local_index, this->n_elements());
-
- return val[local_index];
- }
-
-
-
- template <typename Number>
- inline
- Number &
- ReadWriteVector<Number>::local_element (const size_type local_index)
- {
- AssertIndexRange (local_index, this->n_elements());
-
- return val[local_index];
- }
-
-
-
- template <typename Number>
- inline
- ReadWriteVector<Number> &
- ReadWriteVector<Number>::operator= (const Number s)
- {
- Assert(s==static_cast<Number>(0), ExcMessage("Only 0 can be assigned to a vector."));
- (void)s;
-
- std::fill(begin(),end(),Number());
-
- return *this;
- }
-
-
-
- template <typename Number>
- template <typename Number2>
- inline
- void
- ReadWriteVector<Number>::add (const std::vector<size_type> &indices,
- const std::vector<Number2> &values)
- {
- AssertDimension (indices.size(), values.size());
- add (indices.size(), &indices[0], &values[0]);
- }
-
-
-
- template <typename Number>
- template <typename Number2>
- inline
- void
- ReadWriteVector<Number>::add (const std::vector<size_type> &indices,
- const ReadWriteVector<Number2> &values)
- {
- for (size_type i=0; i<indices.size(); ++i)
- {
- Assert (numbers::is_finite(values[i]),
- ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
- this->operator()(indices[i]) += values[indices[i]];
- }
- }
-
-
-
- template <typename Number>
- template <typename Number2>
- inline
- void
- ReadWriteVector<Number>::add (const size_type n_indices,
- const size_type *indices,
- const Number2 *values)
- {
- for (size_type i=0; i<n_indices; ++i)
- {
- Assert (numbers::is_finite(values[i]),
- ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
- this->operator()(indices[i]) += values[i];
- }
- }
-
-#endif // ifndef DOXYGEN
-
-} // end of namespace LinearAlgebra
-
-
-
-
-/**
- * Global function @p swap which overloads the default implementation of the
- * C++ standard library which uses a temporary object. The function simply
- * exchanges the data of the two vectors.
- *
- * @relates Vector
- */
-template <typename Number>
-inline
-void swap (LinearAlgebra::ReadWriteVector<Number> &u,
- LinearAlgebra::ReadWriteVector<Number> &v)
-{
- u.swap (v);
-}
-
-
-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__parallel_vector_templates_h
-#define dealii__parallel_vector_templates_h
-
-
-#include <deal.II/base/config.h>
-#include <deal.II/lac/read_write_vector.h>
-#include <deal.II/lac/petsc_parallel_vector.h>
-#include <deal.II/lac/trilinos_vector.h>
-
-#ifdef DEAL_II_WITH_TRILINOS
-# include <Epetra_Import.h>
-#endif
-
-DEAL_II_NAMESPACE_OPEN
-
-
-namespace LinearAlgebra
-{
- template <typename Number>
- void
- ReadWriteVector<Number>::resize_val (const size_type new_alloc_size)
- {
- if (new_alloc_size == 0)
- {
- if (val != NULL)
- free(val);
- val = NULL;
- }
- else
- {
- if (val != NULL)
- free(val);
-
- Utilities::System::posix_memalign ((void **)&val, 64, sizeof(Number)*new_alloc_size);
- }
- }
-
-
-
- template <typename Number>
- void
- ReadWriteVector<Number>::reinit (const size_type size,
- const bool omit_zeroing_entries)
- {
- // check whether we need to reallocate
- resize_val(size);
-
- stored_elements = complete_index_set(size);
- stored_elements.compress();
-
- // set entries to zero if so requested
- if (omit_zeroing_entries == false)
- this->operator = (Number());
- }
-
-
-
- template <typename Number>
- template <typename Number2>
- void
- ReadWriteVector<Number>::reinit (const ReadWriteVector<Number2> &v,
- const bool omit_zeroing_entries)
- {
- resize_val(v.n_elements());
-
- stored_elements = v.get_stored_elements();
-
- if (omit_zeroing_entries == false)
- this->operator= (Number());
- }
-
-
-
- template <typename Number>
- void
- ReadWriteVector<Number>::reinit (const IndexSet &locally_stored_indices,
- const bool omit_zeroing_entries)
- {
- stored_elements = locally_stored_indices;
-
- // set vector size and allocate memory
- resize_val(stored_elements.n_elements());
-
- // initialize to zero
- if (omit_zeroing_entries == false)
- this->operator= (Number());
- }
-
-
-
-#ifdef DEAL_II_WITH_PETSC
- namespace internal
- {
- template <typename PETSC_Number, typename Number>
- void copy_petsc_vector (const PETSC_Number *petsc_start_ptr,
- const PETSC_Number *petsc_end_ptr,
- Number *ptr)
- {
- std::copy(petsc_start_ptr, petsc_end_ptr, ptr);
- }
-
- template <typename PETSC_Number, typename Number>
- void copy_petsc_vector (const std::complex<PETSC_Number> *petsc_start_ptr,
- const std::complex<PETSC_Number> *petsc_end_ptr,
- std::complex<Number> *ptr)
- {
- std::copy(petsc_start_ptr, petsc_end_ptr, ptr);
- }
-
- template <typename PETSC_Number, typename Number>
- void copy_petsc_vector (const std::complex<PETSC_Number> *petsc_start_ptr,
- const std::complex<PETSC_Number> *petsc_end_ptr,
- Number *ptr)
- {
- AssertThrow(false, ExcMessage("Tried to copy complex -> real"));
- }
- }
-
-
-
- template <typename Number>
- ReadWriteVector<Number> &
- ReadWriteVector<Number>::operator = (const PETScWrappers::MPI::Vector &petsc_vec)
- {
- //TODO: this works only if no communication is needed.
- Assert(petsc_vec.locally_owned_elements() == stored_elements,
- StandardExceptions::ExcInvalidState());
-
- // get a representation of the vector and copy it
- PetscScalar *start_ptr;
- int ierr = VecGetArray (static_cast<const Vec &>(petsc_vec), &start_ptr);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
-
- const size_type vec_size = petsc_vec.local_size();
- internal::copy_petsc_vector (start_ptr, start_ptr + vec_size, begin());
-
- // restore the representation of the vector
- ierr = VecRestoreArray (static_cast<const Vec &>(petsc_vec), &start_ptr);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
-
- // return a pointer to this object per normal c++ operator overloading
- // semantics
- return *this;
- }
-#endif
-
-
-
-#ifdef DEAL_II_WITH_TRILINOS
- template <typename Number>
- ReadWriteVector<Number> &
- ReadWriteVector<Number>::operator = (const TrilinosWrappers::MPI::Vector &trilinos_vec)
- {
- // Copy the local elements
- std::vector<size_type> trilinos_indices, readwrite_indices;
- trilinos_vec.locally_owned_elements().fill_index_vector(trilinos_indices);
- stored_elements.fill_index_vector(readwrite_indices);
- std::vector<size_type> intersection(trilinos_indices.size());
- std::vector<size_type>::iterator end;
- end = std::set_intersection(trilinos_indices.begin(), trilinos_indices.end(),
- readwrite_indices.begin(), readwrite_indices.end(),
- intersection.begin());
- const size_t intersection_size = end-intersection.begin();
- intersection.resize(intersection_size);
- std::vector<size_t> trilinos_position(intersection_size);
- unsigned int j=0;
- for (size_t i=0; i<intersection_size; ++i)
- {
- while (intersection[i]!=trilinos_indices[j])
- ++j;
- trilinos_position[i] = j;
- }
- double *values = trilinos_vec.trilinos_vector().Values();
- for (size_t i=0; i<trilinos_position.size(); ++i)
- val[global_to_local(intersection[i])] = values[trilinos_position[i]];
-
-
-#ifdef DEAL_II_WITH_MPI
- // Copy the off-processor elements if necessary
- if (intersection_size != n_elements())
- {
- Epetra_BlockMap source_map = trilinos_vec.trilinos_vector().Map();
- std::vector<size_type> off_proc_indices(readwrite_indices.size());
- // Subtract the elements of readwrite that are in intersection.
- end = std::set_difference(readwrite_indices.begin(), readwrite_indices.end(),
- intersection.begin(), intersection.end(), off_proc_indices.begin());
- off_proc_indices.resize(end-off_proc_indices.begin());
- // Cast size_type to TrilinosWrappers::type::int_type
- TrilinosWrappers::types::int_type *off_proc_array =
- new TrilinosWrappers::types::int_type [off_proc_indices.size()];
- for (size_t i=0; i<off_proc_indices.size(); ++i)
- off_proc_array[i] = static_cast<TrilinosWrappers::types::int_type> (
- off_proc_indices[i]);
- Epetra_Map target_map(off_proc_indices.size(),off_proc_indices.size(),
- off_proc_array,0,source_map.Comm());
- delete [] off_proc_array;
- Epetra_Import import(target_map, source_map);
- Epetra_FEVector target_vector(target_map);
- target_vector.Import(trilinos_vec.trilinos_vector(), import, Insert);
- values = target_vector.Values();
- for (unsigned int i=0; i<off_proc_indices.size(); ++i)
- val[global_to_local(off_proc_indices[i])] = values[i];
- }
-#endif
-
- // return a pointer to this object per normal c++ operator overloading
- // semantics
- return *this;
- }
-#endif
-
-
-
- template <typename Number>
- void
- ReadWriteVector<Number>::swap (ReadWriteVector<Number> &v)
- {
- std::swap(stored_elements, v.stored_elements);
- std::swap(val, v.val);
- }
-
-
-
- template <typename Number>
- std::size_t
- ReadWriteVector<Number>::memory_consumption () const
- {
- std::size_t memory = sizeof(*this);
- memory += sizeof (Number) * static_cast<std::size_t>(this->n_elements());
-
- memory += stored_elements.memory_consumption();
-
- return memory;
- }
-
-
-
- template <typename Number>
- void
- ReadWriteVector<Number>::print (std::ostream &out,
- const unsigned int precision,
- const bool scientific) const
- {
- AssertThrow (out, ExcIO());
- std::ios::fmtflags old_flags = out.flags();
- unsigned int old_precision = out.precision (precision);
-
- out.precision (precision);
- if (scientific)
- out.setf (std::ios::scientific, std::ios::floatfield);
- else
- out.setf (std::ios::fixed, std::ios::floatfield);
-
- out << "IndexSet: ";
- stored_elements.print(out);
- out << std::endl;
- for (unsigned int i=0; i<this->n_elements(); ++i)
- out << val[i] << std::endl;
- out << std::flush;
-
-
- AssertThrow (out, ExcIO());
- // reset output format
- out.flags (old_flags);
- out.precision(old_precision);
- }
-} // end of namespace LinearAlgebra
-
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2015 - 2016 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__vector_space_vector_h
-#define dealii__vector_space_vector_h
-
-#include <deal.II/base/config.h>
-#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;
-
-
- /*! @addtogroup Vectors
- *@{
- */
-
- /**
- * VectorSpaceVector is an abstract class that is used to define the
- * interface that vector classes need to implement when they want to
- * implement global operations. This class is complementary of
- * ReadWriteVector which allows the access of individual elements but does
- * not allow global operations.
- *
- * @author Bruno Turcksin, 2015.
- */
- template <typename Number>
- class VectorSpaceVector
- {
- public:
- typedef Number value_type;
- 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 VectorSpaceVector<Number> &operator*= (const Number factor) = 0;
-
- /**
- * Divide the entire vector by a fixed factor.
- */
- virtual VectorSpaceVector<Number> &operator/= (const Number factor) = 0;
-
- /**
- * Add the vector @p V to the present one.
- */
- virtual VectorSpaceVector<Number> &operator+= (const VectorSpaceVector<Number> &V) = 0;
-
- /**
- * Subtract the vector @p V from the present one.
- */
- virtual VectorSpaceVector<Number> &operator-= (const VectorSpaceVector<Number> &V) = 0;
-
- /**
- * Return the scalar product of two vectors.
- */
- virtual Number operator* (const VectorSpaceVector<Number> &V) = 0;
-
- /**
- * Add @p a to all components. Note that @p a is a scalar not a vector.
- */
- virtual void add(const Number a) = 0;
-
- /**
- * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
- */
- virtual void add(const Number a, const VectorSpaceVector<Number> &V) = 0;
-
- /**
- * Multiple addition of scaled vectors, i.e. <tt>*this += a*V+b*W</tt>.
- */
- virtual void add(const Number a, const VectorSpaceVector<Number> &V,
- const Number b, const VectorSpaceVector<Number> &W) = 0;
-
- /**
- * 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,
- const VectorSpaceVector<Number> &V) = 0;
-
- /**
- * 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.
- */
- virtual void scale(const VectorSpaceVector<Number> &scaling_factors) = 0;
-
- /**
- * Assignement <tt>*this = a*V</tt>.
- */
- virtual void equ(const Number a, const VectorSpaceVector<Number> &V) = 0;
-
- /**
- * Return the l<sub>1</sub> norm of the vector (i.e., the sum of the
- * absolute values of all entries among all processors).
- */
- virtual real_type l1_norm() = 0;
-
- /**
- * Return 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).
- */
- virtual real_type l2_norm() = 0;
-
- /**
- * 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;
-
- /**
- * 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
- * this->add(a, V);
- * return_value = *this * W;
- * @endcode
- *
- * The reason this function exists is that this operation involves less
- * memory transfer than calling the two functions separately. This method
- * only needs to load three vectors, @p this, @p V, @p W, whereas calling
- * separate methods means to load the calling vector @p this twice. Since
- * most vector operations are memory transfer limited, this reduces the
- * time by 25\% (or 50\% if @p W equals @p this).
- */
- virtual Number add_and_dot(const Number a,
- const VectorSpaceVector<Number> &V,
- const VectorSpaceVector<Number> &W) = 0;
-
- /**
- * 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;
-
- /**
- * 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 created on only one processor, then the result would
- * satisfy
- * @code
- * vec.locally_owned_elements() == complete_index_set(vec.size())
- * @endcode
- */
- virtual const dealii::IndexSet &locally_owned_elements() const = 0;
-
- /**
- * Print the vector to the output stream @p out.
- */
- virtual void print(std::ostream &out,
- const unsigned int precision=3,
- const bool scientific=true,
- const bool across=true) const = 0;
-
- /**
- * Return the memory consumption of this class in bytes.
- */
- virtual std::size_t memory_consumption() const = 0;
- };
- /*@}*/
-}
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
constraint_matrix.cc
full_matrix.cc
lapack_full_matrix.cc
- la_vector.cc
matrix_lib.cc
matrix_out.cc
parallel_vector.cc
precondition_block.cc
precondition_block_ez.cc
relaxation_block.cc
- read_write_vector.cc
solver.cc
solver_control.cc
sparse_decomposition.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
- read_write_vector.inst.in
solver.inst.in
sparse_matrix_ez.inst.in
sparse_matrix.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 - 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.
-//
-// ---------------------------------------------------------------------
-
-
-
-for (SCALAR : REAL_SCALARS)
- {
- template class Vector<SCALAR>;
- }
-
-
-
-for (SCALAR : COMPLEX_SCALARS)
- {
- template class Vector<SCALAR>;
- }
-
+++ /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/read_write_vector.h>
-#include <deal.II/lac/read_write_vector.templates.h>
-
-DEAL_II_NAMESPACE_OPEN
-
-#include "read_write_vector.inst"
-
-DEAL_II_NAMESPACE_CLOSE
+++ /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.
-//
-// ---------------------------------------------------------------------
-
-
-
-for (SCALAR : REAL_SCALARS)
-{
- namespace LinearAlgebra
- \{
- template class ReadWriteVector<SCALAR>;
- \}
-}
-
-for (S1, S2 : REAL_SCALARS)
-{
- namespace LinearAlgebra
- \{
- template void ReadWriteVector<S1>::reinit<S2> (const ReadWriteVector<S2>&,
- const bool);
- \}
-}
-
-
-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.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 - 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 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);
- // Check that the assignment works
- v1 = 0.;
- 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;
-
- const number constant = 1.;
- v1.add(constant);
- deallog << "Vector add constant: ";
- for (unsigned int i=0; i<size; ++i)
- deallog << v1[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.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::Vector add constant: 1.03 1.03 1.04 1.05 1.06 1.06 1.07 1.08 1.09 1.09 1.10 1.11 1.12 1.12 1.13 1.14 1.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.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 = 1e2*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.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) 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 <fstream>
-#include <string>
-#include <boost/archive/text_oarchive.hpp>
-#include <boost/archive/text_iarchive.hpp>
-
-
-void test()
-{
- const char *filename = "test.txt";
- const unsigned int size(10);
- LinearAlgebra::Vector<double> vec(size);
- for (unsigned int i=0; i<size; ++i)
- vec[i] = i;
-
- // Check block_write and block_read
- std::ofstream file_out(filename);
- // Write the vector in the file
- vec.block_write(file_out);
- file_out.close();
- // Clear the vector
- vec.reinit(0);
- // Load the vector form the file
- std::ifstream file_in(filename);
- vec.block_read(file_in);
- file_in.close();
- // Check the vector
- double eps=1e-12;
- for (unsigned int i=0; i<size; ++i)
- AssertThrow(std::abs(vec[i]-(double)i)<eps,
- ExcMessage("Value in the vector has been changed by block_write or block_read"));
-
-
- // save data to archive
- {
- std::ofstream file_out2(filename);
- boost::archive::text_oarchive oa(file_out2);
- oa<<vec;
- // archive and stream closed when destructors are called
- }
-
- // Clear the vector
- vec.reinit(0);
- {
- std::ifstream file_in2(filename);
- boost::archive::text_iarchive ia(file_in2);
- ia >> vec;
- }
- // Check the vector
- for (unsigned int i=0; i<size; ++i)
- AssertThrow(std::abs(vec[i]-(double)i)<eps,
- ExcMessage("Value in the vector has been changed by boost archive"));
-}
-
-
-int main()
-{
- std::ofstream logfile("output");
- deallog << std::fixed;
- deallog << std::setprecision(2);
- deallog.attach(logfile);
- deallog.threshold_double(1.e-10);
-
- test();
-
- 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., 1., 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.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
+++ /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.
-//
-// ---------------------------------------------------------------------
-
-// Check reinit and add
-
-#include "../tests.h"
-#include <deal.II/base/index_set.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/read_write_vector.h>
-
-#include <fstream>
-#include <vector>
-
-void test()
-{
- unsigned int size(10);
- std::vector<types::global_dof_index> indices(size);
- std::vector<float> float_values(size);
- std::vector<double> double_values(size);
- for (unsigned int i=0; i<size; ++i)
- {
- indices[i] = i;
- float_values[i] = 2.*i;
- double_values[i] = 3.*i;
- }
-
- LinearAlgebra::ReadWriteVector<float> float_vector;
- LinearAlgebra::ReadWriteVector<double> double_vector_1;
- LinearAlgebra::ReadWriteVector<double> double_vector_2;
- IndexSet is(50);
- is.add_range(0,10);
- is.add_index(46);
- is.add_range(11,25);
- is.compress();
-
- float_vector.reinit(size);
- double_vector_1.reinit(float_vector);
- double_vector_2.reinit(is);
-
- IndexSet is_copy(double_vector_2.get_stored_elements());
- is.print(std::cout);
- is_copy.print(std::cout);
-
- float_vector.add(indices,float_values);
- double_vector_1.add(indices,float_vector);
- double_vector_2.add(size,&indices[0],&double_values[0]);
- double_vector_1.swap(double_vector_2);
- LinearAlgebra::ReadWriteVector<double>::iterator val_1 = double_vector_1.begin();
- LinearAlgebra::ReadWriteVector<double>::iterator end_1 = double_vector_1.end();
- LinearAlgebra::ReadWriteVector<double>::iterator val_2 = double_vector_2.begin();
- LinearAlgebra::ReadWriteVector<double>::iterator end_2 = double_vector_2.end();
- deallog << "double_vector_1" << std::endl;
- for (; val_1<end_1; ++val_1)
- deallog << *val_1 <<std::endl;
- deallog << "double_vector_2" << std::endl;
- for (; val_2<end_2; ++val_2)
- deallog << *val_2 <<std::endl;
- double_vector_1.extract_subvector_to(indices,float_values);
- deallog << "subvector" << std::endl;
- for (unsigned int i=0; i<indices.size(); ++i)
- deallog << float_values[i] << std::endl;
-}
-
-int main()
-{
- std::ofstream logfile("output");
- deallog.attach(logfile);
- test();
-
- return 0;
-}
+++ /dev/null
-
-DEAL::double_vector_1
-DEAL::0.00000
-DEAL::3.00000
-DEAL::6.00000
-DEAL::9.00000
-DEAL::12.0000
-DEAL::15.0000
-DEAL::18.0000
-DEAL::21.0000
-DEAL::24.0000
-DEAL::27.0000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::0.00000
-DEAL::double_vector_2
-DEAL::0.00000
-DEAL::2.00000
-DEAL::4.00000
-DEAL::6.00000
-DEAL::8.00000
-DEAL::10.0000
-DEAL::12.0000
-DEAL::14.0000
-DEAL::16.0000
-DEAL::18.0000
-DEAL::subvector
-DEAL::0.00000
-DEAL::3.00000
-DEAL::6.00000
-DEAL::9.00000
-DEAL::12.0000
-DEAL::15.0000
-DEAL::18.0000
-DEAL::21.0000
-DEAL::24.0000
-DEAL::27.0000
+++ /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.
-//
-// ---------------------------------------------------------------------
-
-// Check constructor and operator=
-
-#include "../tests.h"
-#include <deal.II/base/index_set.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/read_write_vector.h>
-
-#include <fstream>
-
-
-void test()
-{
- unsigned int double_size = 2;
- unsigned int float_size = 10;
- IndexSet is(50);
- is.add_range(0,2);
- is.add_index(46);
- is.add_range(10,15);
- LinearAlgebra::ReadWriteVector<double> double_vector(is);
- LinearAlgebra::ReadWriteVector<float> float_vector(float_size);
- deallog << "double_size " << double_vector.n_elements() <<std::endl;
- deallog << "float_size " << float_vector.n_elements() <<std::endl;
-
- double_vector = 0.;
- for (unsigned int i=0; i<double_vector.n_elements(); ++i)
- double_vector.local_element(i) += i;
- for (unsigned int i=0; i<float_vector.n_elements(); ++i)
- float_vector[i] = i;
-
- double_vector.print(deallog.get_file_stream());
- float_vector.print(deallog.get_file_stream());
-
- float_vector = double_vector;
- float_vector.print(deallog.get_file_stream());
-
- LinearAlgebra::ReadWriteVector<double> double_vector2(double_vector);
- double_vector2.print(deallog.get_file_stream());
-
- for (unsigned int i=0; i<double_vector.n_elements(); ++i)
- double_vector2.local_element(i) += i;
- double_vector = double_vector2;
- double_vector.print(deallog.get_file_stream());
-}
-
-int main()
-{
- std::ofstream logfile("output");
- deallog.attach(logfile);
- test();
-
- return 0;
-}
+++ /dev/null
-
-DEAL::double_size 8
-DEAL::float_size 10
-IndexSet: {[0,1], [10,14], 46}
-
-0.000e+00
-1.000e+00
-2.000e+00
-3.000e+00
-4.000e+00
-5.000e+00
-6.000e+00
-7.000e+00
-IndexSet: {[0,9]}
-
-0.000e+00
-1.000e+00
-2.000e+00
-3.000e+00
-4.000e+00
-5.000e+00
-6.000e+00
-7.000e+00
-8.000e+00
-9.000e+00
-IndexSet: {[0,1], [10,14], 46}
-
-0.000e+00
-1.000e+00
-2.000e+00
-3.000e+00
-4.000e+00
-5.000e+00
-6.000e+00
-7.000e+00
-IndexSet: {[0,1], [10,14], 46}
-
-0.000e+00
-1.000e+00
-2.000e+00
-3.000e+00
-4.000e+00
-5.000e+00
-6.000e+00
-7.000e+00
-IndexSet: {[0,1], [10,14], 46}
-
-0.000e+00
-2.000e+00
-4.000e+00
-6.000e+00
-8.000e+00
-1.000e+01
-1.200e+01
-1.400e+01
+++ /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.
-//
-// ---------------------------------------------------------------------
-
-// Check constructor and operator=
-
-#include "../tests.h"
-#include <deal.II/base/index_set.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/read_write_vector.h>
-
-#include <fstream>
-
-
-void test()
-{
- unsigned int double_size = 2;
- unsigned int float_size = 10;
- IndexSet is(50);
- is.add_range(0,2);
- is.add_index(46);
- is.add_range(10,15);
- LinearAlgebra::ReadWriteVector<double> double_vector(is);
- LinearAlgebra::ReadWriteVector<float> float_vector(float_size);
- deallog << "double_size " << double_vector.n_elements() <<std::endl;
- deallog << "float_size " << float_vector.n_elements() <<std::endl;
-
- double_vector = 0.;
- for (unsigned int i=0; i<double_vector.n_elements(); ++i)
- double_vector.local_element(i) += i;
- for (unsigned int i=0; i<float_vector.n_elements(); ++i)
- float_vector[i] = i;
-
- double_vector.print(deallog.get_file_stream());
- float_vector.print(deallog.get_file_stream());
-
- float_vector = double_vector;
- float_vector.print(deallog.get_file_stream());
-
- LinearAlgebra::ReadWriteVector<double> double_vector2(double_vector);
- double_vector2.print(deallog.get_file_stream());
-
- for (unsigned int i=0; i<double_vector.n_elements(); ++i)
- double_vector2.local_element(i) += i;
- double_vector = double_vector2;
- double_vector.print(deallog.get_file_stream());
-}
-
-int main()
-{
- std::ofstream logfile("output");
- deallog.attach(logfile);
- test();
-
- return 0;
-}
+++ /dev/null
-
-DEAL::double_size 8
-DEAL::float_size 10
-IndexSet: {[0,1], [10,14], 46}
-
-0.000e+00
-1.000e+00
-2.000e+00
-3.000e+00
-4.000e+00
-5.000e+00
-6.000e+00
-7.000e+00
-IndexSet: {[0,9]}
-
-0.000e+00
-1.000e+00
-2.000e+00
-3.000e+00
-4.000e+00
-5.000e+00
-6.000e+00
-7.000e+00
-8.000e+00
-9.000e+00
-IndexSet: {[0,1], [10,14], 46}
-
-0.000e+00
-1.000e+00
-2.000e+00
-3.000e+00
-4.000e+00
-5.000e+00
-6.000e+00
-7.000e+00
-IndexSet: {[0,1], [10,14], 46}
-
-0.000e+00
-1.000e+00
-2.000e+00
-3.000e+00
-4.000e+00
-5.000e+00
-6.000e+00
-7.000e+00
-IndexSet: {[0,1], [10,14], 46}
-
-0.000e+00
-2.000e+00
-4.000e+00
-6.000e+00
-8.000e+00
-1.000e+01
-1.200e+01
-1.400e+01
+++ /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.
-//
-// ---------------------------------------------------------------------
-
-// test operator =
-
-#include "../tests.h"
-#include <deal.II/base/index_set.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/read_write_vector.h>
-#include <deal.II/base/utilities.h>
-#include <deal.II/lac/trilinos_vector.h>
-
-#include <fstream>
-#include <vector>
-
-void test()
-{
- IndexSet is(8);
- unsigned int rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
- if (rank==0)
- is.add_range(0,4);
- if (rank==1)
- is.add_range(4,8);
- is.compress();
- TrilinosWrappers::MPI::Vector tril_vector(is);
- Vector<double> tmp(8);
- for (unsigned int i=0; i<8; ++i)
- tmp[i] = i;
- tril_vector = tmp;
-
- tril_vector.compress(VectorOperation::insert);
-
- MPI_Barrier(MPI_COMM_WORLD);
-
- IndexSet readwrite_is(8);
- if (rank==0)
- {
- readwrite_is.add_range(0,2);
- readwrite_is.add_range(6,8);
- }
- if (rank==1)
- readwrite_is.add_range(2,6);
- readwrite_is.compress();
- LinearAlgebra::ReadWriteVector<double> readwrite(readwrite_is);
- readwrite = tril_vector;
- if (rank==0)
- {
- std::vector<double> comp(4);
- comp[0] = 0.;
- comp[1] = 1.;
- comp[2] = 6.;
- comp[3] = 7.;
- for (unsigned int i=0; i<4; ++i)
- AssertThrow(readwrite.local_element(i)==comp[i],
- ExcMessage("Element not copied correctly"));
- }
- MPI_Barrier(MPI_COMM_WORLD);
- if (rank==1)
- {
- std::vector<double> comp(4);
- comp[0] = 2.;
- comp[1] = 3.;
- comp[2] = 4.;
- comp[3] = 5.;
- for (unsigned int i=0; i<4; ++i)
- AssertThrow(readwrite.local_element(i)==comp[i],
- ExcMessage("Element not copied correctly"));
- }
- deallog << "OK" <<std::endl;
-}
-
-int main (int argc, char **argv)
-{
- std::ofstream logfile("output");
- deallog.attach(logfile);
- deallog.threshold_double(1.e-10);
-
- Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads());
- test();
-}
+++ /dev/null
-
-DEAL::OK
projects / dealii.git / commitdiff