* @ingroup Matrix1
* @author Martin Kronbichler, Wolfgang Bangerth, 2008
*/
- class BlockSparseMatrix : public BlockMatrixBase<TrilinosWrappers::SparseMatrix>
+ class BlockSparseMatrix : public BlockMatrixBase<SparseMatrix>
{
public:
/**
* let $dst = M*src$ with $M$
* being this matrix.
*/
+ void vmult (MPI::BlockVector &dst,
+ const MPI::BlockVector &src) const;
+
+
+ /**
+ * Matrix-vector multiplication:
+ * let $dst = M*src$ with $M$
+ * being this matrix, now applied
+ * to localized block vectors
+ * (works only when run on one
+ * processor).
+ */
void vmult (BlockVector &dst,
- const BlockVector &src) const;
+ const BlockVector &src) const;
/**
* Matrix-vector
* applicable if the matrix has
* only one block column.
*/
+ void vmult (MPI::BlockVector &dst,
+ const MPI::Vector &src) const;
+
+ /**
+ * Matrix-vector
+ * multiplication. Just like the
+ * previous function, but only
+ * applicable if the matrix has
+ * only one block column, now
+ * applied to localized vectors
+ * (works only when run on one
+ * processor).
+ */
void vmult (BlockVector &dst,
- const Vector &src) const;
+ const Vector &src) const;
/**
* Matrix-vector
* applicable if the matrix has
* only one block row.
*/
+ void vmult (MPI::Vector &dst,
+ const MPI::BlockVector &src) const;
+
+ /**
+ * Matrix-vector
+ * multiplication. Just like the
+ * previous function, but only
+ * applicable if the matrix has
+ * only one block row, now
+ * applied to localized vectors
+ * (works only when run on one
+ * processor).
+ */
void vmult (Vector &dst,
- const BlockVector &src) const;
+ const BlockVector &src) const;
/**
* Matrix-vector
* applicable if the matrix has
* only one block.
*/
+ void vmult (MPI::Vector &dst,
+ const MPI::Vector &src) const;
+
+ /**
+ * Matrix-vector
+ * multiplication. Just like the
+ * previous function, but only
+ * applicable if the matrix has
+ * only one block, now applied to
+ * localized vectors (works only
+ * when run on one processor).
+ */
void vmult (Vector &dst,
- const Vector &src) const;
-
+ const Vector &src) const;
+
/**
* Matrix-vector multiplication:
* let $dst = M^T*src$ with $M$
* vmult() but takes the
* transposed matrix.
*/
+ void Tvmult (MPI::BlockVector &dst,
+ const MPI::BlockVector &src) const;
+
+ /**
+ * Matrix-vector multiplication:
+ * let $dst = M^T*src$ with $M$
+ * being this matrix. This
+ * function does the same as
+ * vmult() but takes the
+ * transposed matrix, now applied
+ * to localized Trilinos vectors
+ * (works only when run on one
+ * processor).
+ */
void Tvmult (BlockVector &dst,
const BlockVector &src) const;
-
+
/**
* Matrix-vector
* multiplication. Just like the
* applicable if the matrix has
* only one block row.
*/
+ void Tvmult (MPI::BlockVector &dst,
+ const MPI::Vector &src) const;
+
+ /**
+ * Matrix-vector
+ * multiplication. Just like the
+ * previous function, but only
+ * applicable if the matrix has
+ * only one block row, now
+ * applied to localized Trilinos
+ * vectors (works only when run
+ * on one processor).
+ */
void Tvmult (BlockVector &dst,
const Vector &src) const;
* applicable if the matrix has
* only one block column.
*/
+ void Tvmult (MPI::Vector &dst,
+ const MPI::BlockVector &src) const;
+
+ /**
+ * Matrix-vector
+ * multiplication. Just like the
+ * previous function, but only
+ * applicable if the matrix has
+ * only one block column, now
+ * applied to localized Trilinos
+ * vectors (works only when run
+ * on one processor).
+ */
void Tvmult (Vector &dst,
const BlockVector &src) const;
* applicable if the matrix has
* only one block.
*/
+ void Tvmult (MPI::Vector &dst,
+ const MPI::Vector &src) const;
+
+ /**
+ * Matrix-vector
+ * multiplication. Just like the
+ * previous function, but only
+ * applicable if the matrix has
+ * only one block, now applied to
+ * localized Trilinos vectors
+ * (works only when run on one
+ * processor).
+ */
void Tvmult (Vector &dst,
const Vector &src) const;
int, int, int, int,
<< "The blocks [" << arg1 << ',' << arg2 << "] and ["
<< arg3 << ',' << arg4 << "] have differing row numbers.");
- /**
+
+ /**
* Exception
*/
DeclException4 (ExcIncompatibleColNumbers,
// ------------- inline and template functions -----------------
+ inline
+ void
+ BlockSparseMatrix::vmult (MPI::BlockVector &dst,
+ const MPI::BlockVector &src) const
+ {
+ BaseClass::vmult_block_block (dst, src);
+ }
+
+
+
inline
void
BlockSparseMatrix::vmult (BlockVector &dst,
{
BaseClass::vmult_block_block (dst, src);
}
+
+
+
+ inline
+ void
+ BlockSparseMatrix::vmult (MPI::BlockVector &dst,
+ const MPI::Vector &src) const
+ {
+ BaseClass::vmult_block_nonblock (dst, src);
+ }
{
BaseClass::vmult_block_nonblock (dst, src);
}
+
+
+
+ inline
+ void
+ BlockSparseMatrix::vmult (MPI::Vector &dst,
+ const MPI::BlockVector &src) const
+ {
+ BaseClass::vmult_nonblock_block (dst, src);
+ }
{
BaseClass::vmult_nonblock_block (dst, src);
}
-
+
+
+
+ inline
+ void
+ BlockSparseMatrix::vmult (MPI::Vector &dst,
+ const MPI::Vector &src) const
+ {
+ BaseClass::vmult_nonblock_nonblock (dst, src);
+ }
+
inline
}
+
+ inline
+ void
+ BlockSparseMatrix::Tvmult (MPI::BlockVector &dst,
+ const MPI::BlockVector &src) const
+ {
+ BaseClass::Tvmult_block_block (dst, src);
+ }
+
+
+
inline
void
BlockSparseMatrix::Tvmult (BlockVector &dst,
- const BlockVector &src) const
+ const BlockVector &src) const
{
BaseClass::Tvmult_block_block (dst, src);
}
-
+
+
+
+ inline
+ void
+ BlockSparseMatrix::Tvmult (MPI::BlockVector &dst,
+ const MPI::Vector &src) const
+ {
+ BaseClass::Tvmult_block_nonblock (dst, src);
+ }
+
inline
void
BlockSparseMatrix::Tvmult (BlockVector &dst,
- const Vector &src) const
+ const Vector &src) const
{
BaseClass::Tvmult_block_nonblock (dst, src);
}
-
+
+
+
+ inline
+ void
+ BlockSparseMatrix::Tvmult (MPI::Vector &dst,
+ const MPI::BlockVector &src) const
+ {
+ BaseClass::Tvmult_nonblock_block (dst, src);
+ }
+
inline
void
BlockSparseMatrix::Tvmult (Vector &dst,
- const BlockVector &src) const
+ const BlockVector &src) const
{
BaseClass::Tvmult_nonblock_block (dst, src);
}
-
+
+
+
+ inline
+ void
+ BlockSparseMatrix::Tvmult (MPI::Vector &dst,
+ const MPI::Vector &src) const
+ {
+ BaseClass::Tvmult_nonblock_nonblock (dst, src);
+ }
+
inline
void
BlockSparseMatrix::Tvmult (Vector &dst,
- const Vector &src) const
+ const Vector &src) const
{
BaseClass::Tvmult_nonblock_nonblock (dst, src);
}
-
+
}
DEAL_II_NAMESPACE_CLOSE
#include <lac/trilinos_vector.h>
#include <lac/block_indices.h>
#include <lac/block_vector_base.h>
-#include <lac/block_vector.h>
#include <lac/exceptions.h>
#ifdef DEAL_II_USE_TRILINOS
namespace TrilinosWrappers
{
// forward declaration
+ namespace MPI
+ {
+ class BlockVector;
+ }
class BlockVector;
+ namespace MPI
+ {
/**
* An implementation of block vectors based on the vector class
* implemented in TrilinosWrappers. While the base class provides for
* the distribution of data among
* the MPI processes.
*/
- explicit BlockVector (const std::vector<Epetra_Map> &InputMaps);
+ BlockVector (const std::vector<Epetra_Map> &InputMaps);
/**
* Copy-Constructor. Set all the
* given scalar value.
*/
BlockVector &
- operator = (const value_type s);
+ operator = (const value_type s);
/**
* Copy operator for arguments of
* the same type.
*/
BlockVector &
- operator = (const BlockVector &V);
+ operator = (const BlockVector &V);
/**
* Another copy function. This
* possible number type in the
* deal.II vector.
*/
- template <typename Number>
- BlockVector &
- operator = (const ::dealii::BlockVector<Number> &V);
+ template <typename Number>
+ BlockVector &
+ operator = (const ::dealii::BlockVector<Number> &V);
/**
* Reinitialize the BlockVector to
* If <tt>fast==false</tt>, the vector
* is filled with zeros.
*/
- void reinit (const std::vector<Epetra_Map> &input_maps);
+ void reinit (const std::vector<Epetra_Map> &input_maps,
+ const bool fast = false);
/**
* Change the dimension to that
* the wrong block.
*/
void reinit (const BlockVector &V,
- const bool fast=false);
+ const bool fast = false);
/**
* Change the number of blocks to
DeclException0 (ExcNonMatchingBlockVectors);
};
+
+
/*----------------------- Inline functions ----------------------------------*/
+
+/**
+ * Global function 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.
+ *
+ * @ingroup TrilinosWrappers
+ * @relates TrilinosWrappers::MPI::BlockVector
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ */
inline
- BlockVector &
- BlockVector::operator = (const value_type s)
+ void swap (BlockVector &u,
+ BlockVector &v)
{
- BaseClass::operator = (s);
- return *this;
+ u.swap (v);
}
+ } /* end of namespace MPI */
- inline
- BlockVector &
- BlockVector::operator = (const BlockVector &v)
- {
- if (this->n_blocks() != v.n_blocks())
- reinit(v.n_blocks());
- for (unsigned int i=0; i<this->n_blocks(); ++i)
- this->components[i] = v.block(i);
- collect_sizes();
-
- return *this;
- }
+/**
+ * An implementation of block vectors based on the vector class
+ * implemented in TrilinosWrappers. While the base class provides for
+ * most of the interface, this class handles the actual allocation of
+ * vectors and provides functions that are specific to the underlying
+ * vector type.
+ *
+ * In contrast to the class MPI::BlockVector, this class is based on a
+ * localized version of the vectors, which means that the whole vector
+ * is stored on each processor. Note that matrix vector products with
+ * this block vector class do only work in case the program is run on
+ * only one processor, since the Trilinos matrices are inherently
+ * parallel.
+ *
+ * @ingroup Vectors
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, 2008
+ */
+ class BlockVector : public BlockVectorBase<Vector>
+ {
+ public:
+ /**
+ * Typedef the base class for simpler
+ * access to its own typedefs.
+ */
+ typedef BlockVectorBase<Vector> BaseClass;
+
+ /**
+ * Typedef the type of the underlying
+ * vector.
+ */
+ typedef BaseClass::BlockType BlockType;
+ /**
+ * Import the typedefs from the base
+ * class.
+ */
+ typedef BaseClass::value_type value_type;
+ typedef BaseClass::pointer pointer;
+ typedef BaseClass::const_pointer const_pointer;
+ typedef BaseClass::reference reference;
+ typedef BaseClass::const_reference const_reference;
+ typedef BaseClass::size_type size_type;
+ typedef BaseClass::iterator iterator;
+ typedef BaseClass::const_iterator const_iterator;
+ /**
+ * Default constructor. Generate an
+ * empty vector without any blocks.
+ */
+ BlockVector ();
+
+ /**
+ * Constructor. Generate a block
+ * vector with as many blocks as
+ * there are entries in Input_Maps.
+ * Each Epetra_Map already knows
+ * the distribution of data among
+ * the MPI processes.
+ */
+ BlockVector (const std::vector<Epetra_Map> &InputMaps);
- template <typename Number>
- inline
- BlockVector &
- BlockVector::operator = (const ::dealii::BlockVector<Number> &v)
- {
- Assert (n_blocks() == v.n_blocks(),
- ExcDimensionMismatch(n_blocks(),v.n_blocks()));
+ /**
+ * Copy-Constructor. Set all the
+ * properties of the parallel
+ * vector to those of the given
+ * argument and copy the
+ * elements.
+ */
+ BlockVector (const MPI::BlockVector &V);
+
+ /**
+ * Copy-Constructor. Set all the
+ * properties of the parallel
+ * vector to those of the given
+ * argument and copy the
+ * elements.
+ */
+ BlockVector (const BlockVector &V);
+
+ /**
+ * Creates a block vector
+ * consisting of
+ * <tt>num_blocks</tt>
+ * components, but there is no
+ * content in the individual
+ * components and the user has to
+ * fill appropriate data using a
+ * reinit of the blocks.
+ */
+ BlockVector (const unsigned int num_blocks);
+
+ /**
+ * Destructor. Clears memory
+ */
+ ~BlockVector ();
- for (unsigned int i=0; i<this->n_blocks(); ++i)
- this->components[i] = v.block(i);
+ /**
+ * Copy operator: fill all
+ * components of the vector that
+ * are locally stored with the
+ * given scalar value.
+ */
+ BlockVector &
+ operator = (const value_type s);
- return *this;
- }
+ /**
+ * Copy operator for a
+ * distributed Trilinos vector to
+ * a localized one.
+ */
+ BlockVector &
+ operator = (const MPI::BlockVector &V);
+ /**
+ * Copy operator for arguments of
+ * the same type.
+ */
+ BlockVector &
+ operator = (const BlockVector &V);
+ /**
+ * Another copy function. This
+ * one takes a deal.II block
+ * vector and copies it into a
+ * TrilinosWrappers block
+ * vector. Note that the number
+ * of blocks has to be the same
+ * in the vector as in the input
+ * vector. Use the reinit()
+ * command for resizing the
+ * BlockVector or for changing
+ * the internal structure of the
+ * block components.
+ *
+ * Since Trilinos only works on
+ * doubles, this function is
+ * limited to accept only one
+ * possible number type in the
+ * deal.II vector.
+ */
+ template <typename Number>
+ BlockVector &
+ operator = (const ::dealii::BlockVector<Number> &V);
- inline
- BlockVector::~BlockVector ()
- {}
+ /**
+ * Reinitialize the BlockVector to
+ * contain as many blocks as there
+ * are Epetra_Maps given in the input
+ * argument, according to the
+ * parallel distribution of the
+ * individual components described
+ * in the maps.
+ *
+ * If <tt>fast==false</tt>, the vector
+ * is filled with zeros.
+ */
+ void reinit (const std::vector<Epetra_Map> &input_maps,
+ const bool fast = false);
+ /**
+ * Reinit the function
+ * according to a distributed
+ * block vector. The elements
+ * will be copied in this
+ * process.
+ */
+ void reinit (const MPI::BlockVector &V);
+ /**
+ * Change the dimension to that
+ * of the vector <tt>V</tt>. The same
+ * applies as for the other
+ * reinit() function.
+ *
+ * The elements of <tt>V</tt> are not
+ * copied, i.e. this function is
+ * the same as calling <tt>reinit
+ * (V.size(), fast)</tt>.
+ *
+ * Note that you must call this
+ * (or the other reinit()
+ * functions) function, rather
+ * than calling the reinit()
+ * functions of an individual
+ * block, to allow the block
+ * vector to update its caches of
+ * vector sizes. If you call
+ * reinit() on one of the
+ * blocks, then subsequent
+ * actions on this object may
+ * yield unpredictable results
+ * since they may be routed to
+ * the wrong block.
+ */
+ void reinit (const BlockVector &V,
+ const bool fast = false);
- inline
- void
- BlockVector::reinit (const std::vector<Epetra_Map> &input_maps)
- {
- unsigned int no_blocks = input_maps.size();
- std::vector<unsigned int> block_sizes (no_blocks);
+ /**
+ * Change the number of blocks to
+ * <tt>num_blocks</tt>. The individual
+ * blocks will get initialized with
+ * zero size, so it is assumed that
+ * the user resizes the
+ * individual blocks by herself
+ * in an appropriate way, and
+ * calls <tt>collect_sizes</tt>
+ * afterwards.
+ */
+ void reinit (const unsigned int num_blocks);
- for (unsigned int i=0; i<no_blocks; ++i)
- {
- block_sizes[i] = input_maps[i].NumGlobalElements();
- }
+ /**
+ * Swap the contents of this
+ * vector and the other vector
+ * <tt>v</tt>. 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.
+ *
+ * Limitation: right now this
+ * function only works if both
+ * vectors have the same number
+ * of blocks. If needed, the
+ * numbers of blocks should be
+ * exchanged, too.
+ *
+ * This function is analog to the
+ * the swap() function of all C++
+ * standard containers. Also,
+ * there is a global function
+ * swap(u,v) that simply calls
+ * <tt>u.swap(v)</tt>, again in analogy
+ * to standard functions.
+ */
+ void swap (BlockVector &v);
- this->block_indices.reinit (block_sizes);
- if (this->components.size() != this->n_blocks())
- this->components.resize(this->n_blocks());
+ /**
+ * Exception
+ */
+ DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize);
- for (unsigned int i=0; i<this->n_blocks(); ++i)
- this->components[i].reinit(input_maps[i]);
+ /**
+ * Exception
+ */
+ DeclException0 (ExcNonMatchingBlockVectors);
- collect_sizes();
- }
+ /**
+ * Exception
+ */
+ DeclException2 (ExcNonLocalizedMap,
+ int, int,
+ << "For the generation of a localized vector the map has "
+ << "to assign all elements to all vectors! "
+ << "local_size = global_size is a necessary condition, but"
+ << arg1 << " != " << arg2 << " was given!");
+ };
- inline
- void
- BlockVector::reinit (const BlockVector& v,
- const bool fast)
- {
- this->block_indices = v.get_block_indices();
- if (this->components.size() != this->n_blocks())
- this->components.resize(this->n_blocks());
-
- for (unsigned int i=0;i<this->n_blocks();++i)
- block(i).reinit(v.block(i), fast);
-
- collect_sizes();
- }
+/*----------------------- Inline functions ----------------------------------*/
- inline
- void
- BlockVector::reinit (const unsigned int num_blocks)
- {
- std::vector<unsigned int> block_sizes (num_blocks, 0);
- this->block_indices.reinit (block_sizes);
- if (this->components.size() != this->n_blocks())
- this->components.resize(this->n_blocks());
-
- for (unsigned int i=0;i<this->n_blocks();++i)
- block(i).clear();
- }
+ inline
+ BlockVector::BlockVector ()
+ {}
- inline
- void
- BlockVector::compress ()
- {
- for (unsigned int i=0; i<n_blocks(); ++i)
- block(i).compress();
- }
+ inline
+ BlockVector::BlockVector (const std::vector<Epetra_Map> &InputMaps)
+ {
+ for (unsigned int i=0; i<InputMaps.size(); ++i)
+ {
+ Assert (InputMaps[i].NumGlobalElements() == InputMaps[i].NumMyElements(),
+ ExcNonLocalizedMap(InputMaps[i].NumGlobalElements(),
+ InputMaps[i].NumMyElements()));
+ }
+
+ reinit (InputMaps);
+ }
+
+
+
+ inline
+ BlockVector::BlockVector (const unsigned int num_blocks)
+ {
+ reinit (num_blocks);
+ }
+
+
+
+ inline
+ BlockVector::~BlockVector()
+ {}
+
+
+
+ inline
+ BlockVector::BlockVector (const MPI::BlockVector &v)
+ {
+ reinit (v);
+ }
+
+
+
+ inline
+ BlockVector::BlockVector (const BlockVector& v)
+ :
+ BlockVectorBase<Vector > ()
+ {
+ reinit(v);
+ }
- inline
- void
- BlockVector::swap (BlockVector &v)
- {
- Assert (this->n_blocks() == v.n_blocks(),
- ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
-
- for (unsigned int i=0; i<this->n_blocks(); ++i)
- this->components[i].swap (v.components[i]);
- ::dealii::swap (this->block_indices, v.block_indices);
- }
-
/**
*
* @ingroup TrilinosWrappers
* @relates TrilinosWrappers::BlockVector
- * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ * @author Martin Kronbichler, 2008
*/
inline
void swap (BlockVector &u,
namespace TrilinosWrappers
{
- class Vector;
+ class VectorBase;
class SparseMatrix;
/**
/**
* Apply the preconditioner.
*/
- void vmult (Vector &dst,
- const Vector &src) const;
+ void vmult (VectorBase &dst,
+ const VectorBase &src) const;
private:
boost::shared_ptr<Ifpack_Preconditioner> preconditioner;
#include <base/config.h>
#include <lac/exceptions.h>
-#include <lac/trilinos_vector.h>
+#include <lac/trilinos_vector_base.h>
#include <lac/trilinos_sparse_matrix.h>
#include <lac/vector.h>
#include <lac/sparse_matrix.h>
const bool higher_order_elements = false,
const double aggregation_threshold = 1e-4,
const std::vector<std::vector<bool> > &null_space = std::vector<std::vector<bool> > (),
- const bool output_details = false);
+ const bool output_details = false);
/**
* Let Trilinos compute a
* and return it in the
* destination vector.
*/
- void vmult (Vector &dst,
- const Vector &src) const;
+ void vmult (VectorBase &dst,
+ const VectorBase &src) const;
/**
* Do the same as before, but
#include <lac/sparsity_pattern.h>
#include <lac/sparse_matrix.h>
#include <lac/exceptions.h>
-#include <lac/trilinos_vector.h>
+#include <lac/trilinos_vector_base.h>
#include <boost/shared_ptr.hpp>
#include <vector>
namespace TrilinosWrappers
{
// forward declarations
- class Vector;
+ class VectorBase;
class SparseMatrix;
namespace MatrixIterators
* case you work on a distributed
* memory architecture, using the
* interface in the
- * TrilinosWrappers::Vector
- * class.
+ * TrilinosWrappers::VectorBase
+ * class (or one of the two
+ * derived classes Vector and
+ * MPI::Vector).
+ *
+ * In case of a localized Vector,
+ * this function will only work
+ * when running on one processor,
+ * since the matrix object is
+ * inherently
+ * distributed. Otherwise, and
+ * exception will be thrown.
*/
- void vmult (Vector &dst,
- const Vector &src) const;
+ void vmult (VectorBase &dst,
+ const VectorBase &src) const;
/**
* Matrix-vector multiplication:
* as was used for the matrix in
* case you work on a distributed
* memory architecture, using the
- * interface of the
- * TrilinosWrappers::Vector
- * class.
+ * interface in the
+ * TrilinosWrappers::VectorBase
+ * class (or one of the two
+ * derived classes Vector and
+ * MPI::Vector).
+ *
+ * In case of a localized Vector,
+ * this function will only work
+ * when running on one processor,
+ * since the matrix object is
+ * inherently
+ * distributed. Otherwise, and
+ * exception will be thrown.
*/
- void Tvmult (Vector &dst,
- const Vector &src) const;
+ void Tvmult (VectorBase &dst,
+ const VectorBase &src) const;
/**
* Adding Matrix-vector
* as was used for the matrix in
* case you work on a distributed
* memory architecture, using the
- * interface of the
- * TrilinosWrappers::Vector
- * class.
+ * interface in the
+ * TrilinosWrappers::VectorBase
+ * class (or one of the two
+ * derived classes Vector and
+ * MPI::Vector).
+ *
+ * In case of a localized Vector,
+ * this function will only work
+ * when running on one processor,
+ * since the matrix object is
+ * inherently
+ * distributed. Otherwise, and
+ * exception will be thrown.
*/
- void vmult_add (Vector &dst,
- const Vector &src) const;
+ void vmult_add (VectorBase &dst,
+ const VectorBase &src) const;
/**
* Adding Matrix-vector
* as was used for the matrix in
* case you work on a distributed
* memory architecture, using the
- * interface of the
- * TrilinosWrappers::Vector
- * class.
+ * interface in the
+ * TrilinosWrappers::VectorBase
+ * class (or one of the two
+ * derived classes Vector and
+ * MPI::Vector).
+ *
+ * In case of a localized Vector,
+ * this function will only work
+ * when running on one processor,
+ * since the matrix object is
+ * inherently
+ * distributed. Otherwise, and
+ * exception will be thrown.
*/
- void Tvmult_add (Vector &dst,
- const Vector &src) const;
+ void Tvmult_add (VectorBase &dst,
+ const VectorBase &src) const;
/**
* Return the square of the norm
* support this operation and
* needs a temporary vector.
*
- * Note that both vectors have
- * to be distributed vectors
+ * Note that both vectors have to
+ * be distributed vectors
* generated using the same Map
- * as was used for the matrix
- * in case you work on a
- * distributed memory
- * architecture, using the
- * interface of the
- * TrilinosWrappers::Vector
- * class.
+ * as was used for the matrix in
+ * case you work on a distributed
+ * memory architecture, using the
+ * interface in the
+ * TrilinosWrappers::VectorBase
+ * class (or one of the two
+ * derived classes Vector and
+ * MPI::Vector).
+ *
+ * In case of a localized Vector,
+ * this function will only work
+ * when running on one processor,
+ * since the matrix object is
+ * inherently
+ * distributed. Otherwise, and
+ * exception will be thrown.
*/
- TrilinosScalar matrix_norm_square (const Vector &v) const;
+ TrilinosScalar matrix_norm_square (const VectorBase &v) const;
/**
* Compute the matrix scalar
* this operation and needs a
* temporary vector.
*
- * Note that both vectors have
- * to be distributed vectors
+ * Note that both vectors have to
+ * be distributed vectors
* generated using the same Map
- * as was used for the matrix
- * in case you work on a
- * distributed memory
- * architecture, using the
- * interface of the
- * TrilinosWrappers::Vector
- * class.
+ * as was used for the matrix in
+ * case you work on a distributed
+ * memory architecture, using the
+ * interface in the
+ * TrilinosWrappers::VectorBase
+ * class (or one of the two
+ * derived classes Vector and
+ * MPI::Vector).
+ *
+ * In case of a localized Vector,
+ * this function will only work
+ * when running on one processor,
+ * since the matrix object is
+ * inherently
+ * distributed. Otherwise, and
+ * exception will be thrown.
*/
- TrilinosScalar matrix_scalar_product (const Vector &u,
- const Vector &v) const;
+ TrilinosScalar matrix_scalar_product (const VectorBase &u,
+ const VectorBase &v) const;
/**
* Compute the residual of an
* destination <i>dst</i> must
* not be the same vector.
*
- * Note that both vectors have
- * to be distributed vectors
+ * Note that both vectors have to
+ * be distributed vectors
* generated using the same Map
- * as was used for the matrix
- * in case you work on a
- * distributed memory
- * architecture, using the
- * interface of the
- * TrilinosWrappers::Vector
- * class.
+ * as was used for the matrix in
+ * case you work on a distributed
+ * memory architecture, using the
+ * interface in the
+ * TrilinosWrappers::VectorBase
+ * class (or one of the two
+ * derived classes Vector and
+ * MPI::Vector).
+ *
+ * In case of a localized Vector,
+ * this function will only work
+ * when running on one processor,
+ * since the matrix object is
+ * inherently
+ * distributed. Otherwise, and
+ * exception will be thrown.
*/
- TrilinosScalar residual (Vector &dst,
- const Vector &x,
- const Vector &b) const;
+ TrilinosScalar residual (VectorBase &dst,
+ const VectorBase &x,
+ const VectorBase &b) const;
/**
* Add <tt>matrix</tt> scaled
-#ifndef DOXYGEN
// -------------------------- inline and template functions ----------------------
+#ifndef DOXYGEN
+
namespace MatrixIterators
{
#ifndef __deal2__trilinos_vector_h
#define __deal2__trilinos_vector_h
-
#include <base/config.h>
#include <base/subscriptor.h>
#include <lac/exceptions.h>
#include <lac/vector.h>
-
-#include <vector>
-#include <utility>
-#include <memory>
+#include <lac/trilinos_vector_base.h>
#ifdef DEAL_II_USE_TRILINOS
-#define TrilinosScalar double
-# include "Epetra_ConfigDefs.h"
-# ifdef DEAL_II_COMPILER_SUPPORTS_MPI // only if MPI is installed
-# include "mpi.h"
-# include "Epetra_MpiComm.h"
-# else
-# include "Epetra_SerialComm.h"
-# endif
-# include "Epetra_FEVector.h"
# include "Epetra_Map.h"
# include "Epetra_LocalMap.h"
-# include "Epetra_MultiVector.h"
DEAL_II_NAMESPACE_OPEN
- // forward declaration
-template <typename number> class Vector;
+// forward declaration
+template <typename> class Vector;
/**
* @addtogroup TrilinosWrappers
*/
namespace TrilinosWrappers
{
- // forward declaration
- class Vector;
- class LocalizedVector;
-
- /**
- * @cond internal
- */
-
/**
- * A namespace for internal implementation details of the
- * TrilinosWrapper members. @ingroup TrilinosWrappers
+ * Namespace for Trilinos vector classes that work in parallel over
+ * MPI. This namespace is restricted to vectors only, whereas matrices
+ * are always MPI based when run on more than one processor.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
*/
- namespace internal
+ namespace MPI
{
- /**
- * This class implements a
- * wrapper for accessing the
- * Trilinos vector in the same
- * way as we access deal.II
- * objects: it is initialized
- * with a vector and an element
- * within it, and has a
- * conversion operator to
- * extract the scalar value of
- * this element. It also has a
- * variety of assignment
- * operator for writing to this
- * one element. @ingroup
- * TrilinosWrappers
- */
- class VectorReference
- {
- private:
- /**
- * Constructor. It is made
- * private so as to only allow
- * the actual vector class to
- * create it.
- */
- VectorReference (Vector &vector,
- const unsigned int index);
-
- public:
- /**
- * This looks like a copy
- * operator, but does something
- * different than usual. In
- * particular, it does not copy
- * the member variables of this
- * reference. Rather, it
- * handles the situation where
- * we have two vectors @p v and
- * @p w, and assign elements
- * like in
- * <tt>v(i)=w(i)</tt>. Here,
- * both left and right hand
- * side of the assignment have
- * data type VectorReference,
- * but what we really mean is
- * to assign the vector
- * elements represented by the
- * two references. This
- * operator implements this
- * operation. Note also that
- * this allows us to make the
- * assignment operator const.
- */
- const VectorReference & operator = (const VectorReference &r) const;
-
- /**
- * Set the referenced element of the
- * vector to <tt>s</tt>.
- */
- const VectorReference & operator = (const TrilinosScalar &s) const;
-
- /**
- * Add <tt>s</tt> to the
- * referenced element of the
- * vector->
- */
- const VectorReference & operator += (const TrilinosScalar &s) const;
-
- /**
- * Subtract <tt>s</tt> from the
- * referenced element of the
- * vector->
- */
- const VectorReference & operator -= (const TrilinosScalar &s) const;
-
- /**
- * Multiply the referenced
- * element of the vector by
- * <tt>s</tt>.
- */
- const VectorReference & operator *= (const TrilinosScalar &s) const;
-
- /**
- * Divide the referenced
- * element of the vector by
- * <tt>s</tt>.
- */
- const VectorReference & operator /= (const TrilinosScalar &s) const;
-
- /**
- * Convert the reference to an
- * actual value, i.e. return
- * the value of the referenced
- * element of the vector.
- */
- operator TrilinosScalar () const;
-
- /**
- * Exception
- */
- DeclException1 (ExcTrilinosError,
- int,
- << "An error with error number " << arg1
- << " occured while calling a Trilinos function");
-
- /**
- * Exception
- */
- DeclException3 (ExcAccessToNonLocalElement,
- int, int, int,
- << "You tried to access element " << arg1
- << " of a distributed vector, but only elements "
- << arg2 << " through " << arg3
- << " are stored locally and can be accessed.");
-
- private:
- /**
- * Point to the vector we are
- * referencing.
- */
- Vector &vector;
-
- /**
- * Index of the referenced element
- * of the vector.
- */
- const unsigned int index;
-
- /**
- * Make the vector class a
- * friend, so that it can
- * create objects of the
- * present type.
- */
- friend class ::dealii::TrilinosWrappers::Vector;
- };
- }
- /**
- * @endcond
- */
-
/**
* This class implements a wrapper to use the Trilinos distributed
- * vector class Epetra_FEVector. This is precisely the kind of vector we
- * deal with all the time - we probably get it from some assembly
- * process, where also entries not locally owned might need to written
- * and hence need to be forwarded to the owner. This class is designed
- * to be used in a distributed memory architecture with an MPI compiler
- * on the bottom, but works equally well also for serial processes. The
- * only requirement for this class to work is that Trilinos is installed
- * with the respective compiler as a basis.
- *
- * The interface of this class is modeled after the existing Vector
- * class in deal.II. It has almost the same member functions, and is
- * often exchangable. However, since Trilinos only supports a single
- * scalar type (double), it is not templated, and only works with that
- * type.
+ * vector class Epetra_FEVector. This class is derived from the
+ * TrilinosWrappers::VectorBase class and provides all functionality
+ * included there.
*
* Note that Trilinos only guarantees that operations do what you expect
* if the function @p GlobalAssemble has been called after vector
- * assembly in order to distribute the data. Therefore, you need to call
- * Vector::compress() before you actually use the vectors.
+ * assembly in order to distribute the data. This is necessary since
+ * some processes might have accumulated data of elements that are not
+ * owned by themselves, but must be sent to the owning process. In order
+ * to avoid using the wrong data, you need to call Vector::compress()
+ * before you actually use the vectors.
*
* <h3>Parallel communication model</h3>
*
* @see @ref SoftwareTrilinos
* @author Martin Kronbichler, Wolfgang Bangerth, 2008
*/
- class Vector : public Subscriptor
- {
- public:
- /**
- * Declare some of the standard
- * types used in all
- * containers. These types
- * parallel those in the
- * <tt>C</tt> standard libraries
- * <tt>vector<...></tt> class.
- */
- typedef TrilinosScalar value_type;
- typedef TrilinosScalar real_type;
- typedef size_t size_type;
- typedef internal::VectorReference reference;
- typedef const internal::VectorReference const_reference;
-
+ class Vector : public VectorBase
+ {
+ public:
/**
* Default constructor that
* generates an empty (zero size)
* processes in case of an MPI
* run.
*/
- Vector ();
+ Vector ();
/**
* This constructor takes an
* Epetra_Map that already knows
* how to distribute the
* individual components among
- * the MPI processors. It also
- * includes information about the
- * size of the vector.
+ * the MPI processors. Since it
+ * also includes information
+ * about the size of the vector,
+ * this is all we need to
+ * generate a parallel vector.
*/
- Vector (const Epetra_Map &InputMap);
+ Vector (const Epetra_Map &InputMap);
- /**
- * Copy constructor. Sets the
- * dimension to that of the given
- * vector and uses the map of
- * that vector, but does not copy
- * any element. Instead, the
- * memory will remain untouched
- * in case <tt>fast</tt> is false
- * and initialized with zero
- * otherwise.
- */
- Vector (const Vector &v,
- const bool fast = false);
+ /**
+ * Copy constructor using the
+ * given vector.
+ */
+ Vector (const Vector &V);
- /**
- * Destructor
- */
- virtual ~Vector ();
+ /**
+ * Copy constructor from the
+ * TrilinosWrappers vector
+ * class. Since a vector of this
+ * class does not necessarily
+ * need to be distributed among
+ * processes, the user needs to
+ * supply us with an Epetra_Map
+ * that sets the partitioning
+ * details.
+ */
+ explicit Vector (const Epetra_Map &InputMap,
+ const VectorBase &v);
+
+ /**
+ * Copy-constructor from deal.II
+ * vectors. Sets the dimension to that
+ * of the given vector, and copies all
+ * elements.
+ */
+ template <typename Number>
+ explicit Vector (const Epetra_Map &InputMap,
+ const dealii::Vector<Number> &v);
/**
* Reinit functionality. This
* new one based on the input
* map.
*/
- void reinit (const Epetra_Map &input_map);
+ void reinit (const Epetra_Map &input_map,
+ const bool fast = false);
/**
* Reinit functionality. This
- * function copies the vector v
- * to the current one.
+ * function sets the calling
+ * vector to the dimension and
+ * the parallel distribution of
+ * the input vector, but does not
+ * copy the elements in
+ * <tt>v</tt>. If <tt>fast</tt>
+ * is not <tt>true</tt>, the
+ * elements in the vector are
+ * initialized with zero,
+ * otherwise the content will be
+ * left unchanged and the user
+ * has to set all elements.
+ *
+ * This class has a third
+ * possible argument,
+ * <tt>allow_different_maps</tt>,
+ * that allows for an exchange of
+ * data between two equal-sized
+ * vectors (but being distributed
+ * differently among the
+ * processors). A trivial
+ * application of this function
+ * is to generate a replication
+ * of a whole vector on each
+ * machine, when the calling
+ * vector is build according to
+ * the localized vector class
+ * TrilinosWrappers::Vector, and
+ * <tt>v</tt> is a distributed
+ * vector. In this case, the
+ * variable <tt>fast</tt> needs
+ * to be set to <tt>false</tt>,
+ * since it does not make sense
+ * to exchange data between
+ * differently parallelized
+ * vectors without touching the
+ * elements.
*/
- void reinit (const Vector &v,
- const bool fast = false);
-
- /**
- * Release all memory and return
- * to a state just like after
- * having called the default
- * constructor.
- */
- void clear ();
-
- /**
- * Compress the underlying
- * representation of the Trilinos
- * object, i.e. flush the buffers
- * of the vector object if it has
- * any. This function is
- * necessary after writing into a
- * vector element-by-element and
- * before anything else can be
- * done on it.
- */
- void compress ();
-
- /**
- * Set all components of the
- * vector to the given number @p
- * s. Simply pass this down to
- * the Trilinos Epetra object,
- * but we still need to declare
- * this function to make the
- * example given in the
- * discussion about making the
- * constructor explicit work.
- *
- * Since the semantics of
- * assigning a scalar to a vector
- * are not immediately clear,
- * this operator should really
- * only be used if you want to
- * set the entire vector to
- * zero. This allows the
- * intuitive notation
- * <tt>v=0</tt>. Assigning other
- * values is deprecated and may
- * be disallowed in the future.
- */
- Vector &
- operator = (const TrilinosScalar s);
+ void reinit (const VectorBase &v,
+ const bool fast = false,
+ const bool allow_different_maps = false);
+
+ /**
+ * Set all components of the
+ * vector to the given number @p
+ * s. Simply pass this down to
+ * the base class, but we still
+ * need to declare this function
+ * to make the example given in
+ * the discussion about making
+ * the constructor explicit work.
+ */
+ Vector & operator = (const TrilinosScalar s);
/**
* Copy the given vector. Resize
* the present vector if
- * necessary.
+ * necessary. In this case, also
+ * the Epetra_Map that designs
+ * the parallel partitioning is
+ * taken from the input vector.
*/
- Vector &
- operator = (const Vector &V);
+ Vector &
+ operator = (const Vector &V);
/**
* Copy operator from a given
* side vector. Otherwise, an
* exception will be thrown.
*/
- Vector &
- operator = (const LocalizedVector &V);
+ Vector &
+ operator = (const ::dealii::TrilinosWrappers::Vector &V);
/**
* Another copy function. This
* possible number type in the
* deal.II vector.
*/
- Vector &
- operator = (const ::dealii::Vector<double> &v);
-
- /**
- * Test for equality. This
- * function assumes that the
- * present vector and the one to
- * compare with have the same
- * size already, since comparing
- * vectors of different sizes
- * makes not much sense anyway.
- */
- bool operator == (const Vector &v) const;
-
- /**
- * Test for inequality. This
- * function assumes that the
- * present vector and the one to
- * compare with have the same
- * size already, since comparing
- * vectors of different sizes
- * makes not much sense anyway.
- */
- bool operator != (const Vector &v) const;
-
- /**
- * Return the global dimension of
- * the vector.
- */
- unsigned int size () const;
-
- /**
- * Return the local dimension of
- * the vector, i.e. the number of
- * elements stored on the present
- * MPI process. For sequential
- * vectors, this number is the
- * same as size(), but for
- * parallel vectors it may be
- * smaller.
- *
- * To figure out which elements
- * exactly are stored locally,
- * use local_range().
- */
- unsigned int local_size () const;
-
- /**
- * Return a pair of indices
- * indicating which elements of
- * this vector are stored
- * locally. The first number is
- * the index of the first element
- * stored, the second the index
- * of the one past the last one
- * that is stored locally. If
- * this is a sequential vector,
- * then the result will be the
- * pair (0,N), otherwise it will
- * be a pair (i,i+n), where
- * <tt>n=local_size()</tt>.
- */
- std::pair<unsigned int, unsigned int> local_range () const;
-
- /**
- * Return whether @p index is in
- * the local range or not, see
- * also local_range().
- */
- bool in_local_range (const unsigned int index) const;
-
- /**
- * Provide access to a given
- * element, both read and write.
- */
- reference
- operator () (const unsigned int index);
-
- /**
- * Provide read-only access to an
- * element. This is equivalent to
- * * the <code>el()</code>
- * command.
- */
- TrilinosScalar
- operator () (const unsigned int index) const;
-
- /**
- * Return the value of the vector
- * entry <i>i</i>. Note that this
- * function does only work
- * properly when we request a
- * data stored on the local
- * processor. The function will
- * throw an exception in case the
- * elements sits on another
- * process.
- */
- TrilinosScalar el (const unsigned int index) const;
-
- /**
- * A collective set operation:
- * instead of setting individual
- * elements of a vector, this
- * function allows to set a whole
- * set of elements at once. The
- * indices of the elements to be
- * set are stated in the first
- * argument, the corresponding
- * values in the second.
- */
- void set (const std::vector<unsigned int> &indices,
- const std::vector<TrilinosScalar> &values);
-
- /**
- * This is a second collective
- * set operation. As a
- * difference, this function
- * takes a deal.II vector of
- * values.
- */
- void set (const std::vector<unsigned int> &indices,
- const ::dealii::Vector<TrilinosScalar> &values);
-
- /**
- * This collective set operation
- * is of lower level and can
- * handle anything else –
- * the only thing you have to
- * provide is an address where
- * all the indices are stored and
- * the number of elements to be
- * set.
- */
- void set (const unsigned int n_elements,
- const unsigned int *indices,
- const TrilinosScalar *values);
-
- /**
- * A collective add operation:
- * This funnction adds a whole
- * set of values stored in @p
- * values to the vector
- * components specified by @p
- * indices.
- */
- void add (const std::vector<unsigned int> &indices,
- const std::vector<TrilinosScalar> &values);
-
- /**
- * This is a second collective
- * add operation. As a
- * difference, this function
- * takes a deal.II vector of
- * values.
- */
- void add (const std::vector<unsigned int> &indices,
- const ::dealii::Vector<TrilinosScalar> &values);
-
- /**
- * Take an address where
- * <tt>n_elements</tt> are stored
- * contiguously and add them into
- * the vector.
- */
- void add (const unsigned int n_elements,
- const unsigned int *indices,
- const TrilinosScalar *values);
-
- /**
- * Return the scalar (inner)
- * product of two vectors. The
- * vectors must have the same
- * size.
- */
- TrilinosScalar operator * (const Vector &vec) const;
-
- /**
- * Return square of the
- * $l_2$-norm.
- */
- real_type norm_sqr () const;
-
- /**
- * Mean value of the elements of
- * this vector.
- */
- TrilinosScalar mean_value () const;
-
- /**
- * $l_1$-norm of the vector. The
- * sum of the absolute values.
- */
- real_type l1_norm () const;
-
- /**
- * $l_2$-norm of the vector. The
- * square root of the sum of the
- * squares of the elements.
- */
- real_type l2_norm () const;
-
- /**
- * $l_p$-norm of the vector. The
- * <i>p</i>th root of the sum of
- * the <i>p</i>th powers of the
- * absolute values of the
- * elements.
- */
- real_type lp_norm (const TrilinosScalar p) const;
-
- /**
- * Maximum absolute value of the
- * elements.
- */
- real_type linfty_norm () const;
-
- /**
- * Return whether the vector
- * contains only elements with
- * value zero. This function is
- * mainly for internal
- * consistency checks and should
- * seldomly be used when not in
- * debug mode since it uses quite
- * some time.
- */
- bool all_zero () const;
-
- /**
- * Return @p true if the vector
- * has no negative entries,
- * i.e. all entries are zero or
- * positive. This function is
- * used, for example, to check
- * whether refinement indicators
- * are really all positive (or
- * zero).
- */
- bool is_non_negative () const;
-
- /**
- * Multiply the entire vector by
- * a fixed factor.
- */
- Vector & operator *= (const TrilinosScalar factor);
-
- /**
- * Divide the entire vector by a
- * fixed factor.
- */
- Vector & operator /= (const TrilinosScalar factor);
-
- /**
- * Add the given vector to the
- * present one.
- */
- Vector & operator += (const Vector &V);
-
- /**
- * Subtract the given vector from
- * the present one.
- */
- Vector & operator -= (const Vector &V);
-
- /**
- * Addition of @p s to all
- * components. Note that @p s is
- * a scalar and not a vector.
- */
- void add (const TrilinosScalar s);
-
- /**
- * Simple vector addition, equal
- * to the <tt>operator +=</tt>.
- */
- void add (const Vector &V);
-
- /**
- * Simple addition of a multiple
- * of a vector, i.e. <tt>*this =
- * a*V</tt>.
- */
- void add (const TrilinosScalar a, const Vector &V);
-
- /**
- * Multiple addition of scaled
- * vectors, i.e. <tt>*this = a*V
- * + b*W</tt>.
- */
- void add (const TrilinosScalar a, const Vector &V,
- const TrilinosScalar b, const Vector &W);
-
- /**
- * Scaling and simple vector
- * addition, i.e. <tt>*this =
- * s*(*this) + V</tt>.
- */
- void sadd (const TrilinosScalar s,
- const Vector &V);
-
- /**
- * Scaling and simple addition,
- * i.e. <tt>*this = s*(*this) +
- * a*V</tt>.
- */
- void sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const Vector &V);
-
- /**
- * Scaling and multiple addition.
- */
- void sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const Vector &V,
- const TrilinosScalar b,
- const Vector &W);
+ template <typename Number>
+ Vector &
+ operator = (const ::dealii::Vector<Number> &v);
+ private:
/**
- * Scaling and multiple addition.
- * <tt>*this = s*(*this) + a*V +
- * b*W + c*X</tt>.
+ * The Epetra map is used to map
+ * (or rather, partition) vector
+ * data accross multiple
+ * processes. This is the
+ * communicator and data
+ * distribution object common to
+ * all Trilinos objects used by
+ * deal.II.
*/
- void sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const Vector &V,
- const TrilinosScalar b,
- const Vector &W,
- const TrilinosScalar c,
- const Vector &X);
+ Epetra_Map map;
+ };
- /**
- * Scale each element of this
- * vector by the corresponding
- * element in the argument. This
- * function is mostly meant to
- * simulate multiplication (and
- * immediate re-assignment) by a
- * diagonal scaling matrix.
- */
- void scale (const Vector &scaling_factors);
- /**
- * Assignment <tt>*this =
- * a*V</tt>.
- */
- void equ (const TrilinosScalar a, const Vector &V);
- /**
- * Assignment <tt>*this = a*V +
- * b*W</tt>.
- */
- void equ (const TrilinosScalar a, const Vector &V,
- const TrilinosScalar b, const Vector &W);
- /**
- * Compute the elementwise ratio
- * of the two given vectors, that
- * is let <tt>this[i] =
- * a[i]/b[i]</tt>. This is useful
- * for example if you want to
- * compute the cellwise ratio of
- * true to estimated error.
- *
- * This vector is appropriately
- * scaled to hold the result.
- *
- * If any of the <tt>b[i]</tt> is
- * zero, the result is
- * undefined. No attempt is made
- * to catch such situations.
- */
- void ratio (const Vector &a,
- const Vector &b);
+// ------------------- inline and template functions --------------
- /**
- * Output of vector in
- * user-defined format in analogy
- * to the dealii::Vector<number>
- * class.
- */
- void print (const char* format = 0) const;
- /**
- * Print to a stream. @p
- * precision denotes the desired
- * precision with which values
- * shall be printed, @p
- * scientific whether scientific
- * notation shall be used. If @p
- * across is @p true then the
- * vector is printed in a line,
- * while if @p false then the
- * elements are printed on a
- * separate line each.
- */
- void print (std::ostream &out,
- const unsigned int precision = 3,
- const bool scientific = true,
- const bool across = true) const;
+/**
+ * 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.
+ *
+ * @ingroup TrilinosWrappers
+ * @relates TrilinosWrappers::MPI::Vector
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ */
+ inline
+ void swap (Vector &u, Vector &v)
+ {
+ u.swap (v);
+ }
- /**
- * 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. Note that the vectors
- * need to be of the same size
- * and base on the same map.
- *
- * 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 (Vector &v);
- /**
- * Estimate for the memory
- * consumption (not implemented
- * for this class).
- */
- unsigned int memory_consumption () const;
+#ifndef DOXYGEN
- /**
- * Exception
- */
- DeclException1 (ExcTrilinosError,
- int,
- << "An error with error number " << arg1
- << " occured while calling a Trilinos function");
+ template <typename number>
+ Vector::Vector (const Epetra_Map &InputMap,
+ const dealii::Vector<number> &v)
+ :
+ map (InputMap)
+ {
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map, false));
+
+ const int min_my_id = map.MinMyGID();
+ const unsigned int size = map.NumMyElements();
- /**
- * Exception
- */
- DeclException3 (ExcAccessToNonlocalElement,
- int, int, int,
- << "You tried to access element " << arg1
- << " of a distributed vector, but only entries "
- << arg2 << " through " << arg3
- << " are stored locally and can be accessed.");
+ Assert (map.MaxLID() == size-1,
+ ExcDimensionMismatch(map.MaxLID(), size-1));
+ std::vector<int> indices (size);
+ for (unsigned int i=0; i<size; ++i)
+ indices[i] = map.GID(i);
- private:
- /**
- * The Epetra map is used to map
- * (or rather, partition) vector
- * data accross multiple
- * processes. This is the
- * communicator and data
- * distribution object common to
- * all Trilinos objects used by
- * deal.II.
- */
- Epetra_Map map;
- /**
- * Trilinos doesn't allow to
- * mix additions to matrix
- * entries and overwriting them
- * (to make synchronisation of
- * parallel computations
- * simpler). The way we do it
- * is to, for each access
- * operation, store whether it
- * is an insertion or an
- * addition. If the previous
- * one was of different type,
- * then we first have to flush
- * the Trilinos buffers;
- * otherwise, we can simply go
- * on. Luckily, Trilinos has
- * an object for this which
- * does already all the
- * parallel communications in
- * such a case, so we simply
- * use their model, which
- * stores whether the last
- * operation was an addition or
- * an insertion.
- */
- Epetra_CombineMode last_action;
+ const int ierr = vector->ReplaceGlobalValues(size, 0, v.begin()+min_my_id,
+ &indices[0]);
+ AssertThrow (ierr == 0, ExcTrilinosError());
+ }
- public:
- /**
- * An Epetra distibuted vector
- * type. Requires an existing
- * Epetra_Map for storing data.
- * TODO: Should become private
- * at some point.
- */
- std::auto_ptr<Epetra_FEVector> vector;
+
+
+ inline
+ Vector &
+ Vector::operator = (const TrilinosScalar s)
+ {
+ VectorBase::operator = (s);
+ return *this;
+ }
- /**
- * Make the reference class a
- * friend.
- */
- friend class internal::VectorReference;
- friend class LocalizedVector;
+#endif
- };
+ } /* end of namespace MPI */
/**
- * This class implements a localized version of a Trilinos vector. The
- * purpose of this class is to provide a copy interface from the
- * possibly parallel Vector class to a local vector on each processor.
+ * This class is a specialization of a Trilinos vector to a localized
+ * version. The purpose of this class is to provide a copy interface
+ * from the possibly parallel Vector class to a local vector on each
+ * processor, in order to be able to access all elements in the vector
+ * or to apply certain deal.II functions.
*
* @ingroup TrilinosWrappers
* @ingroup Vectors
* @see @ref SoftwareTrilinos
* @author Martin Kronbichler, 2008
*/
- class LocalizedVector : public Vector
+ class Vector : public VectorBase
{
public:
/**
* give the vector the correct
* size.
*/
- LocalizedVector ();
+ Vector ();
/**
* This constructor takes as
* input the number of elements
* in the vector.
*/
- LocalizedVector (const unsigned int n);
+ Vector (const unsigned int n);
+
+ /**
+ * This constructor takes as
+ * input the number of elements
+ * in the vector. If the map is
+ * not localized, i.e., if there
+ * are some elements that are not
+ * present on all processes, only
+ * the global size of the map
+ * will be taken and a localized
+ * map will be generated
+ * internally.
+ */
+ Vector (const Epetra_Map &InputMap);
/**
* This constructor takes a
* localized version of the whole
* content on each processor.
*/
- LocalizedVector (const Vector &V);
+ explicit Vector (const VectorBase &V);
- /**
- * This constructor generates a
- * LocalizedVector based on a
- * LocalizedVector input.
- */
- LocalizedVector (const LocalizedVector &V,
- const bool fast);
+ /**
+ * Copy-constructor from deal.II
+ * vectors. Sets the dimension to that
+ * of the given vector, and copies all
+ * elements.
+ */
+ template <typename Number>
+ explicit Vector (const dealii::Vector<Number> &v);
/**
* Reinit function that resizes
* the vector to the size
* specified by <tt>n</tt>.
*/
- void reinit (unsigned int n);
-
- /**
- * Reinit function. Takes the
- * information of a parallel
- * Trilinos vector and copies
- * everything to this local
- * version.
+ void reinit (const unsigned int n);
+
+ /**
+ * Initialization with an
+ * Epetra_Map. Similar to the
+ * call in the other class
+ * MPI::Vector, with the only
+ * difference that now a copy on
+ * all processes is
+ * generated. The variable
+ * <tt>fast</tt> determines
+ * whether the vector should be
+ * filled with zero or left
+ * untouched.
*/
- void reinit (const Vector &V);
+ void reinit (const Epetra_Map &input_map,
+ const bool fast = false);
/**
* Reinit function. Takes the
* information of a
- * LocalizedVector and copies
+ * Vector and copies
* everything to the calling
* vector.
*/
- void reinit (const LocalizedVector &V);
+ void reinit (const VectorBase &V,
+ const bool fast = false,
+ const bool allow_different_maps = false);
+
+ /**
+ * Set all components of the
+ * vector to the given number @p
+ * s. Simply pass this down to
+ * the base class, but we still
+ * need to declare this function
+ * to make the example given in
+ * the discussion about making
+ * the constructor explicit work.
+ */
+ Vector & operator = (const TrilinosScalar s);
/**
* Sets the left hand argument to
* Vector. Equivalent to the @p
* reinit function.
*/
- LocalizedVector &
- operator = (const Vector &V);
+ Vector &
+ operator = (const MPI::Vector &V);
+
+ /**
+ * Sets the left hand argument to
+ * the deal.II vector.
+ */
+ template <typename Number>
+ Vector &
+ operator = (const ::dealii::Vector<Number> &V);
/**
* Copy operator. Copies both the
* dimension and the content in
* the right hand argument.
*/
- LocalizedVector &
- operator = (const LocalizedVector &V);
+ Vector &
+ operator = (const Vector &V);
private:
Epetra_LocalMap map;
- public:
- std::auto_ptr<Epetra_MultiVector> vector;
-
- friend class Vector;
};
-
-
// ------------------- inline and template functions --------------
+
/**
- * Global function 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.
+ * 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.
*
* @ingroup TrilinosWrappers
* @relates TrilinosWrappers::Vector
* @author Martin Kronbichler, Wolfgang Bangerth, 2008
*/
- inline
- void swap (Vector &u, Vector &v)
- {
- u.swap (v);
- }
-
-#ifndef DOXYGEN
- namespace internal
- {
- inline
- VectorReference::VectorReference (Vector &vector,
- const unsigned int index)
- :
- vector (vector),
- index (index)
- {}
-
-
inline
- const VectorReference &
- VectorReference::operator = (const VectorReference &r) const
+ void swap (Vector &u, Vector &v)
{
- // as explained in the class
- // documentation, this is not the copy
- // operator. so simply pass on to the
- // "correct" assignment operator
- *this = static_cast<TrilinosScalar> (r);
-
- return *this;
+ u.swap (v);
}
+#ifndef DOXYGEN
- inline
- const VectorReference &
- VectorReference::operator = (const TrilinosScalar &value) const
+ template <typename number>
+ Vector::Vector (const dealii::Vector<number> &v)
+ :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ map (v.size(), 0, Epetra_MpiComm(MPI_COMM_WORLD))
+#else
+ map (v.size(), 0, Epetra_SerialComm())
+#endif
{
- vector.set (1, &index, &value);
- return *this;
- }
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map, false));
+ std::vector<int> indices (v.size());
+ for (unsigned int i=0; i<v.size(); ++i)
+ indices[i] = map.GID(i);
-
- inline
- const VectorReference &
- VectorReference::operator += (const TrilinosScalar &value) const
- {
- vector.add (1, &index, &value);
- return *this;
+ const int ierr = vector->ReplaceGlobalValues(v.size(), 0, v.begin(),
+ &indices[0]);
+ AssertThrow (ierr == 0, ExcTrilinosError());
}
-
-
+
+
inline
- const VectorReference &
- VectorReference::operator -= (const TrilinosScalar &value) const
+ Vector &
+ Vector::operator = (const TrilinosScalar s)
{
- TrilinosScalar new_value = -value;
- vector.add (1, &index, &new_value);
- return *this;
- }
-
-
+ VectorBase::operator = (s);
- inline
- const VectorReference &
- VectorReference::operator *= (const TrilinosScalar &value) const
- {
- TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
- vector.set (1, &index, &new_value);
return *this;
}
+#endif
- inline
- const VectorReference &
- VectorReference::operator /= (const TrilinosScalar &value) const
- {
- TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
- vector.set (1, &index, &new_value);
- return *this;
- }
- }
-
-
-
- inline
- bool
- Vector::in_local_range (const unsigned int index) const
- {
- std::pair<unsigned int, unsigned int> range = local_range();
-
- return ((index >= range.first) && (index < range.second));
- }
-
-
-
- inline
- Vector &
- Vector::operator = (const Vector &v)
- {
- // if the vectors have different sizes,
- // then first resize the present one
- if (map.SameAs(v.map) == false)
- {
- map = v.map;
- vector = std::auto_ptr<Epetra_FEVector>
- (new Epetra_FEVector(*v.vector));
- }
- else
- *vector = *v.vector;
-
- last_action = Insert;
-
- return *this;
- }
-
-
-
- inline
- internal::VectorReference
- Vector::operator () (const unsigned int index)
- {
- return internal::VectorReference (*this, index);
- }
-
-
-
- inline
- TrilinosScalar
- Vector::operator () (const unsigned int index) const
- {
- TrilinosScalar value = el(index);
-
- return value;
- }
-
-
-#endif // DOXYGEN
}
+
/*@}*/
DEAL_II_NAMESPACE_CLOSE
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id: trilinos_vector.h 16818 2008-09-12 14:30:24Z kronbichler $
+// Version: $Name$
+//
+// Copyright (C) 2008 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__trilinos_vector_base_h
+#define __deal2__trilinos_vector_base_h
+
+
+#include <base/config.h>
+#include <base/subscriptor.h>
+#include <lac/exceptions.h>
+#include <lac/vector.h>
+
+#include <vector>
+#include <utility>
+#include <memory>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+# define TrilinosScalar double
+# include "Epetra_ConfigDefs.h"
+# ifdef DEAL_II_COMPILER_SUPPORTS_MPI // only if MPI is installed
+# include "mpi.h"
+# include "Epetra_MpiComm.h"
+# else
+# include "Epetra_SerialComm.h"
+# endif
+# include "Epetra_FEVector.h"
+
+DEAL_II_NAMESPACE_OPEN
+
+ // forward declaration
+template <typename number> class Vector;
+
+
+/**
+ * @addtogroup TrilinosWrappers
+ *@{
+ */
+namespace TrilinosWrappers
+{
+ // forward declaration
+ class VectorBase;
+
+
+ /**
+ * @cond internal
+ */
+
+/**
+ * A namespace for internal implementation details of the
+ * TrilinosWrapper members.
+ *
+ * @ingroup TrilinosWrappers
+ */
+ namespace internal
+ {
+ /**
+ * This class implements a
+ * wrapper for accessing the
+ * Trilinos vector in the same
+ * way as we access deal.II
+ * objects: it is initialized
+ * with a vector and an element
+ * within it, and has a
+ * conversion operator to
+ * extract the scalar value of
+ * this element. It also has a
+ * variety of assignment
+ * operator for writing to this
+ * one element. @ingroup
+ * TrilinosWrappers
+ */
+ class VectorReference
+ {
+ private:
+ /**
+ * Constructor. It is made
+ * private so as to only allow
+ * the actual vector class to
+ * create it.
+ */
+ VectorReference (VectorBase &vector,
+ const unsigned int index);
+
+ public:
+ /**
+ * This looks like a copy
+ * operator, but does something
+ * different than usual. In
+ * particular, it does not copy
+ * the member variables of this
+ * reference. Rather, it
+ * handles the situation where
+ * we have two vectors @p v and
+ * @p w, and assign elements
+ * like in
+ * <tt>v(i)=w(i)</tt>. Here,
+ * both left and right hand
+ * side of the assignment have
+ * data type VectorReference,
+ * but what we really mean is
+ * to assign the vector
+ * elements represented by the
+ * two references. This
+ * operator implements this
+ * operation. Note also that
+ * this allows us to make the
+ * assignment operator const.
+ */
+ const VectorReference &
+ operator = (const VectorReference &r) const;
+
+ /**
+ * Set the referenced element of the
+ * vector to <tt>s</tt>.
+ */
+ const VectorReference &
+ operator = (const TrilinosScalar &s) const;
+
+ /**
+ * Add <tt>s</tt> to the
+ * referenced element of the
+ * vector->
+ */
+ const VectorReference &
+ operator += (const TrilinosScalar &s) const;
+
+ /**
+ * Subtract <tt>s</tt> from the
+ * referenced element of the
+ * vector->
+ */
+ const VectorReference &
+ operator -= (const TrilinosScalar &s) const;
+
+ /**
+ * Multiply the referenced
+ * element of the vector by
+ * <tt>s</tt>.
+ */
+ const VectorReference &
+ operator *= (const TrilinosScalar &s) const;
+
+ /**
+ * Divide the referenced
+ * element of the vector by
+ * <tt>s</tt>.
+ */
+ const VectorReference &
+ operator /= (const TrilinosScalar &s) const;
+
+ /**
+ * Convert the reference to an
+ * actual value, i.e. return
+ * the value of the referenced
+ * element of the vector.
+ */
+ operator TrilinosScalar () const;
+
+ /**
+ * Exception
+ */
+ DeclException1 (ExcTrilinosError,
+ int,
+ << "An error with error number " << arg1
+ << " occured while calling a Trilinos function");
+
+ /**
+ * Exception
+ */
+ DeclException3 (ExcAccessToNonLocalElement,
+ int, int, int,
+ << "You tried to access element " << arg1
+ << " of a distributed vector, but only elements "
+ << arg2 << " through " << arg3
+ << " are stored locally and can be accessed.");
+
+ private:
+ /**
+ * Point to the vector we are
+ * referencing.
+ */
+ VectorBase &vector;
+
+ /**
+ * Index of the referenced element
+ * of the vector.
+ */
+ const unsigned int index;
+
+ /**
+ * Make the vector class a
+ * friend, so that it can
+ * create objects of the
+ * present type.
+ */
+ friend class ::dealii::TrilinosWrappers::VectorBase;
+ };
+ }
+ /**
+ * @endcond
+ */
+
+
+/**
+ * Base class for the two types of Trilinos vectors, the distributed
+ * memory vector MPI::Vector and a localized vector Vector. The latter
+ * is designed for use in either serial implementations or as a
+ * localized copy on each processor. The implementation of this class
+ * is based on the Trilinos vector class Epetra_FEVector, the (parallel)
+ * partitioning of which is governed by an Epetra_Map. This means that
+ * the vector type is generic and can be done in this base class, while
+ * the definition of the partition map (and hence, the constructor and
+ * reinit function) will have to be done in the derived classes. The
+ * Epetra_FEVector is precisely the kind of vector we deal with all the
+ * time - we probably get it from some assembly process, where also
+ * entries not locally owned might need to written and hence need to be
+ * forwarded to the owner. The only requirement for this class to work
+ * is that Trilinos is installed with the same compiler as is used for
+ * compilation of deal.II.
+ *
+ * The interface of this class is modeled after the existing Vector
+ * class in deal.II. It has almost the same member functions, and is
+ * often exchangable. However, since Trilinos only supports a single
+ * scalar type (double), it is not templated, and only works with that
+ * type.
+ *
+ * Note that Trilinos only guarantees that operations do what you expect
+ * if the function @p GlobalAssemble has been called after vector
+ * assembly in order to distribute the data. Therefore, you need to call
+ * Vector::compress() before you actually use the vectors.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Vectors
+ * @see @ref SoftwareTrilinos
+ * @author Martin Kronbichler, 2008
+ */
+ class VectorBase : public Subscriptor
+ {
+ public:
+ /**
+ * Declare some of the standard
+ * types used in all
+ * containers. These types
+ * parallel those in the
+ * <tt>C</tt> standard libraries
+ * <tt>vector<...></tt> class.
+ */
+ typedef TrilinosScalar value_type;
+ typedef TrilinosScalar real_type;
+ typedef size_t size_type;
+ typedef internal::VectorReference reference;
+ typedef const internal::VectorReference const_reference;
+
+ /**
+ * Default constructor that
+ * generates an empty (zero size)
+ * vector. The function
+ * <tt>reinit()</tt> will have to
+ * give the vector the correct
+ * size and distribution among
+ * processes in case of an MPI
+ * run.
+ */
+ VectorBase ();
+
+ /**
+ * Copy constructor. Sets the
+ * dimension to that of the given
+ * vector, and copies all the
+ * elements.
+ */
+ VectorBase (const VectorBase &v);
+
+ /**
+ * Destructor
+ */
+ virtual ~VectorBase ();
+
+ /**
+ * Release all memory and return
+ * to a state just like after
+ * having called the default
+ * constructor.
+ */
+ void clear ();
+
+ /**
+ * Compress the underlying
+ * representation of the Trilinos
+ * object, i.e. flush the buffers
+ * of the vector object if it has
+ * any. This function is
+ * necessary after writing into a
+ * vector element-by-element and
+ * before anything else can be
+ * done on it.
+ */
+ void compress ();
+
+ /**
+ * Returns the state of the
+ * vector, i.e., whether
+ * compress() has already been
+ * called after an operation
+ * requiring data exchange. Does
+ * only return non-false values
+ * when used in <tt>debug</tt>
+ * mode.
+ */
+ bool is_compressed () const;
+
+ /**
+ * Set all components of the
+ * vector to the given number @p
+ * s. Simply pass this down to
+ * the Trilinos Epetra object,
+ * but we still need to declare
+ * this function to make the
+ * example given in the
+ * discussion about making the
+ * constructor explicit work.
+ *
+ * Since the semantics of
+ * assigning a scalar to a vector
+ * are not immediately clear,
+ * this operator should really
+ * only be used if you want to
+ * set the entire vector to
+ * zero. This allows the
+ * intuitive notation
+ * <tt>v=0</tt>. Assigning other
+ * values is deprecated and may
+ * be disallowed in the future.
+ */
+ VectorBase &
+ operator = (const TrilinosScalar s);
+
+ /**
+ * Test for equality. This
+ * function assumes that the
+ * present vector and the one to
+ * compare with have the same
+ * size already, since comparing
+ * vectors of different sizes
+ * makes not much sense anyway.
+ */
+ bool operator == (const VectorBase &v) const;
+
+ /**
+ * Test for inequality. This
+ * function assumes that the
+ * present vector and the one to
+ * compare with have the same
+ * size already, since comparing
+ * vectors of different sizes
+ * makes not much sense anyway.
+ */
+ bool operator != (const VectorBase &v) const;
+
+ /**
+ * Return the global dimension of
+ * the vector.
+ */
+ unsigned int size () const;
+
+ /**
+ * Return the local dimension of
+ * the vector, i.e. the number of
+ * elements stored on the present
+ * MPI process. For sequential
+ * vectors, this number is the
+ * same as size(), but for
+ * parallel vectors it may be
+ * smaller.
+ *
+ * To figure out which elements
+ * exactly are stored locally,
+ * use local_range().
+ */
+ unsigned int local_size () const;
+
+ /**
+ * Return a pair of indices
+ * indicating which elements of
+ * this vector are stored
+ * locally. The first number is
+ * the index of the first element
+ * stored, the second the index
+ * of the one past the last one
+ * that is stored locally. If
+ * this is a sequential vector,
+ * then the result will be the
+ * pair (0,N), otherwise it will
+ * be a pair (i,i+n), where
+ * <tt>n=local_size()</tt>.
+ */
+ std::pair<unsigned int, unsigned int> local_range () const;
+
+ /**
+ * Return whether @p index is in
+ * the local range or not, see
+ * also local_range().
+ */
+ bool in_local_range (const unsigned int index) const;
+
+ /**
+ * Provide access to a given
+ * element, both read and write.
+ */
+ reference
+ operator () (const unsigned int index);
+
+ /**
+ * Provide read-only access to an
+ * element. This is equivalent to
+ * the <code>el()</code> command.
+ */
+ TrilinosScalar
+ operator () (const unsigned int index) const;
+
+ /**
+ * Return the value of the vector
+ * entry <i>i</i>. Note that this
+ * function does only work
+ * properly when we request a
+ * data stored on the local
+ * processor. The function will
+ * throw an exception in case the
+ * elements sits on another
+ * process.
+ */
+ TrilinosScalar el (const unsigned int index) const;
+
+ /**
+ * A collective set operation:
+ * instead of setting individual
+ * elements of a vector, this
+ * function allows to set a whole
+ * set of elements at once. The
+ * indices of the elements to be
+ * set are stated in the first
+ * argument, the corresponding
+ * values in the second.
+ */
+ void set (const std::vector<unsigned int> &indices,
+ const std::vector<TrilinosScalar> &values);
+
+ /**
+ * This is a second collective
+ * set operation. As a
+ * difference, this function
+ * takes a deal.II vector of
+ * values.
+ */
+ void set (const std::vector<unsigned int> &indices,
+ const ::dealii::Vector<TrilinosScalar> &values);
+
+ /**
+ * This collective set operation
+ * is of lower level and can
+ * handle anything else —
+ * the only thing you have to
+ * provide is an address where
+ * all the indices are stored and
+ * the number of elements to be
+ * set.
+ */
+ void set (const unsigned int n_elements,
+ const unsigned int *indices,
+ const TrilinosScalar *values);
+
+ /**
+ * A collective add operation:
+ * This funnction adds a whole
+ * set of values stored in @p
+ * values to the vector
+ * components specified by @p
+ * indices.
+ */
+ void add (const std::vector<unsigned int> &indices,
+ const std::vector<TrilinosScalar> &values);
+
+ /**
+ * This is a second collective
+ * add operation. As a
+ * difference, this function
+ * takes a deal.II vector of
+ * values.
+ */
+ void add (const std::vector<unsigned int> &indices,
+ const ::dealii::Vector<TrilinosScalar> &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.
+ */
+ void add (const unsigned int n_elements,
+ const unsigned int *indices,
+ const TrilinosScalar *values);
+
+ /**
+ * Return the scalar (inner)
+ * product of two vectors. The
+ * vectors must have the same
+ * size.
+ */
+ TrilinosScalar operator * (const VectorBase &vec) const;
+
+ /**
+ * Return square of the
+ * $l_2$-norm.
+ */
+ real_type norm_sqr () const;
+
+ /**
+ * Mean value of the elements of
+ * this vector.
+ */
+ TrilinosScalar mean_value () const;
+
+ /**
+ * $l_1$-norm of the vector. The
+ * sum of the absolute values.
+ */
+ real_type l1_norm () const;
+
+ /**
+ * $l_2$-norm of the vector. The
+ * square root of the sum of the
+ * squares of the elements.
+ */
+ real_type l2_norm () const;
+
+ /**
+ * $l_p$-norm of the vector. The
+ * <i>p</i>th root of the sum of
+ * the <i>p</i>th powers of the
+ * absolute values of the
+ * elements.
+ */
+ real_type lp_norm (const TrilinosScalar p) const;
+
+ /**
+ * Maximum absolute value of the
+ * elements.
+ */
+ real_type linfty_norm () const;
+
+ /**
+ * Return whether the vector
+ * contains only elements with
+ * value zero. This function is
+ * mainly for internal
+ * consistency checks and should
+ * seldomly be used when not in
+ * debug mode since it uses quite
+ * some time.
+ */
+ bool all_zero () const;
+
+ /**
+ * Return @p true if the vector
+ * has no negative entries,
+ * i.e. all entries are zero or
+ * positive. This function is
+ * used, for example, to check
+ * whether refinement indicators
+ * are really all positive (or
+ * zero).
+ */
+ bool is_non_negative () const;
+
+ /**
+ * Multiply the entire vector by
+ * a fixed factor.
+ */
+ VectorBase & operator *= (const TrilinosScalar factor);
+
+ /**
+ * Divide the entire vector by a
+ * fixed factor.
+ */
+ VectorBase & operator /= (const TrilinosScalar factor);
+
+ /**
+ * Add the given vector to the
+ * present one.
+ */
+ VectorBase & operator += (const VectorBase &V);
+
+ /**
+ * Subtract the given vector from
+ * the present one.
+ */
+ VectorBase & operator -= (const VectorBase &V);
+
+ /**
+ * Addition of @p s to all
+ * components. Note that @p s is
+ * a scalar and not a vector.
+ */
+ void add (const TrilinosScalar s);
+
+ /**
+ * Simple vector addition, equal
+ * to the <tt>operator
+ * +=</tt>.
+ *
+ * Though, if the second argument
+ * <tt>allow_different_maps</tt>
+ * is set, then it is possible to
+ * add data from a different map.
+ */
+ void add (const VectorBase &V,
+ const bool allow_different_maps = false);
+
+ /**
+ * Simple addition of a multiple
+ * of a vector, i.e. <tt>*this =
+ * a*V</tt>.
+ */
+ void add (const TrilinosScalar a,
+ const VectorBase &V);
+
+ /**
+ * Multiple addition of scaled
+ * vectors, i.e. <tt>*this = a*V
+ * + b*W</tt>.
+ */
+ void add (const TrilinosScalar a,
+ const VectorBase &V,
+ const TrilinosScalar b,
+ const VectorBase &W);
+
+ /**
+ * Scaling and simple vector
+ * addition, i.e. <tt>*this =
+ * s*(*this) + V</tt>.
+ */
+ void sadd (const TrilinosScalar s,
+ const VectorBase &V);
+
+ /**
+ * Scaling and simple addition,
+ * i.e. <tt>*this = s*(*this) +
+ * a*V</tt>.
+ */
+ void sadd (const TrilinosScalar s,
+ const TrilinosScalar a,
+ const VectorBase &V);
+
+ /**
+ * Scaling and multiple addition.
+ */
+ void sadd (const TrilinosScalar s,
+ const TrilinosScalar a,
+ const VectorBase &V,
+ const TrilinosScalar b,
+ const VectorBase &W);
+
+ /**
+ * Scaling and multiple addition.
+ * <tt>*this = s*(*this) + a*V +
+ * b*W + c*X</tt>.
+ */
+ void sadd (const TrilinosScalar s,
+ const TrilinosScalar a,
+ const VectorBase &V,
+ const TrilinosScalar b,
+ const VectorBase &W,
+ const TrilinosScalar c,
+ const VectorBase &X);
+
+ /**
+ * Scale each element of this
+ * vector by the corresponding
+ * element in the argument. This
+ * function is mostly meant to
+ * simulate multiplication (and
+ * immediate re-assignment) by a
+ * diagonal scaling matrix.
+ */
+ void scale (const VectorBase &scaling_factors);
+
+ /**
+ * Assignment <tt>*this =
+ * a*V</tt>.
+ */
+ void equ (const TrilinosScalar a,
+ const VectorBase &V);
+
+ /**
+ * Assignment <tt>*this = a*V +
+ * b*W</tt>.
+ */
+ void equ (const TrilinosScalar a,
+ const VectorBase &V,
+ const TrilinosScalar b,
+ const VectorBase &W);
+
+ /**
+ * Compute the elementwise ratio
+ * of the two given vectors, that
+ * is let <tt>this[i] =
+ * a[i]/b[i]</tt>. This is useful
+ * for example if you want to
+ * compute the cellwise ratio of
+ * true to estimated error.
+ *
+ * This vector is appropriately
+ * scaled to hold the result.
+ *
+ * If any of the <tt>b[i]</tt> is
+ * zero, the result is
+ * undefined. No attempt is made
+ * to catch such situations.
+ */
+ void ratio (const VectorBase &a,
+ const VectorBase &b);
+
+ /**
+ * Output of vector in
+ * user-defined format in analogy
+ * to the dealii::Vector<number>
+ * class.
+ */
+ void print (const char* format = 0) const;
+
+ /**
+ * Print to a stream. @p
+ * precision denotes the desired
+ * precision with which values
+ * shall be printed, @p
+ * scientific whether scientific
+ * notation shall be used. If @p
+ * across is @p true then the
+ * vector is printed in a line,
+ * while if @p false then the
+ * elements are printed on a
+ * separate line each.
+ */
+ void print (std::ostream &out,
+ const unsigned int precision = 3,
+ const bool scientific = true,
+ const bool across = true) const;
+
+ /**
+ * 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. Note that the vectors
+ * need to be of the same size
+ * and base on the same map.
+ *
+ * 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 (VectorBase &v);
+
+ /**
+ * Estimate for the memory
+ * consumption (not implemented
+ * for this class).
+ */
+ unsigned int memory_consumption () const;
+
+ /**
+ * Exception
+ */
+ DeclException1 (ExcTrilinosError,
+ int,
+ << "An error with error number " << arg1
+ << " occured while calling a Trilinos function");
+
+ /**
+ * Exception
+ */
+ DeclException3 (ExcAccessToNonlocalElement,
+ int, int, int,
+ << "You tried to access element " << arg1
+ << " of a distributed vector, but only entries "
+ << arg2 << " through " << arg3
+ << " are stored locally and can be accessed.");
+
+
+ private:
+ /**
+ * Trilinos doesn't allow to
+ * mix additions to matrix
+ * entries and overwriting them
+ * (to make synchronisation of
+ * parallel computations
+ * simpler). The way we do it
+ * is to, for each access
+ * operation, store whether it
+ * is an insertion or an
+ * addition. If the previous
+ * one was of different type,
+ * then we first have to flush
+ * the Trilinos buffers;
+ * otherwise, we can simply go
+ * on. Luckily, Trilinos has
+ * an object for this which
+ * does already all the
+ * parallel communications in
+ * such a case, so we simply
+ * use their model, which
+ * stores whether the last
+ * operation was an addition or
+ * an insertion.
+ */
+ Epetra_CombineMode last_action;
+
+ /**
+ * A boolean variable to hold
+ * information on whether the
+ * vector is compressed or not.
+ */
+ bool compressed;
+
+ public:
+ /**
+ * An Epetra distibuted vector
+ * type. Requires an existing
+ * Epetra_Map for storing data.
+ * TODO: Should become private
+ * at some point.
+ */
+ std::auto_ptr<Epetra_FEVector> vector;
+
+
+ /**
+ * Make the reference class a
+ * friend.
+ */
+ friend class internal::VectorReference;
+ friend class Vector;
+ friend class MPI::Vector;
+ };
+
+
+
+
+// ------------------- inline and template functions --------------
+
+/**
+ * Global function 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.
+ *
+ * @ingroup TrilinosWrappers
+ * @relates TrilinosWrappers::VectorBase
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ */
+ inline
+ void swap (VectorBase &u, VectorBase &v)
+ {
+ u.swap (v);
+ }
+
+
+#ifndef DOXYGEN
+
+ namespace internal
+ {
+ inline
+ VectorReference::VectorReference (VectorBase &vector,
+ const unsigned int index)
+ :
+ vector (vector),
+ index (index)
+ {}
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator = (const VectorReference &r) const
+ {
+ // as explained in the class
+ // documentation, this is not the copy
+ // operator. so simply pass on to the
+ // "correct" assignment operator
+ *this = static_cast<TrilinosScalar> (r);
+
+ return *this;
+ }
+
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator = (const TrilinosScalar &value) const
+ {
+ vector.set (1, &index, &value);
+ return *this;
+ }
+
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator += (const TrilinosScalar &value) const
+ {
+ vector.add (1, &index, &value);
+ return *this;
+ }
+
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator -= (const TrilinosScalar &value) const
+ {
+ TrilinosScalar new_value = -value;
+ vector.add (1, &index, &new_value);
+ return *this;
+ }
+
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator *= (const TrilinosScalar &value) const
+ {
+ TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) * value;
+ vector.set (1, &index, &new_value);
+ return *this;
+ }
+
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator /= (const TrilinosScalar &value) const
+ {
+ TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
+ vector.set (1, &index, &new_value);
+ return *this;
+ }
+ }
+
+
+
+ inline
+ bool
+ VectorBase::is_compressed () const
+ {
+ return compressed;
+ }
+
+
+
+ inline
+ bool
+ VectorBase::in_local_range (const unsigned int index) const
+ {
+ std::pair<unsigned int, unsigned int> range = local_range();
+
+ return ((index >= range.first) && (index < range.second));
+ }
+
+
+
+ inline
+ internal::VectorReference
+ VectorBase::operator () (const unsigned int index)
+ {
+ return internal::VectorReference (*this, index);
+ }
+
+
+#endif // DOXYGEN
+
+}
+
+/*@}*/
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS
+
+/*---------------------------- trilinos_vector_base.h ---------------------------*/
+
+#endif
+/*---------------------------- trilinos_vector_base.h ---------------------------*/
#ifdef DEAL_II_USE_TRILINOS
namespace TrilinosWrappers
{
+ namespace MPI
+ {
+ class Vector;
+ }
class Vector;
}
#endif
#ifdef DEAL_II_USE_TRILINOS
/**
- * Another copy constructor: copy the
- * values from a Trilinos wrapper
- * vector class. This copy constructor is
- * only available if Trilinos was detected
- * during configuration time.
+ * Another copy constructor: copy
+ * the values from a Trilinos
+ * wrapper vector. This copy
+ * constructor is only available if
+ * Trilinos was detected during
+ * configuration time.
*
- * Note that due to the communication
- * model used in MPI, this operation can
- * only succeed if all processes do it at
- * the same time. I.e., it is not
- * possible for only one process to
- * obtain a copy of a parallel vector
- * while the other jobs do something
- * else.
+ * Note that due to the
+ * communication model used in MPI,
+ * this operation can only succeed
+ * if all processes do it at the
+ * same time. This means that it is
+ * not possible for only one
+ * process to obtain a copy of a
+ * parallel vector while the other
+ * jobs do something else. This
+ * call will rather result in a
+ * copy of the vector on all
+ * processors.
+ */
+ Vector (const TrilinosWrappers::MPI::Vector &v);
+
+ /**
+ * Another copy constructor: copy
+ * the values from a localized
+ * Trilinos wrapper vector. This
+ * copy constructor is only
+ * available if Trilinos was
+ * detected during configuration
+ * time.
*/
Vector (const TrilinosWrappers::Vector &v);
#endif
#ifdef DEAL_II_USE_PETSC
/**
- * Another copy operator: copy the values
- * from a sequential PETSc wrapper vector
- * class. This operator is only available
- * if PETSc was detected during
+ * Another copy operator: copy the
+ * values from a sequential PETSc
+ * wrapper vector class. This
+ * operator is only available if
+ * PETSc was detected during
* configuration time.
*/
Vector<Number> &
operator = (const PETScWrappers::Vector &v);
/**
- * Another copy operator: copy the values
- * from a parallel PETSc wrapper vector
- * class. This operator is only available
- * if PETSc was detected during
+ * Another copy operator: copy the
+ * values from a parallel PETSc
+ * wrapper vector class. This
+ * operator is only available if
+ * PETSc was detected during
* configuration time.
*
- * Note that due to the communication
- * model used in MPI, this operation can
- * only succeed if all processes do it at
- * the same time. I.e., it is not
+ * Note that due to the
+ * communication model used in MPI,
+ * this operation can only succeed
+ * if all processes do it at the
+ * same time. I.e., it is not
* possible for only one process to
- * obtain a copy of a parallel vector
- * while the other jobs do something
- * else.
+ * obtain a copy of a parallel
+ * vector while the other jobs do
+ * something else.
*/
Vector<Number> &
operator = (const PETScWrappers::MPI::Vector &v);
* Trilinos was detected during
* configuration time.
*
- * Note that due to the communication
- * model used in MPI, this operation can
- * only succeed if all processes do it at
- * the same time. I.e., it is not
+ * Note that due to the
+ * communication model used in MPI,
+ * this operation can only succeed
+ * if all processes do it at the
+ * same time. I.e., it is not
* possible for only one process to
- * obtain a copy of a parallel vector
- * while the other jobs do something
- * else.
+ * obtain a copy of a parallel
+ * vector while the other jobs do
+ * something else.
+ */
+ Vector<Number> &
+ operator = (const TrilinosWrappers::MPI::Vector &v);
+
+ /**
+ * Another copy operator: copy the
+ * values from a sequential
+ * Trilinos wrapper vector
+ * class. This operator is only
+ * available if Trilinos was
+ * detected during configuration
+ * time.
*/
Vector<Number> &
operator = (const TrilinosWrappers::Vector &v);
/**
* Test for equality. This function
- * assumes that the present vector and
- * the one to compare with have the same
- * size already, since comparing vectors
- * of different sizes makes not much
- * sense anyway.
+ * assumes that the present vector
+ * and the one to compare with have
+ * the same size already, since
+ * comparing vectors of different
+ * sizes makes not much sense
+ * anyway.
*/
template <typename Number2>
bool operator == (const Vector<Number2> &v) const;
#ifdef DEAL_II_USE_TRILINOS
template <typename Number>
-Vector<Number>::Vector (const TrilinosWrappers::Vector &v)
+Vector<Number>::Vector (const TrilinosWrappers::MPI::Vector &v)
:
Subscriptor(),
vec_size(v.size()),
// be a better solution than
// this, but it has not yet been
// found.
- TrilinosWrappers::LocalizedVector localized_vector (v);
+ TrilinosWrappers::Vector localized_vector (v);
// get a representation of the vector
// and copy it
}
}
+
+
+template <typename Number>
+Vector<Number>::Vector (const TrilinosWrappers::Vector &v)
+ :
+ Subscriptor(),
+ vec_size(v.size()),
+ max_vec_size(v.size()),
+ val(0)
+{
+ if (vec_size != 0)
+ {
+ val = new Number[max_vec_size];
+ Assert (val != 0, ExcOutOfMemory());
+
+ // get a representation of the vector
+ // and copy it
+ TrilinosScalar **start_ptr;
+
+ int ierr = v.vector->ExtractView (&start_ptr);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ std::copy (start_ptr[0], start_ptr[0]+vec_size, begin());
+ }
+}
+
#endif
template <typename Number>
#ifdef DEAL_II_USE_TRILINOS
+template <typename Number>
+Vector<Number> &
+Vector<Number>::operator = (const TrilinosWrappers::MPI::Vector &v)
+{
+ // Generate a localized version
+ // of the Trilinos vectors and
+ // then call the other =
+ // operator.
+ TrilinosWrappers::Vector localized_vector (v);
+ *this = localized_vector;
+ return *this;
+}
+
+
+
template <typename Number>
Vector<Number> &
Vector<Number>::operator = (const TrilinosWrappers::Vector &v)
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id: trilinos_block_vector.cc 15631 2008-01-17 23:47:31Z bangerth $
+// Version: $Name$
+//
+// Copyright (C) 2008 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <lac/trilinos_block_vector.h>
+
+
+#ifdef DEAL_II_USE_TRILINOS
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+ namespace MPI
+ {
+ BlockVector &
+ BlockVector::operator = (const value_type s)
+ {
+ BaseClass::operator = (s);
+ return *this;
+ }
+
+
+
+ BlockVector &
+ BlockVector::operator = (const BlockVector &v)
+ {
+ if (this->n_blocks() != v.n_blocks())
+ reinit(v.n_blocks());
+
+ for (unsigned int i=0; i<this->n_blocks(); ++i)
+ this->components[i] = v.block(i);
+
+ collect_sizes();
+
+ return *this;
+ }
+
+
+
+ template <typename Number>
+ BlockVector &
+ BlockVector::operator = (const ::dealii::BlockVector<Number> &v)
+ {
+ Assert (n_blocks() == v.n_blocks(),
+ ExcDimensionMismatch(n_blocks(),v.n_blocks()));
+
+ for (unsigned int i=0; i<this->n_blocks(); ++i)
+ this->components[i] = v.block(i);
+
+ return *this;
+ }
+
+
+
+ BlockVector::~BlockVector ()
+ {}
+
+
+
+ void
+ BlockVector::reinit (const std::vector<Epetra_Map> &input_maps,
+ const bool fast)
+ {
+ unsigned int no_blocks = input_maps.size();
+ std::vector<unsigned int> block_sizes (no_blocks);
+
+ for (unsigned int i=0; i<no_blocks; ++i)
+ {
+ block_sizes[i] = input_maps[i].NumGlobalElements();
+ }
+
+ this->block_indices.reinit (block_sizes);
+ if (components.size() != n_blocks())
+ components.resize(n_blocks());
+
+ for (unsigned int i=0; i<n_blocks(); ++i)
+ components[i].reinit(input_maps[i], fast);
+
+ collect_sizes();
+ }
+
+
+
+ void
+ BlockVector::reinit (const BlockVector& v,
+ const bool fast)
+ {
+ block_indices = v.get_block_indices();
+ if (components.size() != n_blocks())
+ components.resize(n_blocks());
+
+ for (unsigned int i=0;i<n_blocks();++i)
+ components[i].reinit(v.block(i), fast);
+
+ collect_sizes();
+ }
+
+
+
+ void
+ BlockVector::reinit (const unsigned int num_blocks)
+ {
+ std::vector<unsigned int> block_sizes (num_blocks, 0);
+ this->block_indices.reinit (block_sizes);
+ if (this->components.size() != this->n_blocks())
+ this->components.resize(this->n_blocks());
+
+ for (unsigned int i=0;i<this->n_blocks();++i)
+ block(i).clear();
+ }
+
+
+
+ void
+ BlockVector::compress ()
+ {
+ for (unsigned int i=0; i<n_blocks(); ++i)
+ components[i].compress();
+ }
+
+
+
+ void
+ BlockVector::swap (BlockVector &v)
+ {
+ Assert (this->n_blocks() == v.n_blocks(),
+ ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
+
+ for (unsigned int i=0; i<this->n_blocks(); ++i)
+ this->components[i].swap (v.components[i]);
+ ::dealii::swap (this->block_indices, v.block_indices);
+ }
+
+ } /* end of namespace MPI */
+
+
+
+
+
+
+ BlockVector &
+ BlockVector::operator = (const value_type s)
+ {
+ BaseClass::operator = (s);
+ return *this;
+ }
+
+
+
+ void
+ BlockVector::reinit (const std::vector<Epetra_Map> &input_maps,
+ const bool fast)
+ {
+ unsigned int no_blocks = input_maps.size();
+ std::vector<unsigned int> block_sizes (no_blocks);
+
+ for (unsigned int i=0; i<no_blocks; ++i)
+ block_sizes[i] = input_maps[i].NumGlobalElements();
+
+
+ this->block_indices.reinit (block_sizes);
+ if (components.size() != n_blocks())
+ components.resize(n_blocks());
+
+ for (unsigned int i=0; i<n_blocks(); ++i)
+ components[i].reinit(input_maps[i], fast);
+
+ collect_sizes();
+ }
+
+
+
+ void
+ BlockVector::reinit (const MPI::BlockVector &v)
+ {
+ block_indices = v.get_block_indices();
+ if (components.size() != n_blocks())
+ components.resize(n_blocks());
+
+ for (unsigned int i=0;i<n_blocks();++i)
+ components[i].reinit(v.block(i));
+
+ collect_sizes();
+ }
+
+
+
+ void
+ BlockVector::reinit (const unsigned int num_blocks)
+ {
+ std::vector<unsigned int> block_sizes (num_blocks, 0);
+ this->block_indices.reinit (block_sizes);
+ if (this->components.size() != this->n_blocks())
+ this->components.resize(this->n_blocks());
+
+ for (unsigned int i=0;i<this->n_blocks();++i)
+ block(i).clear();
+ }
+
+
+
+ void
+ BlockVector::reinit (const BlockVector &v,
+ const bool fast)
+ {
+ block_indices = v.get_block_indices();
+ if (components.size() != n_blocks())
+ components.resize(n_blocks());
+
+ for (unsigned int i=0;i<n_blocks();++i)
+ components[i].reinit(v.block(i), fast);
+
+ collect_sizes();
+ }
+
+
+
+ BlockVector &
+ BlockVector::operator = (const MPI::BlockVector &v)
+ {
+ reinit (v);
+
+ return *this;
+ }
+
+
+
+ BlockVector &
+ BlockVector::operator = (const BlockVector &v)
+ {
+ Assert (n_blocks() == v.n_blocks(),
+ ExcDimensionMismatch(n_blocks(),v.n_blocks()));
+
+ for (unsigned int i=0; i<this->n_blocks(); ++i)
+ this->components[i] = v.block(i);
+
+ return *this;
+ }
+
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
#include <lac/trilinos_precondition.h>
-#include <lac/trilinos_vector.h>
+#include <lac/trilinos_vector_base.h>
#include <lac/trilinos_sparse_matrix.h>
#ifdef DEAL_II_USE_TRILINOS
void
- PreconditionSSOR::vmult (Vector &dst,
- const Vector &src) const
+ PreconditionSSOR::vmult (VectorBase &dst,
+ const VectorBase &src) const
{
preconditioner->ApplyInverse (*src.vector, *dst.vector);
}
- void PreconditionAMG::vmult (Vector &dst,
- const Vector &src) const
+ void PreconditionAMG::vmult (VectorBase &dst,
+ const VectorBase &src) const
{
const int ierr = multigrid_operator->ApplyInverse (*src.vector,
*dst.vector);
//
//---------------------------------------------------------------------------
-#include <lac/trilinos_vector.h>
+#include <lac/trilinos_vector_base.h>
#include <lac/trilinos_sparse_matrix.h>
#ifdef DEAL_II_USE_TRILINOS
ierr = matrix->GlobalAssemble(col_map, row_map, false);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- last_action = Insert;
}
+ if (last_action != Insert)
+ last_action = Insert;
+
int trilinos_i = i;
int trilinos_j = j;
ierr = matrix->GlobalAssemble(col_map, row_map, false);
AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- last_action = Add;
}
+ if (last_action != Add)
+ last_action = Add;
+
// we have to do above actions in any
// case to be consistent with the MPI
// communication model (see the
}
+
void
- SparseMatrix::vmult (Vector &dst,
- const Vector &src) const
+ SparseMatrix::vmult (VectorBase &dst,
+ const VectorBase &src) const
{
Assert (&src != &dst, ExcSourceEqualsDestination());
void
- SparseMatrix::Tvmult (Vector &dst,
- const Vector &src) const
+ SparseMatrix::Tvmult (VectorBase &dst,
+ const VectorBase &src) const
{
Assert (&src != &dst, ExcSourceEqualsDestination());
void
- SparseMatrix::vmult_add (Vector &dst,
- const Vector &src) const
+ SparseMatrix::vmult_add (VectorBase &dst,
+ const VectorBase &src) const
{
Assert (&src != &dst, ExcSourceEqualsDestination());
- Vector tmp;
- tmp = dst;
+ VectorBase tmp(dst);
vmult (dst, src);
dst += tmp;
}
void
- SparseMatrix::Tvmult_add (Vector &dst,
- const Vector &src) const
+ SparseMatrix::Tvmult_add (VectorBase &dst,
+ const VectorBase &src) const
{
Assert (&src != &dst, ExcSourceEqualsDestination());
- Vector tmp;
- tmp = dst;
+ VectorBase tmp(dst);
vmult (dst, src);
dst += tmp;
}
TrilinosScalar
- SparseMatrix::matrix_norm_square (const Vector &v) const
+ SparseMatrix::matrix_norm_square (const VectorBase &v) const
{
- Vector tmp(v,true);
+ VectorBase tmp (v);
vmult (tmp, v);
return tmp*v;
}
TrilinosScalar
- SparseMatrix::matrix_scalar_product (const Vector &u,
- const Vector &v) const
+ SparseMatrix::matrix_scalar_product (const VectorBase &u,
+ const VectorBase &v) const
{
- Vector tmp(v,true);
+ VectorBase tmp (v);
vmult (tmp, v);
return u*tmp;
}
TrilinosScalar
- SparseMatrix::residual (Vector &dst,
- const Vector &x,
- const Vector &b) const
+ SparseMatrix::residual (VectorBase &dst,
+ const VectorBase &x,
+ const VectorBase &b) const
{
vmult (dst, x);
dst -= b;
namespace TrilinosWrappers
{
- namespace internal
+ namespace MPI
{
- VectorReference::operator TrilinosScalar () const
- {
- Assert (index < vector.size(),
- ExcIndexRange (index, 0, vector.size()));
-
- // Trilinos allows for vectors
- // to be referenced by the []
- // or () operators but only ()
- // checks index bounds Also,
- // can only get local values
-
- AssertThrow ((static_cast<signed int>(index) >= vector.map.MinMyGID()) &&
- (static_cast<signed int>(index) <= vector.map.MaxMyGID()),
- ExcAccessToNonLocalElement (index, vector.map.MinMyGID(),
- vector.map.MaxMyGID()));
-
- return (*(vector.vector))[0][index];
- }
- }
-
- Vector::Vector ()
- :
+ Vector::Vector ()
+ :
#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
+ map (0,0,Epetra_MpiComm(MPI_COMM_WORLD))
#else
- map (0,0,Epetra_SerialComm()),
+ map (0,0,Epetra_SerialComm())
#endif
- last_action (Insert),
- vector(std::auto_ptr<Epetra_FEVector>
- (new Epetra_FEVector(map)))
- {}
+ {
+ last_action = Zero;
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map));
+ }
- Vector::Vector (const Epetra_Map &InputMap)
- :
- map (InputMap),
- last_action (Insert),
- vector (std::auto_ptr<Epetra_FEVector>
- (new Epetra_FEVector(map)))
- {}
+ Vector::Vector (const Epetra_Map &InputMap)
+ :
+ map (InputMap)
+ {
+ last_action = Zero;
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map));
+ }
- Vector::Vector (const Vector &v,
- const bool fast)
- :
- Subscriptor(),
- map (v.map),
- last_action (Insert),
- vector(std::auto_ptr<Epetra_FEVector>
- (new Epetra_FEVector(map,!fast)))
- {}
-
-
-
- Vector::~Vector ()
- {
- }
-
-
-
- void
- Vector::reinit (const Epetra_Map &input_map)
- {
- vector.reset();
- map = input_map;
-
- vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map));
- last_action = Insert;
- }
-
-
-
- void
- Vector::reinit (const Vector &v,
- const bool fast)
- {
- vector.reset();
- if (map.SameAs(v.map) == false)
- map = v.map;
-
- vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map,!fast));
- last_action = Insert;
- }
-
-
+ Vector::Vector (const Vector &v)
+ :
+ VectorBase(),
+ map (v.map)
+ {
+ last_action = Zero;
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(*v.vector));
+ }
- void
- Vector::clear ()
- {
- // When we clear the vector,
- // reset the pointer and generate
- // an empty vector.
- vector.reset();
-#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- map = Epetra_Map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD));
-#else
- map = Epetra_Map (0, 0, Epetra_SerialComm());
-#endif
- vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map));
- last_action = Insert;
- }
+ Vector::Vector (const Epetra_Map &InputMap,
+ const VectorBase &v)
+ :
+ VectorBase(),
+ map (InputMap)
+ {
+ AssertThrow (map.NumGlobalElements() == v.vector->Map().NumGlobalElements(),
+ ExcDimensionMismatch (map.NumGlobalElements(),
+ v.vector->Map().NumGlobalElements()));
+
+ last_action = Zero;
+
+ if (map.SameAs(v.vector->Map()) == true)
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(*v.vector));
+ else
+ {
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map, false));
+ reinit (v, false, true);
+ }
+ }
- void
- Vector::compress ()
- {
- // Now pass over the information
- // about what we did last to Trilinos.
- const int ierr = vector->GlobalAssemble(last_action);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
+ void
+ Vector::reinit (const Epetra_Map &input_map,
+ const bool fast)
+ {
+ vector.reset();
+ map = input_map;
+
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map,!fast));
+ last_action = Zero;
+ }
- Vector &
- Vector::operator = (const TrilinosScalar s)
- {
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but either "
- "infinite or Not A Number (NaN)"));
+ void
+ Vector::reinit (const VectorBase &v,
+ const bool fast,
+ const bool allow_different_maps)
+ {
+ // In case we do not allow to
+ // have different maps, this
+ // call means that we have to
+ // reset the vector. So clear
+ // the vector, initialize our
+ // map with the map in v, and
+ // generate the vector.
+ if (allow_different_maps == false)
+ {
+ vector.reset();
+ if (map.SameAs(v.vector->Map()) == false)
+ map = Epetra_Map(v.vector->GlobalLength(),0,v.vector->Comm());
+
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map,!fast));
+ last_action = Zero;
+ }
+
+ // Otherwise, we have to check
+ // that the two vectors are
+ // already of the same size,
+ // create an object for the data
+ // exchange and then insert all
+ // the data.
+ else
+ {
+ Assert (fast == false,
+ ExcMessage ("It is not possible to exchange data with the "
+ "option fast set, which would not write "
+ "elements."));
+
+ AssertThrow (size() == v.size(),
+ ExcDimensionMismatch (size(), v.size()));
+
+ Epetra_Import data_exchange (vector->Map(), v.vector->Map());
+
+ const int ierr = vector->Import(*v.vector, data_exchange, Insert);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
- const int ierr = vector->PutScalar(s);
+ }
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- return *this;
- }
+ Vector &
+ Vector::operator = (const Vector &v)
+ {
+ if (vector->Map().SameAs(v.vector->Map()) == true)
+ *vector = *v.vector;
+ else
+ {
+ vector.reset();
+ map = v.map;
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(*v.vector));
+ }
+
+ return *this;
+ }
- Vector &
- Vector::operator = (const LocalizedVector &v)
- {
- Assert (size() == v.size(), ExcDimensionMismatch(size(), v.size()));
- Epetra_Import import_information (map, v.map);
- const int ierr = vector->Import(*v.vector, import_information, Insert);
+ Vector &
+ Vector::operator = (const TrilinosWrappers::Vector &v)
+ {
+ Assert (size() == v.size(), ExcDimensionMismatch(size(), v.size()));
- Assert (ierr == 0, ExcTrilinosError(ierr));
+ Epetra_Import data_exchange (vector->Map(), v.vector->Map());
+ const int ierr = vector->Import(*v.vector, data_exchange, Insert);
- return *this;
- }
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ return *this;
+ }
- Vector &
- Vector::operator = (const ::dealii::Vector<double> &v)
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(),v.size()));
- const unsigned int n_local_elements = local_size();
- const unsigned int first_local_element = local_range().first;
+ template <typename Number>
+ Vector &
+ Vector::operator = (const ::dealii::Vector<Number> &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
- std::vector<int> local_indices (n_local_elements);
- for (unsigned int i=0; i<n_local_elements; ++i)
- local_indices[i] = map.GID(i);
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map, false));
+
+ const int min_my_id = map.MinMyGID();
+ const unsigned int size = map.NumMyElements();
- const int ierr = vector->ReplaceGlobalValues(n_local_elements,
- &local_indices[0],
- v.begin()+first_local_element);
+ Assert (map.MaxLID() == size-1,
+ ExcDimensionMismatch(map.MaxLID(), size-1));
- Assert (ierr==0, ExcTrilinosError(ierr));
+ std::vector<int> indices (size);
+ for (unsigned int i=0; i<size; ++i)
+ indices[i] = map.GID(i);
- return *this;
- }
+ const int ierr = vector->ReplaceGlobalValues(size, 0, v.begin()+min_my_id,
+ &indices[0]);
+ AssertThrow (ierr == 0, ExcTrilinosError());
+ return *this;
+ }
- bool
- Vector::operator == (const Vector &v) const
- {
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
- unsigned int i;
- for (i=0; i<size(); i++)
- if ((*(v.vector))[0][i]!=(*vector)[0][i]) return false;
+ } /* end of namespace MPI */
- return true;
- }
- bool
- Vector::operator != (const Vector &v) const
+ Vector::Vector ()
+ :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD))
+#else
+ map (0, 0, Epetra_SerialComm())
+#endif
{
- Assert (size() == v.size(),
- ExcDimensionMismatch(size(), v.size()));
-
- return (!(*this==v));
+ last_action = Zero;
+ vector = std::auto_ptr<Epetra_FEVector>
+ (new Epetra_FEVector(map,false));
}
- unsigned int
- Vector::size () const
+ Vector::Vector (const unsigned int n)
+ :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ map ((int)n, 0, Epetra_MpiComm(MPI_COMM_WORLD))
+#else
+ map ((int)n, 0, Epetra_SerialComm())
+#endif
{
- return (unsigned int) vector->Map().NumGlobalElements();
+ reinit (n);
}
- unsigned int
- Vector::local_size () const
+ Vector::Vector (const Epetra_Map &InputMap)
+ :
+ map (InputMap.NumGlobalElements(), 0, InputMap.Comm())
{
- return (unsigned int) vector->Map().NumMyElements();
+ last_action = Zero;
+ vector = std::auto_ptr<Epetra_FEVector>
+ (new Epetra_FEVector(map,true));
}
- std::pair<unsigned int, unsigned int>
- Vector::local_range () const
+ Vector::Vector (const VectorBase &v)
+ :
+ map (v.vector->Map().NumGlobalElements(), 0, v.vector->Comm())
{
- int begin, end;
- begin = vector->Map().MinMyGID();
- end = vector->Map().MaxMyGID()+1;
- return std::make_pair (begin, end);
- }
-
+ last_action = Zero;
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map,false));
-
- TrilinosScalar
- Vector::el (const unsigned int index) const
- {
- // Extract local indices in
- // the vector.
- int trilinos_i = vector->Map().LID(index);
- TrilinosScalar value = 0.;
-
- // If the element is not
- // present on the current
- // processor, we can't
- // continue. Just print out 0.
-
- // TODO: Is this reasonable?
- if (trilinos_i == -1 )
- {
- return 0.;
- //Assert (false, ExcAccessToNonlocalElement(index, local_range().first,
- // local_range().second));
- }
+ if (vector->Map().SameAs(v.vector->Map()) == true)
+ *vector = *v.vector;
else
- value = (*vector)[0][trilinos_i];
-
- return value;
- }
+ reinit (v, false, true);
-
-
- void
- Vector::set (const std::vector<unsigned int> &indices,
- const std::vector<TrilinosScalar> &values)
- {
- Assert (indices.size() == values.size(),
- ExcDimensionMismatch(indices.size(),values.size()));
-
- set (indices.size(), &indices[0], &values[0]);
}
void
- Vector::set (const std::vector<unsigned int> &indices,
- const ::dealii::Vector<TrilinosScalar> &values)
+ Vector::reinit (const unsigned int n)
{
- Assert (indices.size() == values.size(),
- ExcDimensionMismatch(indices.size(),values.size()));
-
- set (indices.size(), &indices[0], values.begin());
- }
-
-
+ vector.reset();
- void
- Vector::set (const unsigned int n_elements,
- const unsigned int *indices,
- const TrilinosScalar *values)
- {
- if (last_action == Add)
+ if (map.NumGlobalElements() != (int)n)
{
- vector->GlobalAssemble(Add);
- last_action = Insert;
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ map = Epetra_LocalMap ((int)n, 0, Epetra_MpiComm(MPI_COMM_WORLD));
+#else
+ map = Epetra_LocalMap ((int)n, 0, Epetra_SerialComm());
+#endif
}
- const int ierr= vector->ReplaceGlobalValues (n_elements,
- (int*)(const_cast<unsigned int*>(indices)),
- const_cast<TrilinosScalar*>(values));
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
+ last_action = Zero;
-
-
- void
- Vector::add (const std::vector<unsigned int> &indices,
- const std::vector<TrilinosScalar> &values)
- {
- Assert (indices.size() == values.size(),
- ExcDimensionMismatch(indices.size(),values.size()));
-
- add (indices.size(), &indices[0], &values[0]);
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector (map,true));
}
void
- Vector::add (const std::vector<unsigned int> &indices,
- const ::dealii::Vector<TrilinosScalar> &values)
+ Vector::reinit (const Epetra_Map &input_map,
+ const bool fast)
{
- Assert (indices.size() == values.size(),
- ExcDimensionMismatch(indices.size(),values.size()));
-
- add (indices.size(), &indices[0], values.begin());
- }
-
-
+ vector.reset();
- void
- Vector::add (const unsigned int n_elements,
- const unsigned int *indices,
- const TrilinosScalar *values)
- {
- if (last_action == Insert)
+ if (map.NumGlobalElements() != input_map.NumGlobalElements())
{
- vector->GlobalAssemble(Insert);
- last_action = Add;
+ map = Epetra_LocalMap (input_map.NumGlobalElements(),0,
+ input_map.Comm());
}
- const int ierr= vector->SumIntoGlobalValues (n_elements,
- (int*)(const_cast<unsigned int*>(indices)),
- const_cast<TrilinosScalar*>(values));
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- TrilinosScalar
- Vector::operator * (const Vector &vec) const
- {
- Assert (size() == vec.size(),
- ExcDimensionMismatch(size(), vec.size()));
-
- TrilinosScalar result;
-
- const int ierr = vector->Dot(*(vec.vector), &result);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return result;
- }
-
-
-
- Vector::real_type
- Vector::norm_sqr () const
- {
- const TrilinosScalar d = l2_norm();
- return d*d;
- }
-
-
-
- TrilinosScalar
- Vector::mean_value () const
- {
- TrilinosScalar mean;
-
- const int ierr = vector->MeanValue (&mean);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return mean;
- }
-
-
-
- Vector::real_type
- Vector::l1_norm () const
- {
- TrilinosScalar d;
-
- const int ierr = vector->Norm1 (&d);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ last_action = Zero;
- return d;
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector (map, !fast));
}
- Vector::real_type
- Vector::l2_norm () const
- {
- TrilinosScalar d;
-
- const int ierr = vector->Norm2 (&d);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return d;
- }
-
-
-
- Vector::real_type
- Vector::lp_norm (const TrilinosScalar p) const
- {
- // get a representation of the
- // vector and loop over all
- // the elements
- TrilinosScalar *start_ptr;
- int leading_dimension;
- int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- TrilinosScalar norm = 0;
- TrilinosScalar sum=0;
-
- const TrilinosScalar * ptr = start_ptr;
-
- // add up elements
- while (ptr != start_ptr+size())
- sum += std::pow(std::fabs(*ptr++), p);
-
- norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
-
- return norm;
- }
-
-
-
- Vector::real_type
- Vector::linfty_norm () const
- {
- TrilinosScalar d;
-
- const int ierr = vector->NormInf (&d);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return d;
- }
-
-
-
- bool
- Vector::all_zero () const
- {
- // get a representation of the vector and
- // loop over all the elements
- TrilinosScalar *start_ptr;
- int leading_dimension;
- int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- const TrilinosScalar *ptr = start_ptr,
- *eptr = start_ptr + size();
- bool flag = true;
- while (ptr != eptr)
+ void
+ Vector::reinit (const VectorBase &v,
+ const bool fast,
+ const bool allow_different_maps)
+ {
+ // In case we do not allow to
+ // have different maps, this
+ // call means that we have to
+ // reset the vector. So clear
+ // the vector, initialize our
+ // map with the map in v, and
+ // generate the vector.
+ if (allow_different_maps == false)
{
- if (*ptr != 0)
- {
- flag = false;
- break;
- }
- ++ptr;
- }
-
- return flag;
- }
+ vector.reset();
+ if (map.SameAs(v.vector->Map()) == false)
+ map = Epetra_LocalMap (v.vector->GlobalLength(),0,v.vector->Comm());
-
-
- bool
- Vector::is_non_negative () const
- {
- // get a representation of the vector and
- // loop over all the elements
- TrilinosScalar *start_ptr;
- int leading_dimension;
- int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- const TrilinosScalar *ptr = start_ptr,
- *eptr = start_ptr + size();
- bool flag = true;
- while (ptr != eptr)
- {
- if (*ptr < 0.0)
- {
- flag = false;
- break;
- }
- ++ptr;
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map,!fast));
+ last_action = Zero;
}
- return flag;
- }
-
-
+ // Otherwise, we have to check
+ // that the two vectors are
+ // already of the same size,
+ // create an object for the data
+ // exchange and then insert all
+ // the data.
+ else
+ {
+ Assert (fast == false,
+ ExcMessage ("It is not possible to exchange data with the "
+ "option fast set, which would not write "
+ "elements."));
- Vector &
- Vector::operator *= (const TrilinosScalar a)
- {
+ AssertThrow (size() == v.size(),
+ ExcDimensionMismatch (size(), v.size()));
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
+ Epetra_Import data_exchange (vector->Map(), v.vector->Map());
- const int ierr = vector->Scale(a);
- Assert (ierr == 0, ExcTrilinosError(ierr));
+ const int ierr = vector->Import(*v.vector, data_exchange, Insert);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
- return *this;
}
Vector &
- Vector::operator /= (const TrilinosScalar a)
+ Vector::operator = (const MPI::Vector &v)
{
+ if (vector->Map().SameAs(v.vector->Map()) == false)
+ {
+ map = Epetra_LocalMap (v.vector->Map().NumGlobalElements(), 0,
+ v.vector->Comm());
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map,false));
+ }
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const TrilinosScalar factor = 1./a;
-
- Assert (numbers::is_finite(factor),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Scale(factor);
- Assert (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- Vector &
- Vector::operator += (const Vector &v)
- {
- const int ierr = vector->Update (1.0, *(v.vector), 1.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ reinit (v, false, true);
return *this;
}
Vector &
- Vector::operator -= (const Vector &v)
+ Vector::operator = (const Vector &v)
{
- const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- return *this;
- }
-
-
-
- void
- Vector::add (const TrilinosScalar s)
- {
-
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- unsigned int n_local = local_size();
- int ierr;
-
- for (unsigned int i=0; i<n_local; i++)
+ if (vector->Map().SameAs(v.vector->Map()) == false)
{
- ierr = vector->SumIntoMyValue(i,0,s);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ map = Epetra_LocalMap (v.vector->Map().NumGlobalElements(), 0,
+ v.vector->Comm());
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map,false));
}
- }
-
-
-
- void
- Vector::add (const Vector &v)
- {
- *this += v;
- }
-
-
-
- void
- Vector::add (const TrilinosScalar a,
- const Vector &v)
- {
-
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Update(a, *(v.vector), 1.);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- Vector::add (const TrilinosScalar a,
- const Vector &v,
- const TrilinosScalar b,
- const Vector &w)
- {
-
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(b),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- Vector::sadd (const TrilinosScalar s,
- const Vector &v)
- {
-
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Update(1., *(v.vector), s);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- Vector::sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const Vector &v)
- {
-
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Update(a, *(v.vector), s);
-
- AssertThrow(ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- Vector::sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const Vector &v,
- const TrilinosScalar b,
- const Vector &w)
- {
-
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(b),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- Vector::sadd (const TrilinosScalar s,
- const TrilinosScalar a,
- const Vector &v,
- const TrilinosScalar b,
- const Vector &w,
- const TrilinosScalar c,
- const Vector &x)
- {
-
- Assert (numbers::is_finite(s),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(b),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(c),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- // Update member can only
- // input two other vectors so
- // do it in two steps
- const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
-
- const int jerr = vector->Update(c, *(x.vector), 1.);
- AssertThrow (jerr == 0, ExcTrilinosError(jerr));
- }
-
-
-
- void
- Vector::scale (const Vector &factors)
- {
- const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- Vector::equ (const TrilinosScalar a,
- const Vector &v)
- {
-
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- *vector = *v.vector;
- map = v.map;
-
- *this *= a;
- }
-
-
-
- void
- Vector::equ (const TrilinosScalar a,
- const Vector &v,
- const TrilinosScalar b,
- const Vector &w)
- {
-
- Assert (numbers::is_finite(a),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
- Assert (numbers::is_finite(b),
- ExcMessage("The given value is not finite but "
- "either infinite or Not A Number (NaN)"));
-
- Assert (v.size() == w.size(),
- ExcDimensionMismatch (v.size(), w.size()));
*vector = *v.vector;
- map = v.map;
- sadd (a, b, w);
- }
-
-
-
- void
- Vector::ratio (const Vector &v,
- const Vector &w)
- {
- Assert (v.size() == w.size(),
- ExcDimensionMismatch (v.size(), w.size()));
-
- Assert (size() == w.size(),
- ExcDimensionMismatch (size(), w.size()));
-
- const int ierr = vector->ReciprocalMultiply(1.0, *(w.vector),
- *(v.vector), 0.0);
-
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- // TODO: up to now only local
- // data printed out! Find a
- // way to neatly output
- // distributed data...
- void
- Vector::print (const char *format) const
- {
- Assert (vector->GlobalLength()!=0, ExcEmptyObject());
-
- for (unsigned int j=0; j<size(); ++j)
- {
- double t = (*vector)[0][j];
-
- if (format != 0)
- std::printf (format, t);
- else
- std::printf (" %5.2f", double(t));
- }
- std::printf ("\n");
- }
-
-
-
- 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 TODO: up to
- // now only local data printed
- // out! Find a way to neatly
- // output distributed data...
- TrilinosScalar *val;
- int leading_dimension;
- int ierr = vector->ExtractView (&val, &leading_dimension);
-
- AssertThrow (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 (unsigned int i=0; i<size(); ++i)
- out << static_cast<double>(val[i]) << ' ';
- else
- for (unsigned int i=0; i<size(); ++i)
- out << static_cast<double>(val[i]) << std::endl;
- out << std::endl;
-
- // restore the representation
- // of the vector
- AssertThrow (out, ExcIO());
- }
-
-
-
- void
- Vector::swap (Vector &v)
- {
- // Just swap the pointers to
- // the two Epetra vectors that
- // hold all the data.
- Vector *p_v = &v, *p_this = this;
- Vector* tmp = p_v;
- p_v = p_this;
- p_this = tmp;
- }
-
-
-
- unsigned int
- Vector::memory_consumption () const
- {
- AssertThrow(false, ExcNotImplemented() );
- return 0;
- }
-
-
-
- LocalizedVector::LocalizedVector ()
- :
-#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
-#else
- map (0, 0, Epetra_SerialComm()),
-#endif
- vector (std::auto_ptr<Epetra_MultiVector>
- (new Epetra_MultiVector(map,1)))
- {}
-
-
- LocalizedVector::LocalizedVector (const unsigned int n)
- :
-#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- map ((int)n, 0, Epetra_MpiComm(MPI_COMM_WORLD))
-#else
- map ((int)n, 0, Epetra_SerialComm())
-#endif
- {
- reinit (n);
- }
-
-
-
- LocalizedVector::LocalizedVector (const Vector &v)
- :
-#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD))
-#else
- map (0, 0, Epetra_SerialComm())
-#endif
- {
- reinit (v);
+ return *this;
}
- LocalizedVector::LocalizedVector (const LocalizedVector &v,
- const bool fast)
- :
- map (v.map),
- vector (std::auto_ptr<Epetra_MultiVector>
- (new Epetra_MultiVector(map,1,!fast)))
- {}
-
-
-
- void
- LocalizedVector::reinit (unsigned int n)
+ template <typename Number>
+ Vector &
+ Vector::operator = (const ::dealii::Vector<Number> &v)
{
-
- if (map.NumGlobalElements() != (int)n)
+ if (size() != v.size())
{
-#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
- map = Epetra_LocalMap ((int)n, 0, Epetra_MpiComm(MPI_COMM_WORLD));
-#else
- map = Epetra_LocalMap ((int)n, 0, Epetra_SerialComm());
-#endif
+ map = LocalMap (v.size(), 0, vector->Comm());
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map,false));
}
- vector = std::auto_ptr<Epetra_MultiVector>
- (new Epetra_MultiVector (map,1,true));
- }
-
- void
- LocalizedVector::reinit (const Vector &v)
- {
- map = Epetra_LocalMap (v.size(),0,v.vector->Comm());
-
- Epetra_Import import_information (map, v.map);
-
- vector = std::auto_ptr<Epetra_MultiVector>
- (new Epetra_MultiVector (map,1,false));
-
- int ierr = vector->Import(*v.vector, import_information, Insert);
- Assert (ierr == 0, ExcTrilinosError(ierr));
- }
-
-
-
- void
- LocalizedVector::reinit (const LocalizedVector &v)
- {
- map = v.map;
- *vector = *v.vector;
- }
-
-
+ std::vector<int> indices (v.size());
+ for (unsigned int i=0; i<v.size(); ++i)
+ indices[i] = map.GID(i);
- LocalizedVector &
- LocalizedVector::operator = (const Vector &v)
- {
- reinit (v);
- return *this;
- }
+ const int ierr = vector->ReplaceGlobalValues(v.size(), 0, v.begin(),
+ &indices[0]);
+ AssertThrow (ierr == 0, ExcTrilinosError());
-
-
- LocalizedVector &
- LocalizedVector::operator = (const LocalizedVector &v)
- {
- reinit (v);
return *this;
}
--- /dev/null
+//---------------------------------------------------------------------------
+// $Id: trilinos_vector.cc 16818 2008-09-12 14:30:24Z kronbichler $
+// Version: $Name$
+//
+// Copyright (C) 2008 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------------------------------------------------------
+
+
+#include <lac/trilinos_vector_base.h>
+
+#include <cmath>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+#include <Epetra_Import.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+ namespace internal
+ {
+ VectorReference::operator TrilinosScalar () const
+ {
+ Assert (index < vector.size(),
+ ExcIndexRange (index, 0, vector.size()));
+
+ // Trilinos allows for vectors
+ // to be referenced by the [] or
+ // () operators but only ()
+ // checks index bounds. We check
+ // these bounds by ourselves, so
+ // we can use []. Note that we
+ // can only get local values.
+ AssertThrow (vector.in_local_range(index),
+ ExcAccessToNonLocalElement (index,
+ vector.vector->Map().MinMyGID(),
+ vector.vector->Map().MaxMyGID()));
+
+ return (*(vector.vector))[0][index];
+ }
+ }
+
+
+
+ VectorBase::VectorBase ()
+ :
+ last_action (Zero),
+ compressed (true),
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ vector(std::auto_ptr<Epetra_FEVector>
+ (new Epetra_FEVector(
+ Epetra_Map(0,0,Epetra_MpiComm(MPI_COMM_WORLD)))))
+#else
+ vector(std::auto_ptr<Epetra_FEVector>
+ (new Epetra_FEVector(
+ Epetra_Map(0,0,Epetra_SerialComm()))))
+#endif
+ {}
+
+
+
+ VectorBase::VectorBase (const VectorBase &v)
+ :
+ Subscriptor(),
+ last_action (Zero),
+ compressed (true),
+ vector(std::auto_ptr<Epetra_FEVector>
+ (new Epetra_FEVector(*v.vector)))
+ {}
+
+
+
+ VectorBase::~VectorBase ()
+ {}
+
+
+
+ void
+ VectorBase::clear ()
+ {
+ // When we clear the vector,
+ // reset the pointer and generate
+ // an empty vector.
+ vector.reset();
+
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+ Epetra_Map map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD));
+#else
+ Epetra_Map map (0, 0, Epetra_SerialComm());
+#endif
+
+ vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map));
+ last_action = Zero;
+ }
+
+
+
+ void
+ VectorBase::compress ()
+ {
+ // Now pass over the information about
+ // what we did last to the vector.
+ const int ierr = vector->GlobalAssemble(last_action);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ compressed = true;
+ }
+
+
+
+ VectorBase &
+ VectorBase::operator = (const TrilinosScalar s)
+ {
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but either "
+ "infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->PutScalar(s);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ bool
+ VectorBase::operator == (const VectorBase &v) const
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+ unsigned int i;
+ for (i=0; i<size(); i++)
+ if ((*(v.vector))[0][i]!=(*vector)[0][i]) return false;
+
+ return true;
+ }
+
+
+
+ bool
+ VectorBase::operator != (const VectorBase &v) const
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ return (!(*this==v));
+ }
+
+
+
+ unsigned int
+ VectorBase::size () const
+ {
+ return (unsigned int) vector->Map().NumGlobalElements();
+ }
+
+
+
+ unsigned int
+ VectorBase::local_size () const
+ {
+ return (unsigned int) vector->Map().NumMyElements();
+ }
+
+
+
+ std::pair<unsigned int, unsigned int>
+ VectorBase::local_range () const
+ {
+ int begin, end;
+ begin = vector->Map().MinMyGID();
+ end = vector->Map().MaxMyGID()+1;
+ return std::make_pair (begin, end);
+ }
+
+
+
+ TrilinosScalar
+ VectorBase::el (const unsigned int index) const
+ {
+ // Extract local indices in
+ // the vector.
+ int trilinos_i = vector->Map().LID(index);
+ TrilinosScalar value = 0.;
+
+ // If the element is not
+ // present on the current
+ // processor, we can't
+ // continue. Just print out 0.
+
+ // TODO: Is this reasonable?
+ if (trilinos_i == -1 )
+ {
+ return 0.;
+ //Assert (false, ExcAccessToNonlocalElement(index, local_range().first,
+ // local_range().second));
+ }
+ else
+ value = (*vector)[0][trilinos_i];
+
+ return value;
+ }
+
+
+
+ TrilinosScalar
+ VectorBase::operator () (const unsigned int index) const
+ {
+ // Extract local indices in
+ // the vector.
+ int trilinos_i = vector->Map().LID(index);
+ TrilinosScalar value = 0.;
+
+ // If the element is not present
+ // on the current processor, we
+ // can't continue. This is the
+ // main difference to the el()
+ // function.
+ if (trilinos_i == -1 )
+ {
+ Assert (false, ExcAccessToNonlocalElement(index, local_range().first,
+ local_range().second));
+ }
+ else
+ value = (*vector)[0][trilinos_i];
+
+ return value;
+ }
+
+
+ void
+ VectorBase::set (const std::vector<unsigned int> &indices,
+ const std::vector<TrilinosScalar> &values)
+ {
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(),values.size()));
+
+ set (indices.size(), &indices[0], &values[0]);
+ }
+
+
+
+ void
+ VectorBase::set (const std::vector<unsigned int> &indices,
+ const ::dealii::Vector<TrilinosScalar> &values)
+ {
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(),values.size()));
+
+ set (indices.size(), &indices[0], values.begin());
+ }
+
+
+
+ void
+ VectorBase::set (const unsigned int n_elements,
+ const unsigned int *indices,
+ const TrilinosScalar *values)
+ {
+ if (last_action == Add)
+ vector->GlobalAssemble(Add);
+
+ if (last_action != Insert)
+ last_action = Insert;
+
+#ifdef DEBUG
+ for (unsigned int i=0; i<n_elements; ++i)
+ if (in_local_range (indices[i]) == false)
+ {
+ compressed = false;
+ break;
+ }
+#endif
+
+ const int ierr= vector->ReplaceGlobalValues (n_elements,
+ (int*)(const_cast<unsigned int*>(indices)),
+ const_cast<TrilinosScalar*>(values));
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ void
+ VectorBase::add (const std::vector<unsigned int> &indices,
+ const std::vector<TrilinosScalar> &values)
+ {
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(),values.size()));
+
+ add (indices.size(), &indices[0], &values[0]);
+ }
+
+
+
+ void
+ VectorBase::add (const std::vector<unsigned int> &indices,
+ const ::dealii::Vector<TrilinosScalar> &values)
+ {
+ Assert (indices.size() == values.size(),
+ ExcDimensionMismatch(indices.size(),values.size()));
+
+ add (indices.size(), &indices[0], values.begin());
+ }
+
+
+
+ void
+ VectorBase::add (const unsigned int n_elements,
+ const unsigned int *indices,
+ const TrilinosScalar *values)
+ {
+ if (last_action == Insert)
+ vector->GlobalAssemble(Insert);
+
+ if (last_action != Add)
+ last_action = Add;
+
+#ifdef DEBUG
+ for (unsigned int i=0; i<n_elements; ++i)
+ if (in_local_range (indices[i]) == false)
+ {
+ compressed = false;
+ break;
+ }
+#endif
+
+ const int ierr = vector->SumIntoGlobalValues (n_elements,
+ (int*)(const_cast<unsigned int*>(indices)),
+ const_cast<TrilinosScalar*>(values));
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ TrilinosScalar
+ VectorBase::operator * (const VectorBase &vec) const
+ {
+ Assert (size() == vec.size(),
+ ExcDimensionMismatch(size(), vec.size()));
+
+ TrilinosScalar result;
+
+ const int ierr = vector->Dot(*(vec.vector), &result);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return result;
+ }
+
+
+
+ VectorBase::real_type
+ VectorBase::norm_sqr () const
+ {
+ const TrilinosScalar d = l2_norm();
+ return d*d;
+ }
+
+
+
+ TrilinosScalar
+ VectorBase::mean_value () const
+ {
+ TrilinosScalar mean;
+
+ const int ierr = vector->MeanValue (&mean);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return mean;
+ }
+
+
+
+ VectorBase::real_type
+ VectorBase::l1_norm () const
+ {
+ TrilinosScalar d;
+
+ const int ierr = vector->Norm1 (&d);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return d;
+ }
+
+
+
+ VectorBase::real_type
+ VectorBase::l2_norm () const
+ {
+ TrilinosScalar d;
+
+ const int ierr = vector->Norm2 (&d);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return d;
+ }
+
+
+
+ VectorBase::real_type
+ VectorBase::lp_norm (const TrilinosScalar p) const
+ {
+ // get a representation of the
+ // vector and loop over all
+ // the elements
+ TrilinosScalar *start_ptr;
+ int leading_dimension;
+ int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ TrilinosScalar norm = 0;
+ TrilinosScalar sum=0;
+
+ const TrilinosScalar * ptr = start_ptr;
+
+ // add up elements
+ // TODO: This
+ // won't work in parallel like
+ // this. Find out a better way to
+ // this in that case.
+ while (ptr != start_ptr+size())
+ sum += std::pow(std::fabs(*ptr++), p);
+
+ norm = std::pow(sum, static_cast<TrilinosScalar>(1./p));
+
+ return norm;
+ }
+
+
+
+ VectorBase::real_type
+ VectorBase::linfty_norm () const
+ {
+ TrilinosScalar d;
+
+ const int ierr = vector->NormInf (&d);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return d;
+ }
+
+
+
+ bool
+ VectorBase::all_zero () const
+ {
+ // get a representation of the vector and
+ // loop over all the elements
+ TrilinosScalar *start_ptr;
+ int leading_dimension;
+ int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ // TODO: This
+ // won't work in parallel like
+ // this. Find out a better way to
+ // this in that case.
+ const TrilinosScalar *ptr = start_ptr,
+ *eptr = start_ptr + size();
+ bool flag = true;
+ while (ptr != eptr)
+ {
+ if (*ptr != 0)
+ {
+ flag = false;
+ break;
+ }
+ ++ptr;
+ }
+
+ return flag;
+ }
+
+
+
+ bool
+ VectorBase::is_non_negative () const
+ {
+ // get a representation of the vector and
+ // loop over all the elements
+ TrilinosScalar *start_ptr;
+ int leading_dimension;
+ int ierr = vector->ExtractView (&start_ptr, &leading_dimension);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ // TODO: This
+ // won't work in parallel like
+ // this. Find out a better way to
+ // this in that case.
+ const TrilinosScalar *ptr = start_ptr,
+ *eptr = start_ptr + size();
+ bool flag = true;
+ while (ptr != eptr)
+ {
+ if (*ptr < 0.0)
+ {
+ flag = false;
+ break;
+ }
+ ++ptr;
+ }
+
+ return flag;
+ }
+
+
+
+ VectorBase &
+ VectorBase::operator *= (const TrilinosScalar a)
+ {
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Scale(a);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ VectorBase &
+ VectorBase::operator /= (const TrilinosScalar a)
+ {
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const TrilinosScalar factor = 1./a;
+
+ Assert (numbers::is_finite(factor),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Scale(factor);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ VectorBase &
+ VectorBase::operator += (const VectorBase &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ const int ierr = vector->Update (1.0, *(v.vector), 1.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ VectorBase &
+ VectorBase::operator -= (const VectorBase &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ const int ierr = vector->Update (-1.0, *(v.vector), 1.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ return *this;
+ }
+
+
+
+ void
+ VectorBase::add (const TrilinosScalar s)
+ {
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ unsigned int n_local = local_size();
+ int ierr;
+
+ for (unsigned int i=0; i<n_local; i++)
+ {
+ ierr = vector->SumIntoMyValue(i,0,s);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+ }
+
+
+
+ void
+ VectorBase::add (const VectorBase &v,
+ const bool allow_different_maps)
+ {
+ if (allow_different_maps == false)
+ *this += v;
+ else
+ {
+ AssertThrow (size() == v.size(),
+ ExcDimensionMismatch (size(), v.size()));
+
+ Epetra_Import data_exchange (vector->Map(), v.vector->Map());
+
+ int ierr = vector->Import(*v.vector, data_exchange, Add);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+ }
+
+
+
+ void
+ VectorBase::add (const TrilinosScalar a,
+ const VectorBase &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Update(a, *(v.vector), 1.);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ void
+ VectorBase::add (const TrilinosScalar a,
+ const VectorBase &v,
+ const TrilinosScalar b,
+ const VectorBase &w)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+ Assert (size() == w.size(),
+ ExcDimensionMismatch(size(), w.size()));
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(b),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ void
+ VectorBase::sadd (const TrilinosScalar s,
+ const VectorBase &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Update(1., *(v.vector), s);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ void
+ VectorBase::sadd (const TrilinosScalar s,
+ const TrilinosScalar a,
+ const VectorBase &v)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Update(a, *(v.vector), s);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ void
+ VectorBase::sadd (const TrilinosScalar s,
+ const TrilinosScalar a,
+ const VectorBase &v,
+ const TrilinosScalar b,
+ const VectorBase &w)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+ Assert (size() == w.size(),
+ ExcDimensionMismatch(size(), w.size()));
+
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(b),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ void
+ VectorBase::sadd (const TrilinosScalar s,
+ const TrilinosScalar a,
+ const VectorBase &v,
+ const TrilinosScalar b,
+ const VectorBase &w,
+ const TrilinosScalar c,
+ const VectorBase &x)
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+ Assert (size() == w.size(),
+ ExcDimensionMismatch(size(), w.size()));
+ Assert (size() == x.size(),
+ ExcDimensionMismatch(size(), x.size()));
+
+ Assert (numbers::is_finite(s),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(b),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(c),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ // Update member can only
+ // input two other vectors so
+ // do it in two steps
+ const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+ const int jerr = vector->Update(c, *(x.vector), 1.);
+ Assert (jerr == 0, ExcTrilinosError(jerr));
+ }
+
+
+
+ void
+ VectorBase::scale (const VectorBase &factors)
+ {
+ Assert (size() == factors.size(),
+ ExcDimensionMismatch(size(), factors.size()));
+
+ const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0);
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ void
+ VectorBase::equ (const TrilinosScalar a,
+ const VectorBase &v)
+ {
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ *vector = *v.vector;
+
+ *this *= a;
+ }
+
+
+
+ void
+ VectorBase::equ (const TrilinosScalar a,
+ const VectorBase &v,
+ const TrilinosScalar b,
+ const VectorBase &w)
+ {
+
+ Assert (v.size() == w.size(),
+ ExcDimensionMismatch (v.size(), w.size()));
+
+ Assert (numbers::is_finite(a),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+ Assert (numbers::is_finite(b),
+ ExcMessage("The given value is not finite but "
+ "either infinite or Not A Number (NaN)"));
+
+ if (vector->Map().SameAs(v.vector->Map()) == false)
+ *vector = *v.vector;
+
+ sadd (a, b, w);
+ }
+
+
+
+ void
+ VectorBase::ratio (const VectorBase &v,
+ const VectorBase &w)
+ {
+ Assert (v.size() == w.size(),
+ ExcDimensionMismatch (v.size(), w.size()));
+
+ Assert (size() == w.size(),
+ ExcDimensionMismatch (size(), w.size()));
+
+ const int ierr = vector->ReciprocalMultiply(1.0, *(w.vector),
+ *(v.vector), 0.0);
+
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ // TODO: up to now only local
+ // data printed out! Find a
+ // way to neatly output
+ // distributed data...
+ void
+ VectorBase::print (const char *format) const
+ {
+ Assert (vector->GlobalLength()!=0, ExcEmptyObject());
+
+ for (unsigned int j=0; j<size(); ++j)
+ {
+ double t = (*vector)[0][j];
+
+ if (format != 0)
+ std::printf (format, t);
+ else
+ std::printf (" %5.2f", double(t));
+ }
+ std::printf ("\n");
+ }
+
+
+
+ void
+ VectorBase::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 TODO: up to
+ // now only local data printed
+ // out! Find a way to neatly
+ // output distributed data...
+ TrilinosScalar *val;
+ int leading_dimension;
+ int ierr = vector->ExtractView (&val, &leading_dimension);
+
+ AssertThrow (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 (unsigned int i=0; i<size(); ++i)
+ out << static_cast<double>(val[i]) << ' ';
+ else
+ for (unsigned int i=0; i<size(); ++i)
+ out << static_cast<double>(val[i]) << std::endl;
+ out << std::endl;
+
+ // restore the representation
+ // of the vector
+ AssertThrow (out, ExcIO());
+ }
+
+
+
+ void
+ VectorBase::swap (VectorBase &v)
+ {
+ // Just swap the pointers to
+ // the two Epetra vectors that
+ // hold all the data.
+ VectorBase *p_v = &v, *p_this = this;
+ VectorBase* tmp = p_v;
+ p_v = p_this;
+ p_this = tmp;
+ }
+
+
+
+ unsigned int
+ VectorBase::memory_consumption () const
+ {
+ AssertThrow(false, ExcNotImplemented() );
+ return 0;
+ }
+
+} /* end of namespace TrilinosWrappers */
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS