* vec.locally_owned_elements() == complete_index_set(vec.size())
* @endcode
*/
- const dealii::IndexSet &locally_owned_elements() const;
+ const dealii::IndexSet locally_owned_elements() const;
/**
* Prints the vector to the output stream @p out.
template <typename Number>
inline
- const dealii::IndexSet &Vector<Number>::locally_owned_elements() const
+ const dealii::IndexSet Vector<Number>::locally_owned_elements() const
{
return ReadWriteVector<Number>::get_stored_elements();
}
#include <cstring>
#include <iomanip>
+#ifdef DEAL_II_WITH_TRILINOS
+#include "Epetra_MultiVector.h"
+#endif
DEAL_II_NAMESPACE_OPEN
class Vector;
}
}
+
+namespace LinearAlgebra
+{
+ namespace EpetraWrappers
+ {
+ class Vector;
+ }
+}
#endif
namespace LinearAlgebra
*/
ReadWriteVector<Number> &
operator= (const TrilinosWrappers::MPI::Vector &trilinos_vec);
+
+ /**
+ * Imports all the elements present in the vector's IndexSet from the input
+ * vector @p epetra_vec.
+ */
+ ReadWriteVector<Number> &
+ operator= (const EpetraWrappers::Vector &epetra_vec);
#endif
/**
//@}
protected:
+#ifdef DEAL_II_WITH_TRILINOS
+ /**
+ * Imports all the elements present in the vector's IndexSet from the input
+ * vector @p multivector.
+ */
+ ReadWriteVector<Number> &import(const Epetra_MultiVector &multivector,
+ const IndexSet &locally_owned_elements);
+#endif
/**
* Return the local position of @p global_index.
#include <deal.II/lac/read_write_vector.h>
#include <deal.II/lac/petsc_parallel_vector.h>
#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_epetra_vector.h>
#ifdef DEAL_II_WITH_TRILINOS
# include <Epetra_Import.h>
#ifdef DEAL_II_WITH_TRILINOS
template <typename Number>
ReadWriteVector<Number> &
- ReadWriteVector<Number>::operator = (const TrilinosWrappers::MPI::Vector &trilinos_vec)
+ ReadWriteVector<Number>::import(const Epetra_MultiVector &multivector,
+ const IndexSet &locally_owned_elements)
{
// Copy the local elements
std::vector<size_type> trilinos_indices, readwrite_indices;
- trilinos_vec.locally_owned_elements().fill_index_vector(trilinos_indices);
+ 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;
++j;
trilinos_position[i] = j;
}
- double *values = trilinos_vec.trilinos_vector().Values();
+ double *values = multivector.Values();
for (size_t i=0; i<trilinos_position.size(); ++i)
val[global_to_local(intersection[i])] = values[trilinos_position[i]];
// Copy the off-processor elements if necessary
if (intersection_size != n_elements())
{
- Epetra_BlockMap source_map = trilinos_vec.trilinos_vector().Map();
+ Epetra_BlockMap source_map = multivector.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(),
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);
+ target_vector.Import(multivector, 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
+ template <typename Number>
+ ReadWriteVector<Number> &
+ ReadWriteVector<Number>::operator = (const TrilinosWrappers::MPI::Vector &trilinos_vec)
+ {
+ return import(trilinos_vec.trilinos_vector(), trilinos_vec.locally_owned_elements());
+ }
+
+
+ template <typename Number>
+ ReadWriteVector<Number> &
+ ReadWriteVector<Number>::operator = (const LinearAlgebra::EpetraWrappers::Vector &trilinos_vec)
+ {
+ return import(trilinos_vec.trilinos_vector(), trilinos_vec.locally_owned_elements());
+ }
+
template <typename Number>
void
--- /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__trilinos_epetra_vector_h
+#define dealii__trilinos_epetra_vector_h
+
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+#ifdef DEAL_II_WITH_MPI
+
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/std_cxx11/shared_ptr.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/vector_space_vector.h>
+#include "mpi.h"
+#include "Epetra_FEVector.h"
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+ // Forward declaration
+ template <typename Number>
+ class ReadWriteVector;
+
+ namespace EpetraWrappers
+ {
+ /**
+ * This class implements a wrapper to the the Trilinos distributed vector
+ * class Epetra_FEVector. This class is derived from the
+ * LinearAlgebra::VectorSpaceVector class. Note however that Epetra only
+ * works with Number = double.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Vectors
+ * @author Bruno Turcksin, 2015
+ */
+ class Vector : public VectorSpaceVector<double>
+ {
+ public:
+ /**
+ * Constructor. Create a vector of dimension zero.
+ */
+ Vector();
+
+ /**
+ * Copy constructor. Sets the dimension and the partitioning to that of
+ * the given vector and copies all elements.
+ */
+ Vector(const Vector &V);
+
+ /**
+ * This constructor takes an IndexSet that defines how to distribute the
+ * individual components among the MPI processors. Since it also
+ * includes information abtou the size of the vector, this is all we
+ * need to genereate a %parallel vector.
+ */
+ Vector(const IndexSet ¶llel_partitioner,
+ const MPI_Comm &communicator = MPI_COMM_WORLD);
+
+ /**
+ * Reinit functionality. This function destroys the old vector content
+ * and generates a new one based on the input partitioning. The flag
+ * <tt>omit_zeroing_entries</tt> determines whether the vector should be
+ * filled with zero (false) or left untouched (true).
+ */
+ void reinit (const IndexSet ¶llel_partitioner,
+ const MPI_Comm &communicator = MPI_COMM_WORLD,
+ const bool omit_zeroing_entries = false);
+
+ /**
+ * Copy function. This function takes a Vector and copies all the
+ * elements. The Vector will have the same parallel distribution as @p
+ * V.
+ */
+ Vector &operator= (const Vector &V);
+
+ /**
+ * Copy function. Copy the elements from the ReadWriteVector @p V in the
+ * Vector. If elements if @p V are not locally owned, there is a
+ * communication.
+ */
+ Vector &operator= (const ::dealii::LinearAlgebra::ReadWriteVector<double> &V);
+
+ /**
+ * Multiply the entire vector by a fixed factor.
+ */
+ virtual VectorSpaceVector<double> &operator*= (const double factor);
+
+ /**
+ * Divide the entire vector by a fixed factor.
+ */
+ virtual VectorSpaceVector<double> &operator/= (const double factor);
+
+ /**
+ * Add the vector @p V to the present one.
+ */
+ virtual VectorSpaceVector<double> &operator+= (const VectorSpaceVector<double> &V);
+
+ /**
+ * Substract the vector @p V from the present one.
+ */
+ virtual VectorSpaceVector<double> &operator-= (const VectorSpaceVector<double> &V);
+
+ /**
+ * Return the scalar product of two vectors. The vectors need to have the
+ * same layout.
+ */
+ virtual double operator* (const VectorSpaceVector<double> &V);
+
+ /**
+ * Add @p a to all components. Note that @p is a scalar not a vector.
+ */
+ virtual void add(const double a);
+
+ /**
+ * Simple addition of a multiple of a vector, i.e. <tt>*this +=
+ * a*V</tt>. The vectors need to have the same layout.
+ */
+ virtual void add(const double a, const VectorSpaceVector<double> &V);
+
+ /**
+ * Multiple addition of multiple of a vector, i.e. <tt>*this> +=
+ * a*V+b*W</tt>. The vectors need to have the same layout.
+ */
+ virtual void add(const double a, const VectorSpaceVector<double> &V,
+ const double b, const VectorSpaceVector<double> &W);
+
+ /**
+ * Scaling and simple addition of a multiple of a vector, i.e. <tt>*this
+ * = s*(*this)+a*V</tt>.
+ */
+ virtual void sadd(const double s, const double a,
+ const VectorSpaceVector<double> &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. The
+ * vectors need to have the same layout.
+ */
+ virtual void scale(const VectorSpaceVector<double> &scaling_factors);
+
+ /**
+ * Assignement <tt>*this = a*V</tt>.
+ */
+ virtual void equ(const double a, const VectorSpaceVector<double> &V);
+
+
+ /**
+ * Returns the l<sub>1</sub> norm of the vector (i.e., the sum of the
+ * absolute values of all entries among all processors).
+ */
+ virtual double l1_norm();
+
+ /**
+ * Returns the l<sub>2</sub> norm of the vector (i.e., the square root of
+ * the sum of the square of all entries among all processors).
+ */
+ virtual double l2_norm();
+
+ /**
+ * Returns the maximum norm of the vector (i.e., the maximum absolute value
+ * among all entries and among all processors).
+ */
+ virtual double linfty_norm();
+
+ /**
+ * Performs a combined operation of a vector addition and a subsequent
+ * inner product, returning the value of the inner product. In other
+ * words, the result of this function is the same as if the user called
+ * @code
+ * this->add(a, V);
+ * return_value = *this * W;
+ * @endcode
+ *
+ * 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 double add_and_dot(const double a,
+ const VectorSpaceVector<double> &V,
+ const VectorSpaceVector<double> &W);
+
+ /**
+ * Returns 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;
+
+ /**
+ * Return the MPI communicator object in use with this object.
+ */
+ MPI_Comm get_mpi_communicator() 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
+ */
+ virtual const ::dealii::IndexSet locally_owned_elements() const;
+
+ /**
+ * Return a const reference to the underlying Trilinos
+ * Epetra_FEVector class.
+ */
+ const Epetra_FEVector &trilinos_vector() const;
+
+ /**
+ * Return a (modifyable) reference to the underlying Trilinos
+ * Epetra_FEVector class.
+ */
+ Epetra_FEVector &trilinos_vector();
+
+ /**
+ * 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;
+
+ /**
+ * Returns the memory consumption of this class in bytes.
+ */
+ virtual std::size_t memory_consumption() const;
+
+ /**
+ * The vectors have different partitioning, i.e. they have use different
+ * IndexSet.
+ */
+ DeclException0(ExcDifferentParallelPartitioning);
+
+ /**
+ * Attempt to perform an operation between two incompatible vector types.
+ *
+ * @ingroup Exceptions
+ */
+ DeclException0(ExcVectorTypeNotCompatible);
+
+ /**
+ * Exception thrown by an error in Trilinos.
+ *
+ * @ingroup Exceptions
+ */
+ DeclException1 (ExcTrilinosError,
+ int,
+ << "An error with error number " << arg1
+ << " occurred while calling a Trilinos function");
+
+ private:
+ /**
+ * Pointer to the actual Epetra vector object.
+ */
+ std_cxx11::shared_ptr<Epetra_FEVector> vector;
+ };
+ }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+
+#endif
+
+#endif
* vec.locally_owned_elements() == complete_index_set(vec.size())
* @endcode
*/
- virtual const dealii::IndexSet &locally_owned_elements() const = 0;
+ virtual const dealii::IndexSet locally_owned_elements() const = 0;
/**
* Print the vector to the output stream @p out.
${_src}
trilinos_block_sparse_matrix.cc
trilinos_block_vector.cc
+ trilinos_epetra_vector.cc
trilinos_parallel_block_vector.cc
trilinos_precondition.cc
trilinos_precondition_ml.cc
--- /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/trilinos_epetra_vector.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+#ifdef DEAL_II_WITH_MPI
+
+#include <deal.II/base/index_set.h>
+#include "Epetra_Import.h"
+#include "Epetra_Map.h"
+#include "Epetra_MpiComm.h"
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+ namespace EpetraWrappers
+ {
+ Vector::Vector()
+ :
+ vector(new Epetra_FEVector(Epetra_Map(0,1,0,Epetra_MpiComm(MPI_COMM_SELF))))
+ {}
+
+
+
+ Vector::Vector(const Vector &V)
+ :
+ vector(new Epetra_FEVector(V.trilinos_vector()))
+ {}
+
+
+
+ Vector::Vector(const IndexSet ¶llel_partitioner,
+ const MPI_Comm &communicator)
+ :
+ vector(new Epetra_FEVector(parallel_partitioner.make_trilinos_map(communicator,false)))
+ {}
+
+
+
+ void Vector::reinit(const IndexSet ¶llel_partitioner,
+ const MPI_Comm &communicator,
+ const bool omit_zeroing_entries)
+ {
+ Epetra_Map input_map = parallel_partitioner.make_trilinos_map(communicator,false);
+ if (vector->Map().SameAs(input_map)==false)
+ vector.reset(new Epetra_FEVector(input_map));
+ else if (omit_zeroing_entries==false)
+ {
+ const int ierr = vector->PutScalar(0.);
+ (void) ierr;
+ Assert(ierr==0, ExcTrilinosError(ierr));
+ }
+ }
+
+
+
+ Vector &Vector::operator= (const Vector &V)
+ {
+ // Distinguish three cases:
+ // - First case: both vectors have the same layout.
+ // - Second case: both vectors have the same size but different layout.
+ // - Third case: the vectors have different size.
+ if (vector->Map().SameAs(V.trilinos_vector().Map()))
+ *vector = V.trilinos_vector();
+ else
+ {
+ if (size()==V.size())
+ {
+ Epetra_Import data_exchange(vector->Map(), V.trilinos_vector().Map());
+
+ const int ierr = vector->Import(V.trilinos_vector(), data_exchange, Insert);
+ Assert(ierr==0, ExcTrilinosError(ierr));
+ }
+ else
+ vector.reset(new Epetra_FEVector(V.trilinos_vector()));
+ }
+
+ return *this;
+ }
+
+
+
+ Vector &Vector::operator= (const ::dealii::LinearAlgebra::ReadWriteVector<double> &V)
+ {
+ IndexSet elements_to_write(V.get_stored_elements());
+ IndexSet locally_owned_elements(this->locally_owned_elements());
+ IndexSet off_proc_elements(elements_to_write);
+ off_proc_elements.subtract_set(locally_owned_elements);
+ IndexSet on_proc_elements(elements_to_write);
+ on_proc_elements.subtract_set(off_proc_elements);
+
+ // Write the elements that locally owned.
+ IndexSet::ElementIterator elem = on_proc_elements.begin(),
+ elem_end = on_proc_elements.end();
+ for (; elem!=elem_end; ++elem)
+ {
+ const TrilinosWrappers::types::int_type local_index =
+ vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(*elem));
+ (*vector)[0][local_index] = V[*elem];
+ }
+
+
+ // Write the elements off-processor. Cannot use Import function of Trilinos
+ // because V is not a Trilinos object. The most straightforward to send
+ // the off-processor to the other processors is to AllGather. This way
+ // every processor has access to all the off-processor elements. The
+ // problem with this method is that if every element in V is
+ // off-processor, it is possible that the buffers used to receive the data
+ // has twice the global size of the current Vector.
+ // TODO This won't scale past a few 10,000's of processors.
+
+ // Get recv_count and the size of the buffer needed to receive the data.
+ const Epetra_MpiComm *epetra_comm
+ = dynamic_cast<const Epetra_MpiComm *>(&(vector->Comm()));
+ MPI_Comm comm = epetra_comm->GetMpiComm();
+ int send_buffer_size(off_proc_elements.n_elements());
+ int n_procs(0);
+ MPI_Comm_size(comm, &n_procs);
+ std::vector<int> recv_count(n_procs);
+ MPI_Allgather(&send_buffer_size, 1, MPI_INT, &recv_count, 1,
+ MPI_INT, comm);
+ std::vector<int> displacement(n_procs,0);
+ for (int i=1; i<n_procs; ++i)
+ displacement[i] = displacement[i-1] + recv_count[i-1];
+ const unsigned int recv_buffer_size = displacement[n_procs-1] +
+ recv_count[n_procs-1];
+
+ // Write the indices in the send buffer
+ std::vector<types::global_dof_index> send_index_buffer(send_buffer_size);
+ elem = off_proc_elements.begin();
+ for (int i=0; i<send_buffer_size; ++i)
+ {
+ send_index_buffer[i] = *elem;
+ ++elem;
+ }
+ // Send the index of the off-processor elements
+ std::vector<types::global_dof_index> recv_index_buffer(recv_buffer_size);
+ MPI_Allgatherv(&send_index_buffer[0], send_buffer_size, DEAL_II_DOF_INDEX_MPI_TYPE,
+ &recv_index_buffer[0], &recv_count[0], &displacement[0],
+ DEAL_II_DOF_INDEX_MPI_TYPE, comm);
+
+ // Write the values in the send buffer
+ std::vector<double> send_value_buffer(send_buffer_size);
+ elem = off_proc_elements.begin();
+ for (int i=0; i<send_buffer_size; ++i)
+ {
+ send_value_buffer[i] = V[*elem];
+ ++elem;
+ }
+ // Send the value of the off-processor elements
+ std::vector<double> recv_value_buffer(recv_buffer_size);
+ MPI_Allgatherv(&send_value_buffer[0], send_buffer_size, MPI_DOUBLE,
+ &recv_value_buffer[0], &recv_count[0], &displacement[0], MPI_DOUBLE, comm);
+
+ // Write the data in the vector
+ for (unsigned int i=0; i<recv_buffer_size; ++i)
+ {
+ if (locally_owned_elements.is_element(recv_index_buffer[i]))
+ {
+ const TrilinosWrappers::types::int_type local_index =
+ vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(
+ recv_index_buffer[i]));
+ (*vector)[0][local_index] = recv_value_buffer[i];
+ }
+ }
+
+ return *this;
+ }
+
+
+
+ VectorSpaceVector<double> &Vector::operator*= (const double factor)
+ {
+ AssertIsFinite(factor);
+ vector->Scale(factor);
+
+ return *this;
+ }
+
+
+
+ VectorSpaceVector<double> &Vector::operator/= (const double factor)
+ {
+ AssertIsFinite(factor);
+ Assert(factor!=0., ExcZero());
+ *this *= 1./factor;
+
+ return *this;
+ }
+
+
+
+ VectorSpaceVector<double> &Vector::operator+= (const VectorSpaceVector<double> &V)
+ {
+ // Check that casting will work.
+ Assert(dynamic_cast<const Vector *>(&V)!=NULL, ExcVectorTypeNotCompatible());
+
+ // Downcast V. If fails, throws an exception.
+ const Vector &down_V = dynamic_cast<const Vector &>(V);
+ // If the maps are the same we can Update right away.
+ if (vector->Map().SameAs(down_V.trilinos_vector().Map()))
+ vector->Update(1., down_V.trilinos_vector(), 1.);
+ else
+ {
+ Assert(this->size()==down_V.size(),
+ ExcDimensionMismatch(this->size(), down_V.size()));
+
+#if DEAL_II_TRILINOS_VERSION_GTE(11,11,0)
+ Epetra_Import data_exchange (vector->Map(), down_V.trilinos_vector().Map());
+ int ierr = vector->Import(down_V.trilinos_vector(), data_exchange, Epetra_AddLocalAlso);
+ (void) ierr;
+ Assert(ierr==0, ExcTrilinosError(ierr));
+#else
+ // In versions older than 11.11 the Import function is broken for adding
+ // Hence, we provide a workaround in this case
+
+ Epetra_MultiVector dummy(vector->Map(), 1, false);
+ Epetra_Import data_exchange(dummy.Map(), down_V.trilinos_vector().Map());
+
+ int ierr = dummy.Import(down_V.trilinos_vector(), data_exchange, Insert);
+ (void) ierr;
+ Assert(ierr==0, ExcTrilinosError(ierr));
+
+ ierr = vector->Update(1.0, dummy, 1.0);
+ (void) ierr;
+ Assert(ierr==0, ExcTrilinosError(ierr));
+#endif
+ }
+
+ return *this;
+ }
+
+
+
+ VectorSpaceVector<double> &Vector::operator-= (const VectorSpaceVector<double> &V)
+ {
+ this->add(-1.,V);
+
+ return *this;
+ }
+
+
+
+ double Vector::operator* (const VectorSpaceVector<double> &V)
+ {
+ // Check that casting will work.
+ Assert(dynamic_cast<const Vector *>(&V)!=NULL,
+ ExcVectorTypeNotCompatible());
+
+ // Downcast V. If fails, throws an exception.
+ const Vector &down_V = dynamic_cast<const Vector &>(V);
+ Assert(this->size()==down_V.size(),
+ ExcDimensionMismatch(this->size(), down_V.size()));
+ Assert(vector->Map().SameAs(down_V.trilinos_vector().Map()),
+ ExcDifferentParallelPartitioning());
+
+ double result(0.);
+ const int ierr = vector->Dot(down_V.trilinos_vector(), &result);
+ (void) ierr;
+ Assert(ierr==0, ExcTrilinosError(ierr));
+
+ return result;
+ }
+
+
+
+ void Vector::add(const double a)
+ {
+ AssertIsFinite(a);
+ const unsigned local_size(vector->MyLength());
+ for (unsigned int i=0; i<local_size; ++i)
+ (*vector)[0][i] += a;
+ }
+
+
+
+ void Vector::add(const double a, const VectorSpaceVector<double> &V)
+ {
+ // Check that casting will work.
+ Assert(dynamic_cast<const Vector *>(&V)!=NULL,
+ ExcVectorTypeNotCompatible());
+
+ // Downcast V. If fails, throws an exception.
+ const Vector &down_V = dynamic_cast<const Vector &>(V);
+ AssertIsFinite(a);
+ Assert(vector->Map().SameAs(down_V.trilinos_vector().Map()),
+ ExcDifferentParallelPartitioning());
+
+ const int ierr = vector->Update(a, down_V.trilinos_vector(), 1.);
+ Assert(ierr==0, ExcTrilinosError(ierr));
+ }
+
+
+
+ void Vector::add(const double a, const VectorSpaceVector<double> &V,
+ const double b, const VectorSpaceVector<double> &W)
+ {
+ // Check that casting will work.
+ Assert(dynamic_cast<const Vector *>(&V)!=NULL,
+ ExcVectorTypeNotCompatible());
+ // Check that casting will work.
+ Assert(dynamic_cast<const Vector *>(&W)!=NULL,
+ ExcVectorTypeNotCompatible());
+
+ // Downcast V. If fails, throws an exception.
+ const Vector &down_V = dynamic_cast<const Vector &>(V);
+ // Downcast W. If fails, throws an exception.
+ const Vector &down_W = dynamic_cast<const Vector &>(W);
+ Assert(vector->Map().SameAs(down_V.trilinos_vector().Map()),
+ ExcDifferentParallelPartitioning());
+ Assert(vector->Map().SameAs(down_W.trilinos_vector().Map()),
+ ExcDifferentParallelPartitioning());
+ AssertIsFinite(a);
+ AssertIsFinite(b);
+
+ const int ierr = vector->Update(a, down_V.trilinos_vector(), b,
+ down_W.trilinos_vector(), 1.);
+ Assert(ierr==0, ExcTrilinosError(ierr));
+ }
+
+
+
+ void Vector::sadd(const double s, const double a,
+ const VectorSpaceVector<double> &V)
+ {
+ // Check that casting will work.
+ Assert(dynamic_cast<const Vector *>(&V)!=NULL,
+ ExcVectorTypeNotCompatible());
+
+ *this *= s;
+ // Downcast V. It fails, throws an exception.
+ const Vector &down_V = dynamic_cast<const Vector &>(V);
+ Vector tmp(down_V);
+ tmp *= a;
+ *this += tmp;
+ }
+
+
+
+ void Vector::scale(const VectorSpaceVector<double> &scaling_factors)
+ {
+ // Check that casting will work.
+ Assert(dynamic_cast<const Vector *>(&scaling_factors)!=NULL,
+ ExcVectorTypeNotCompatible());
+
+ // Downcast scaling_factors. If fails, throws an exception.
+ const Vector &down_scaling_factors =
+ dynamic_cast<const Vector &>(scaling_factors);
+ Assert(vector->Map().SameAs(down_scaling_factors.trilinos_vector().Map()),
+ ExcDifferentParallelPartitioning());
+
+ const int ierr = vector->Multiply(1.0, down_scaling_factors.trilinos_vector(),
+ *vector, 0.0);
+ Assert(ierr==0, ExcTrilinosError(ierr));
+ }
+
+
+
+ void Vector::equ(const double a, const VectorSpaceVector<double> &V)
+ {
+ // Check that casting will work.
+ Assert(dynamic_cast<const Vector *>(&V)!=NULL,
+ ExcVectorTypeNotCompatible());
+
+ // Downcast V. If fails, throws an exception.
+ const Vector &down_V = dynamic_cast<const Vector &>(V);
+ // If we don't have the same map, copy.
+ if (vector->Map().SameAs(down_V.trilinos_vector().Map())==false)
+ this->sadd(0., a, V);
+ else
+ {
+ // Otherwise, just update
+ int ierr = vector->Update(a, down_V.trilinos_vector(), 0.);
+ Assert(ierr==0, ExcTrilinosError(ierr));
+ }
+ }
+
+
+
+ double Vector::l1_norm()
+ {
+ double norm(0.);
+ int ierr = vector->Norm1(&norm);
+ Assert(ierr==0, ExcTrilinosError(ierr));
+
+ return norm;
+ }
+
+
+
+ double Vector::l2_norm()
+ {
+ double norm(0.);
+ int ierr = vector->Norm2(&norm);
+ Assert(ierr==0, ExcTrilinosError(ierr));
+
+ return norm;
+ }
+
+
+
+ double Vector::linfty_norm()
+ {
+ double norm(0.);
+ int ierr = vector->NormInf(&norm);
+ Assert(ierr==0, ExcTrilinosError(ierr));
+
+ return norm;
+ }
+
+
+
+ double Vector::add_and_dot(const double a,
+ const VectorSpaceVector<double> &V,
+ const VectorSpaceVector<double> &W)
+ {
+ this->add(a, V);
+
+ return *this * W;
+ }
+
+
+
+ Vector::size_type Vector::size() const
+ {
+#ifndef DEAL_II_WITH_64BIT_INDICES
+ return vector->GlobalLength();
+#else
+ return vector->GlobalLength64();
+#endif
+ }
+
+
+
+ MPI_Comm Vector::get_mpi_communicator() const
+ {
+ const Epetra_MpiComm *epetra_comm
+ = dynamic_cast<const Epetra_MpiComm *>(&(vector->Comm()));
+ return epetra_comm->GetMpiComm();
+ }
+
+
+
+ const ::dealii::IndexSet Vector::locally_owned_elements() const
+ {
+ const Epetra_Map *map = dynamic_cast<const Epetra_Map *>(&(vector->Map()));
+ return IndexSet(*map);
+ }
+
+
+
+ const Epetra_FEVector &Vector::trilinos_vector() const
+ {
+ return *vector;
+ }
+
+
+
+ Epetra_FEVector &Vector::trilinos_vector()
+ {
+ return *vector;
+ }
+
+
+
+ void Vector::print(std::ostream &out,
+ const unsigned int precision,
+ const bool scientific,
+ const bool across) const
+ {
+ AssertThrow(out, ExcIO());
+
+ // Get a representation of the vector and loop over all
+ // the elements
+ double *val;
+ int leading_dimension;
+ int ierr = vector->ExtractView(&val, &leading_dimension);
+
+ Assert(ierr==0, ExcTrilinosError(ierr));
+ out.precision (precision);
+ if (scientific)
+ out.setf(std::ios::scientific, std::ios::floatfield);
+ else
+ out.setf(std::ios::fixed, std::ios::floatfield);
+
+ if (across)
+ for (size_type i=0; i<size(); ++i)
+ out << val[i] << ' ';
+ else
+ for (size_type i=0; i<size(); ++i)
+ out << val[i] << std::endl;
+ out << std::endl;
+
+ // restore the representation
+ // of the vector
+ AssertThrow(out, ExcIO());
+ }
+
+
+
+ std::size_t Vector::memory_consumption() const
+ {
+ return sizeof(*this)
+ + vector->MyLength()*(sizeof(double)+
+ sizeof(TrilinosWrappers::types::int_type));
+ }
+ }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+
+#endif
// ---------------------------------------------------------------------
-// checks that the l1, l2, lp norm is computed correctly for deal.II vectors
+// checks that the l1 and l2 norm are 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)
--- /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 "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/trilinos_epetra_vector.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+// Check LinearAlgebra::EpetraWrappers::Vector assignement and the operator
+// overloading.
+
+
+void test()
+{
+ IndexSet parallel_partitioner_1(10);
+ IndexSet parallel_partitioner_2(10);
+ unsigned int rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+ if (rank==0)
+ {
+ parallel_partitioner_1.add_range(0,5);
+ parallel_partitioner_2.add_range(0,3);
+ }
+ else
+ {
+ parallel_partitioner_1.add_range(5,10);
+ parallel_partitioner_2.add_range(3,10);
+ }
+ parallel_partitioner_1.compress();
+ parallel_partitioner_2.compress();
+ LinearAlgebra::EpetraWrappers::Vector a;
+ LinearAlgebra::EpetraWrappers::Vector b(parallel_partitioner_1);
+ LinearAlgebra::EpetraWrappers::Vector c(b);
+
+ AssertThrow(a.size()==0, ExcMessage("Vector has the wrong size."));
+ AssertThrow(b.size()==10, ExcMessage("Vector has the wrong size."));
+ AssertThrow(c.size()==10, ExcMessage("Vector has the wrong size."));
+
+ a.reinit(parallel_partitioner_2);
+ AssertThrow(a.size()==10, ExcMessage("Vector has the wrong size."));
+
+ AssertThrow(parallel_partitioner_1==b.locally_owned_elements(),
+ ExcMessage("IndexSet has been modified."));
+ AssertThrow(parallel_partitioner_2==a.locally_owned_elements(),
+ ExcMessage("IndexSet has been modified."));
+
+ IndexSet read_write_index_set(10);
+ if (rank==0)
+ read_write_index_set.add_range(0,5);
+ else
+ read_write_index_set.add_range(5,10);
+ read_write_index_set.compress();
+
+ LinearAlgebra::ReadWriteVector<double> read_write_1(read_write_index_set);
+ LinearAlgebra::ReadWriteVector<double> read_write_2(read_write_index_set);
+ LinearAlgebra::ReadWriteVector<double> read_write_3(read_write_index_set);
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ {
+ read_write_1[i] = i;
+ read_write_2[i] = 5.+i;
+ }
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ {
+ read_write_1[i] = i;
+ read_write_2[i] = 5.+i;
+ }
+ }
+
+ a = read_write_2;
+ b = read_write_1;
+ c = read_write_2;
+
+ read_write_3 = a;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(read_write_2[i]==read_write_3[i],
+ ExcMessage("Vector a has been modified."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(read_write_2[i]==read_write_3[i],
+ ExcMessage("Vector a has been modified."));
+ }
+
+ read_write_3 = b;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(read_write_1[i]==read_write_3[i],
+ ExcMessage("Vector b has been modified."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(read_write_1[i]==read_write_3[i],
+ ExcMessage("Vector b has been modified."));
+ }
+
+ read_write_3 = c;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(read_write_2[i]==read_write_3[i],
+ ExcMessage("Vector c has been modified."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(read_write_2[i]==read_write_3[i],
+ ExcMessage("Vector c has been modified."));
+ }
+
+
+ a *= 2;
+ read_write_3 = a;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(2.*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in operator *=."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(2.*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in operator *=."));
+ }
+
+ c /= 2.;
+ read_write_3 = c;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(0.5*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in operator /=."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(0.5*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in operator /=."));
+ }
+
+ b += a;
+ read_write_3 = b;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(2.*read_write_2[i]+read_write_1[i]==read_write_3[i],
+ ExcMessage("Problem in operator +=."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(2.*read_write_2[i]+read_write_1[i]==read_write_3[i],
+ ExcMessage("Problem in operator +=."));
+ }
+
+ b -= c;
+ read_write_3 = b;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(1.5*read_write_2[i]+read_write_1[i]==read_write_3[i],
+ ExcMessage("Problem in operator -=."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(1.5*read_write_2[i]+read_write_1[i]==read_write_3[i],
+ ExcMessage("Problem in operator -=."));
+ }
+
+ b = read_write_1;
+ c = read_write_1;
+ const double val = b*c;
+ AssertThrow(val==285., ExcMessage("Problem in operator *."));
+}
+
+
+int main(int argc, char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+ deallog << "OK" <<std::endl;
+
+ return 0;
+}
--- /dev/null
+
+DEAL::OK
--- /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 "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/trilinos_epetra_vector.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test()
+{
+ IndexSet parallel_partitioner_1(10);
+ IndexSet parallel_partitioner_2(10);
+ unsigned int rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+ if (rank==0)
+ {
+ parallel_partitioner_1.add_range(0,5);
+ parallel_partitioner_2.add_range(0,3);
+ }
+ else
+ {
+ parallel_partitioner_1.add_range(5,10);
+ parallel_partitioner_2.add_range(3,10);
+ }
+ parallel_partitioner_1.compress();
+ parallel_partitioner_2.compress();
+ LinearAlgebra::EpetraWrappers::Vector a(parallel_partitioner_1);
+ LinearAlgebra::EpetraWrappers::Vector b(parallel_partitioner_1);
+ LinearAlgebra::EpetraWrappers::Vector c(parallel_partitioner_2);
+
+ IndexSet read_write_index_set(10);
+ if (rank==0)
+ read_write_index_set.add_range(0,5);
+ else
+ read_write_index_set.add_range(5,10);
+ read_write_index_set.compress();
+
+ LinearAlgebra::ReadWriteVector<double> read_write_1(read_write_index_set);
+ LinearAlgebra::ReadWriteVector<double> read_write_2(read_write_index_set);
+ LinearAlgebra::ReadWriteVector<double> read_write_3(read_write_index_set);
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ {
+ read_write_1[i] = i;
+ read_write_2[i] = 5.+i;
+ }
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ {
+ read_write_1[i] = i;
+ read_write_2[i] = 5.+i;
+ }
+ }
+
+ a = read_write_1;
+ b = read_write_2;
+ c = read_write_2;
+
+ a.add(1.);
+ read_write_3 = a;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(1.+read_write_1[i]==read_write_3[i],
+ ExcMessage("Problem in add(scalar)."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(1.+read_write_1[i]==read_write_3[i],
+ ExcMessage("Problem in add(scalar)."));
+ }
+
+ a.add(2.,b);
+ read_write_3 = a;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(1.+read_write_1[i]+2.*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in add(scalar,Vector)."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(1.+read_write_1[i]+2.*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in add(scalar,Vector)."));
+ }
+
+
+ LinearAlgebra::EpetraWrappers::Vector d(a);
+ a.add(2.,b,3.,d);
+ read_write_3 = a;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(3.+read_write_1[i]+8.*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in add(scalar,Vector,scalar,Vector)."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(3.+read_write_1[i]+8.*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in add(scalar,Vector,scalar,Vector)."));
+ }
+
+
+ a = read_write_1;
+ a.sadd(3.,2.,c);
+ read_write_3 = a;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(3.+read_write_1[i]+2.*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in sadd(scalar,scalar,Vector)."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(3.+read_write_1[i]+2.*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in sadd(scalar,scalar,Vector)."));
+ }
+
+
+ a = read_write_1;
+ a.scale(b);
+ read_write_3 = a;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(read_write_1[i]*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in scale."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(read_write_1[i]*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in scale."));
+ }
+
+
+ a.equ(2.,c);
+ read_write_3 = a;
+ if (rank==0)
+ {
+ for (unsigned int i=0; i<5; ++i)
+ AssertThrow(2.*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in scale."));
+ }
+ else
+ {
+ for (unsigned int i=5; i<10; ++i)
+ AssertThrow(2.*read_write_2[i]==read_write_3[i],
+ ExcMessage("Problem in scale."));
+ }
+
+
+ AssertThrow(b.l1_norm()==45., ExcMessage("Problem in l1_norm."));
+
+ const double eps=1e-6;
+ AssertThrow(std::fabs(b.l2_norm()-16.881943016)<eps,
+ ExcMessage("Problem in l2_norm"));
+
+ AssertThrow(b.linfty_norm()==9., ExcMessage("Problem in linfty_norm."));
+
+ const double val = a.add_and_dot(2.,c,b);
+ AssertThrow(val==1530., ExcMessage("Problem in add_and_dot"));
+}
+
+
+int main(int argc, char **argv)
+{
+ std::ofstream logfile("output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+ deallog << "OK" <<std::endl;
+
+ return 0;
+}
--- /dev/null
+
+DEAL::OK