From 27ade334f5b2e92c52337880c8913e07d3860e1b Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 15 Sep 2008 18:11:22 +0000 Subject: [PATCH] Changed implementation of Trilinos vectors in a way similar to the implementation of the PETSc wrapper classes, i.e., put parallel vectors in an MPI namespace and use a base class. Need to work on parallel performance a bit further. Still need to enable the proper operability of Trilinos vectors outside the lac/ directory. git-svn-id: https://svn.dealii.org/trunk@16829 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/trilinos_block_sparse_matrix.h | 219 ++- .../lac/include/lac/trilinos_block_vector.h | 430 ++++-- .../lac/include/lac/trilinos_precondition.h | 6 +- .../include/lac/trilinos_precondition_amg.h | 8 +- .../lac/include/lac/trilinos_sparse_matrix.h | 178 ++- deal.II/lac/include/lac/trilinos_vector.h | 1217 ++++------------- .../lac/include/lac/trilinos_vector_base.h | 1015 ++++++++++++++ deal.II/lac/include/lac/vector.h | 117 +- deal.II/lac/include/lac/vector.templates.h | 45 +- deal.II/lac/source/trilinos_block_vector.cc | 255 ++++ deal.II/lac/source/trilinos_precondition.cc | 6 +- .../lac/source/trilinos_precondition_amg.cc | 4 +- deal.II/lac/source/trilinos_sparse_matrix.cc | 50 +- deal.II/lac/source/trilinos_vector.cc | 1108 ++++----------- deal.II/lac/source/trilinos_vector_base.cc | 938 +++++++++++++ 15 files changed, 3555 insertions(+), 2041 deletions(-) create mode 100644 deal.II/lac/include/lac/trilinos_vector_base.h create mode 100644 deal.II/lac/source/trilinos_block_vector.cc create mode 100644 deal.II/lac/source/trilinos_vector_base.cc diff --git a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h index 18998dbefb..e50998fc2f 100644 --- a/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_block_sparse_matrix.h @@ -60,7 +60,7 @@ namespace TrilinosWrappers * @ingroup Matrix1 * @author Martin Kronbichler, Wolfgang Bangerth, 2008 */ - class BlockSparseMatrix : public BlockMatrixBase + class BlockSparseMatrix : public BlockMatrixBase { public: /** @@ -240,8 +240,20 @@ namespace TrilinosWrappers * 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 @@ -250,8 +262,21 @@ namespace TrilinosWrappers * 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 @@ -260,8 +285,21 @@ namespace TrilinosWrappers * 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 @@ -270,9 +308,21 @@ namespace TrilinosWrappers * 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$ @@ -281,9 +331,23 @@ namespace TrilinosWrappers * 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 @@ -291,6 +355,19 @@ namespace TrilinosWrappers * 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; @@ -301,6 +378,19 @@ namespace TrilinosWrappers * 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; @@ -311,6 +401,19 @@ namespace TrilinosWrappers * 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; @@ -332,7 +435,8 @@ namespace TrilinosWrappers int, int, int, int, << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3 << ',' << arg4 << "] have differing row numbers."); - /** + + /** * Exception */ DeclException4 (ExcIncompatibleColNumbers, @@ -348,6 +452,16 @@ namespace TrilinosWrappers // ------------- 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, @@ -355,6 +469,16 @@ namespace TrilinosWrappers { BaseClass::vmult_block_block (dst, src); } + + + + inline + void + BlockSparseMatrix::vmult (MPI::BlockVector &dst, + const MPI::Vector &src) const + { + BaseClass::vmult_block_nonblock (dst, src); + } @@ -365,6 +489,16 @@ namespace TrilinosWrappers { BaseClass::vmult_block_nonblock (dst, src); } + + + + inline + void + BlockSparseMatrix::vmult (MPI::Vector &dst, + const MPI::BlockVector &src) const + { + BaseClass::vmult_nonblock_block (dst, src); + } @@ -375,7 +509,17 @@ namespace TrilinosWrappers { 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 @@ -387,44 +531,85 @@ namespace TrilinosWrappers } + + 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 diff --git a/deal.II/lac/include/lac/trilinos_block_vector.h b/deal.II/lac/include/lac/trilinos_block_vector.h index 2be6d4507d..e4f1d5a946 100644 --- a/deal.II/lac/include/lac/trilinos_block_vector.h +++ b/deal.II/lac/include/lac/trilinos_block_vector.h @@ -18,7 +18,6 @@ #include #include #include -#include #include #ifdef DEAL_II_USE_TRILINOS @@ -35,9 +34,15 @@ template class BlockVector; 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 @@ -99,7 +104,7 @@ namespace TrilinosWrappers * the distribution of data among * the MPI processes. */ - explicit BlockVector (const std::vector &InputMaps); + BlockVector (const std::vector &InputMaps); /** * Copy-Constructor. Set all the @@ -133,14 +138,14 @@ namespace TrilinosWrappers * 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 @@ -162,9 +167,9 @@ namespace TrilinosWrappers * possible number type in the * deal.II vector. */ - template - BlockVector & - operator = (const ::dealii::BlockVector &V); + template + BlockVector & + operator = (const ::dealii::BlockVector &V); /** * Reinitialize the BlockVector to @@ -178,7 +183,8 @@ namespace TrilinosWrappers * If fast==false, the vector * is filled with zeros. */ - void reinit (const std::vector &input_maps); + void reinit (const std::vector &input_maps, + const bool fast = false); /** * Change the dimension to that @@ -207,7 +213,7 @@ namespace TrilinosWrappers * the wrong block. */ void reinit (const BlockVector &V, - const bool fast=false); + const bool fast = false); /** * Change the number of blocks to @@ -272,6 +278,8 @@ namespace TrilinosWrappers DeclException0 (ExcNonMatchingBlockVectors); }; + + /*----------------------- Inline functions ----------------------------------*/ @@ -311,133 +319,351 @@ namespace TrilinosWrappers + +/** + * 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; in_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 + { + public: + /** + * Typedef the base class for simpler + * access to its own typedefs. + */ + typedef BlockVectorBase 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 &InputMaps); - template - inline - BlockVector & - BlockVector::operator = (const ::dealii::BlockVector &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 + * num_blocks + * 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; in_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 + BlockVector & + operator = (const ::dealii::BlockVector &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 fast==false, the vector + * is filled with zeros. + */ + void reinit (const std::vector &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 V. The same + * applies as for the other + * reinit() function. + * + * The elements of V are not + * copied, i.e. this function is + * the same as calling reinit + * (V.size(), fast). + * + * 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 &input_maps) - { - unsigned int no_blocks = input_maps.size(); - std::vector block_sizes (no_blocks); + /** + * Change the number of blocks to + * num_blocks. 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 collect_sizes + * afterwards. + */ + void reinit (const unsigned int num_blocks); - for (unsigned int i=0; iv. 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 + * u.swap(v), 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; in_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;in_blocks();++i) - block(i).reinit(v.block(i), fast); - - collect_sizes(); - } +/*----------------------- Inline functions ----------------------------------*/ - inline - void - BlockVector::reinit (const unsigned int num_blocks) - { - std::vector 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;in_blocks();++i) - block(i).clear(); - } + inline + BlockVector::BlockVector () + {} - inline - void - BlockVector::compress () - { - for (unsigned int i=0; i &InputMaps) + { + for (unsigned int i=0; i () + { + 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; in_blocks(); ++i) - this->components[i].swap (v.components[i]); - ::dealii::swap (this->block_indices, v.block_indices); - } - /** @@ -447,7 +673,7 @@ namespace TrilinosWrappers * * @ingroup TrilinosWrappers * @relates TrilinosWrappers::BlockVector - * @author Martin Kronbichler, Wolfgang Bangerth, 2008 + * @author Martin Kronbichler, 2008 */ inline void swap (BlockVector &u, diff --git a/deal.II/lac/include/lac/trilinos_precondition.h b/deal.II/lac/include/lac/trilinos_precondition.h index c056fbd33a..aebf1abc02 100755 --- a/deal.II/lac/include/lac/trilinos_precondition.h +++ b/deal.II/lac/include/lac/trilinos_precondition.h @@ -33,7 +33,7 @@ DEAL_II_NAMESPACE_OPEN namespace TrilinosWrappers { - class Vector; + class VectorBase; class SparseMatrix; /** @@ -77,8 +77,8 @@ namespace TrilinosWrappers /** * Apply the preconditioner. */ - void vmult (Vector &dst, - const Vector &src) const; + void vmult (VectorBase &dst, + const VectorBase &src) const; private: boost::shared_ptr preconditioner; diff --git a/deal.II/lac/include/lac/trilinos_precondition_amg.h b/deal.II/lac/include/lac/trilinos_precondition_amg.h index 0b953b9011..2b27fb80fe 100755 --- a/deal.II/lac/include/lac/trilinos_precondition_amg.h +++ b/deal.II/lac/include/lac/trilinos_precondition_amg.h @@ -16,7 +16,7 @@ #include #include -#include +#include #include #include #include @@ -109,7 +109,7 @@ namespace TrilinosWrappers const bool higher_order_elements = false, const double aggregation_threshold = 1e-4, const std::vector > &null_space = std::vector > (), - const bool output_details = false); + const bool output_details = false); /** * Let Trilinos compute a @@ -164,8 +164,8 @@ namespace TrilinosWrappers * 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 diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index da91f09ee9..4df28a109d 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -19,7 +19,7 @@ #include #include #include -#include +#include #include #include @@ -45,7 +45,7 @@ DEAL_II_NAMESPACE_OPEN namespace TrilinosWrappers { // forward declarations - class Vector; + class VectorBase; class SparseMatrix; namespace MatrixIterators @@ -943,11 +943,21 @@ namespace TrilinosWrappers * 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: @@ -967,12 +977,22 @@ namespace TrilinosWrappers * 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 @@ -990,12 +1010,22 @@ namespace TrilinosWrappers * 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 @@ -1016,12 +1046,22 @@ namespace TrilinosWrappers * 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 @@ -1052,18 +1092,27 @@ namespace TrilinosWrappers * 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 @@ -1080,19 +1129,28 @@ namespace TrilinosWrappers * 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 @@ -1108,20 +1166,29 @@ namespace TrilinosWrappers * destination dst 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 matrix scaled @@ -1337,10 +1404,11 @@ namespace TrilinosWrappers -#ifndef DOXYGEN // -------------------------- inline and template functions ---------------------- +#ifndef DOXYGEN + namespace MatrixIterators { diff --git a/deal.II/lac/include/lac/trilinos_vector.h b/deal.II/lac/include/lac/trilinos_vector.h index 0dc2a13555..6d3079a853 100755 --- a/deal.II/lac/include/lac/trilinos_vector.h +++ b/deal.II/lac/include/lac/trilinos_vector.h @@ -13,36 +13,22 @@ #ifndef __deal2__trilinos_vector_h #define __deal2__trilinos_vector_h - #include #include #include #include - -#include -#include -#include +#include #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 class Vector; +// forward declaration +template class Vector; /** * @addtogroup TrilinosWrappers @@ -50,183 +36,30 @@ template class Vector; */ 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 - * v(i)=w(i). 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 s. - */ - const VectorReference & operator = (const TrilinosScalar &s) const; - - /** - * Add s to the - * referenced element of the - * vector-> - */ - const VectorReference & operator += (const TrilinosScalar &s) const; - - /** - * Subtract s from the - * referenced element of the - * vector-> - */ - const VectorReference & operator -= (const TrilinosScalar &s) const; - - /** - * Multiply the referenced - * element of the vector by - * s. - */ - const VectorReference & operator *= (const TrilinosScalar &s) const; - - /** - * Divide the referenced - * element of the vector by - * s. - */ - 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. * *

Parallel communication model

* @@ -325,23 +158,9 @@ namespace TrilinosWrappers * @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 - * C standard libraries - * vector<...> 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) @@ -352,37 +171,50 @@ namespace TrilinosWrappers * 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 fast 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 + explicit Vector (const Epetra_Map &InputMap, + const dealii::Vector &v); /** * Reinit functionality. This @@ -391,70 +223,75 @@ namespace TrilinosWrappers * 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 + * v. If fast + * is not true, 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, + * allow_different_maps, + * 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 + * v is a distributed + * vector. In this case, the + * variable fast needs + * to be set to false, + * 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 - * v=0. 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 @@ -469,8 +306,8 @@ namespace TrilinosWrappers * side vector. Otherwise, an * exception will be thrown. */ - Vector & - operator = (const LocalizedVector &V); + Vector & + operator = (const ::dealii::TrilinosWrappers::Vector &V); /** * Another copy function. This @@ -495,534 +332,101 @@ namespace TrilinosWrappers * possible number type in the * deal.II vector. */ - Vector & - operator = (const ::dealii::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. - */ - 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 - * n=local_size(). - */ - std::pair 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 el() - * command. - */ - TrilinosScalar - operator () (const unsigned int index) const; - - /** - * Return the value of the vector - * entry 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 &indices, - const std::vector &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 &indices, - const ::dealii::Vector &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 &indices, - const std::vector &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 &indices, - const ::dealii::Vector &values); - - /** - * Take an address where - * n_elements 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 - * pth root of the sum of - * the pth 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 operator +=. - */ - void add (const Vector &V); - - /** - * Simple addition of a multiple - * of a vector, i.e. *this = - * a*V. - */ - void add (const TrilinosScalar a, const Vector &V); - - /** - * Multiple addition of scaled - * vectors, i.e. *this = a*V - * + b*W. - */ - void add (const TrilinosScalar a, const Vector &V, - const TrilinosScalar b, const Vector &W); - - /** - * Scaling and simple vector - * addition, i.e. *this = - * s*(*this) + V. - */ - void sadd (const TrilinosScalar s, - const Vector &V); - - /** - * Scaling and simple addition, - * i.e. *this = s*(*this) + - * a*V. - */ - 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 + Vector & + operator = (const ::dealii::Vector &v); + private: /** - * Scaling and multiple addition. - * *this = s*(*this) + a*V + - * b*W + c*X. + * 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 *this = - * a*V. - */ - void equ (const TrilinosScalar a, const Vector &V); - /** - * Assignment *this = a*V + - * b*W. - */ - 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 this[i] = - * a[i]/b[i]. 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 b[i] 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 - * 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 - * swap(u,v) that simply - * calls u.swap(v), - * 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 + Vector::Vector (const Epetra_Map &InputMap, + const dealii::Vector &v) + : + map (InputMap) + { + vector = std::auto_ptr (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 indices (size); + for (unsigned int i=0; iReplaceGlobalValues(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 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: /** @@ -1033,14 +437,28 @@ namespace TrilinosWrappers * 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 @@ -1049,40 +467,62 @@ namespace TrilinosWrappers * 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 + explicit Vector (const dealii::Vector &v); /** * Reinit function that resizes * the vector to the size * specified by n. */ - 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 + * fast 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 @@ -1090,183 +530,90 @@ namespace TrilinosWrappers * 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 + Vector & + operator = (const ::dealii::Vector &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 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 (r); - - return *this; + u.swap (v); } +#ifndef DOXYGEN - inline - const VectorReference & - VectorReference::operator = (const TrilinosScalar &value) const + template + Vector::Vector (const dealii::Vector &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 (new Epetra_FEVector(map, false)); + std::vector indices (v.size()); + for (unsigned int i=0; iReplaceGlobalValues(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(*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(*this) / value; - vector.set (1, &index, &new_value); - return *this; - } - } - - - - inline - bool - Vector::in_local_range (const unsigned int index) const - { - std::pair 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 - (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 diff --git a/deal.II/lac/include/lac/trilinos_vector_base.h b/deal.II/lac/include/lac/trilinos_vector_base.h new file mode 100644 index 0000000000..fa8a0fdb0c --- /dev/null +++ b/deal.II/lac/include/lac/trilinos_vector_base.h @@ -0,0 +1,1015 @@ +//--------------------------------------------------------------------------- +// $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 +#include +#include +#include + +#include +#include +#include + +#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 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 + * v(i)=w(i). 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 s. + */ + const VectorReference & + operator = (const TrilinosScalar &s) const; + + /** + * Add s to the + * referenced element of the + * vector-> + */ + const VectorReference & + operator += (const TrilinosScalar &s) const; + + /** + * Subtract s from the + * referenced element of the + * vector-> + */ + const VectorReference & + operator -= (const TrilinosScalar &s) const; + + /** + * Multiply the referenced + * element of the vector by + * s. + */ + const VectorReference & + operator *= (const TrilinosScalar &s) const; + + /** + * Divide the referenced + * element of the vector by + * s. + */ + 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 + * C standard libraries + * vector<...> 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 + * reinit() 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 debug + * 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 + * v=0. 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 + * n=local_size(). + */ + std::pair 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 el() command. + */ + TrilinosScalar + operator () (const unsigned int index) const; + + /** + * Return the value of the vector + * entry 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 &indices, + const std::vector &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 &indices, + const ::dealii::Vector &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 &indices, + const std::vector &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 &indices, + const ::dealii::Vector &values); + + /** + * Take an address where + * n_elements are stored + * contiguously and add them into + * the vector. Handles all cases + * which are not covered by the + * other two add() + * 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 + * pth root of the sum of + * the pth 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 operator + * +=. + * + * Though, if the second argument + * allow_different_maps + * 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. *this = + * a*V. + */ + void add (const TrilinosScalar a, + const VectorBase &V); + + /** + * Multiple addition of scaled + * vectors, i.e. *this = a*V + * + b*W. + */ + void add (const TrilinosScalar a, + const VectorBase &V, + const TrilinosScalar b, + const VectorBase &W); + + /** + * Scaling and simple vector + * addition, i.e. *this = + * s*(*this) + V. + */ + void sadd (const TrilinosScalar s, + const VectorBase &V); + + /** + * Scaling and simple addition, + * i.e. *this = s*(*this) + + * a*V. + */ + 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. + * *this = s*(*this) + a*V + + * b*W + c*X. + */ + 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 *this = + * a*V. + */ + void equ (const TrilinosScalar a, + const VectorBase &V); + + /** + * Assignment *this = a*V + + * b*W. + */ + 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 this[i] = + * a[i]/b[i]. 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 b[i] 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 + * 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 + * swap(u,v) that simply + * calls u.swap(v), + * 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 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 (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(*this) * value; + vector.set (1, &index, &new_value); + return *this; + } + + + + inline + const VectorReference & + VectorReference::operator /= (const TrilinosScalar &value) const + { + TrilinosScalar new_value = static_cast(*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 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 ---------------------------*/ diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h index dcb6f37068..7dbab9911a 100644 --- a/deal.II/lac/include/lac/vector.h +++ b/deal.II/lac/include/lac/vector.h @@ -37,6 +37,10 @@ namespace PETScWrappers #ifdef DEAL_II_USE_TRILINOS namespace TrilinosWrappers { + namespace MPI + { + class Vector; + } class Vector; } #endif @@ -190,20 +194,36 @@ class Vector : public Subscriptor #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 @@ -371,30 +391,33 @@ class Vector : public Subscriptor #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 & 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 & operator = (const PETScWrappers::MPI::Vector &v); @@ -412,14 +435,27 @@ class Vector : public Subscriptor * 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 & + 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 & operator = (const TrilinosWrappers::Vector &v); @@ -427,11 +463,12 @@ class Vector : public Subscriptor /** * 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 bool operator == (const Vector &v) const; diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h index e3d336e386..b3f6c2515e 100644 --- a/deal.II/lac/include/lac/vector.templates.h +++ b/deal.II/lac/include/lac/vector.templates.h @@ -172,7 +172,7 @@ Vector::Vector (const PETScWrappers::MPI::Vector &v) #ifdef DEAL_II_USE_TRILINOS template -Vector::Vector (const TrilinosWrappers::Vector &v) +Vector::Vector (const TrilinosWrappers::MPI::Vector &v) : Subscriptor(), vec_size(v.size()), @@ -190,7 +190,7 @@ Vector::Vector (const TrilinosWrappers::Vector &v) // 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 @@ -203,6 +203,32 @@ Vector::Vector (const TrilinosWrappers::Vector &v) } } + + +template +Vector::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 @@ -671,6 +697,21 @@ Vector::operator = (const PETScWrappers::MPI::Vector &v) #ifdef DEAL_II_USE_TRILINOS +template +Vector & +Vector::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 Vector & Vector::operator = (const TrilinosWrappers::Vector &v) diff --git a/deal.II/lac/source/trilinos_block_vector.cc b/deal.II/lac/source/trilinos_block_vector.cc new file mode 100644 index 0000000000..b54adda76d --- /dev/null +++ b/deal.II/lac/source/trilinos_block_vector.cc @@ -0,0 +1,255 @@ +//--------------------------------------------------------------------------- +// $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 + + +#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; in_blocks(); ++i) + this->components[i] = v.block(i); + + collect_sizes(); + + return *this; + } + + + + template + BlockVector & + BlockVector::operator = (const ::dealii::BlockVector &v) + { + Assert (n_blocks() == v.n_blocks(), + ExcDimensionMismatch(n_blocks(),v.n_blocks())); + + for (unsigned int i=0; in_blocks(); ++i) + this->components[i] = v.block(i); + + return *this; + } + + + + BlockVector::~BlockVector () + {} + + + + void + BlockVector::reinit (const std::vector &input_maps, + const bool fast) + { + unsigned int no_blocks = input_maps.size(); + std::vector block_sizes (no_blocks); + + for (unsigned int i=0; iblock_indices.reinit (block_sizes); + if (components.size() != n_blocks()) + components.resize(n_blocks()); + + for (unsigned int i=0; i 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;in_blocks();++i) + block(i).clear(); + } + + + + void + BlockVector::compress () + { + for (unsigned int i=0; in_blocks() == v.n_blocks(), + ExcDimensionMismatch(this->n_blocks(), v.n_blocks())); + + for (unsigned int i=0; in_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 &input_maps, + const bool fast) + { + unsigned int no_blocks = input_maps.size(); + std::vector block_sizes (no_blocks); + + for (unsigned int i=0; iblock_indices.reinit (block_sizes); + if (components.size() != n_blocks()) + components.resize(n_blocks()); + + for (unsigned int i=0; i 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;in_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;in_blocks(); ++i) + this->components[i] = v.block(i); + + return *this; + } + +} + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/deal.II/lac/source/trilinos_precondition.cc b/deal.II/lac/source/trilinos_precondition.cc index a079024d52..0d9eb8768f 100755 --- a/deal.II/lac/source/trilinos_precondition.cc +++ b/deal.II/lac/source/trilinos_precondition.cc @@ -13,7 +13,7 @@ #include -#include +#include #include #ifdef DEAL_II_USE_TRILINOS @@ -59,8 +59,8 @@ namespace TrilinosWrappers void - PreconditionSSOR::vmult (Vector &dst, - const Vector &src) const + PreconditionSSOR::vmult (VectorBase &dst, + const VectorBase &src) const { preconditioner->ApplyInverse (*src.vector, *dst.vector); } diff --git a/deal.II/lac/source/trilinos_precondition_amg.cc b/deal.II/lac/source/trilinos_precondition_amg.cc index c265a76b12..ea61218f94 100755 --- a/deal.II/lac/source/trilinos_precondition_amg.cc +++ b/deal.II/lac/source/trilinos_precondition_amg.cc @@ -147,8 +147,8 @@ namespace TrilinosWrappers - 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); diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 6e70110c5b..5881615a48 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -11,7 +11,7 @@ // //--------------------------------------------------------------------------- -#include +#include #include #ifdef DEAL_II_USE_TRILINOS @@ -432,9 +432,11 @@ namespace TrilinosWrappers 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; @@ -467,10 +469,11 @@ namespace TrilinosWrappers 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 @@ -838,9 +841,10 @@ namespace TrilinosWrappers } + void - SparseMatrix::vmult (Vector &dst, - const Vector &src) const + SparseMatrix::vmult (VectorBase &dst, + const VectorBase &src) const { Assert (&src != &dst, ExcSourceEqualsDestination()); @@ -859,8 +863,8 @@ namespace TrilinosWrappers void - SparseMatrix::Tvmult (Vector &dst, - const Vector &src) const + SparseMatrix::Tvmult (VectorBase &dst, + const VectorBase &src) const { Assert (&src != &dst, ExcSourceEqualsDestination()); @@ -879,13 +883,12 @@ namespace TrilinosWrappers 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; } @@ -893,13 +896,12 @@ namespace TrilinosWrappers 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; } @@ -907,9 +909,9 @@ namespace TrilinosWrappers 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; } @@ -917,10 +919,10 @@ namespace TrilinosWrappers 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; } @@ -928,9 +930,9 @@ namespace TrilinosWrappers 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; diff --git a/deal.II/lac/source/trilinos_vector.cc b/deal.II/lac/source/trilinos_vector.cc index cfc62b919f..adfe635e7a 100755 --- a/deal.II/lac/source/trilinos_vector.cc +++ b/deal.II/lac/source/trilinos_vector.cc @@ -25,988 +25,388 @@ DEAL_II_NAMESPACE_OPEN 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(index) >= vector.map.MinMyGID()) && - (static_cast(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 - (new Epetra_FEVector(map))) - {} + { + last_action = Zero; + vector = std::auto_ptr (new Epetra_FEVector(map)); + } - Vector::Vector (const Epetra_Map &InputMap) - : - map (InputMap), - last_action (Insert), - vector (std::auto_ptr - (new Epetra_FEVector(map))) - {} + Vector::Vector (const Epetra_Map &InputMap) + : + map (InputMap) + { + last_action = Zero; + vector = std::auto_ptr (new Epetra_FEVector(map)); + } - Vector::Vector (const Vector &v, - const bool fast) - : - Subscriptor(), - map (v.map), - last_action (Insert), - vector(std::auto_ptr - (new Epetra_FEVector(map,!fast))) - {} - - - - Vector::~Vector () - { - } - - - - void - Vector::reinit (const Epetra_Map &input_map) - { - vector.reset(); - map = input_map; - - vector = std::auto_ptr (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 (new Epetra_FEVector(map,!fast)); - last_action = Insert; - } - - + Vector::Vector (const Vector &v) + : + VectorBase(), + map (v.map) + { + last_action = Zero; + vector = std::auto_ptr (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 (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 (new Epetra_FEVector(*v.vector)); + else + { + vector = std::auto_ptr (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 (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 (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 (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 &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 + Vector & + Vector::operator = (const ::dealii::Vector &v) + { + Assert (size() == v.size(), + ExcDimensionMismatch(size(), v.size())); - std::vector local_indices (n_local_elements); - for (unsigned int i=0; i (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 indices (size); + for (unsigned int i=0; iReplaceGlobalValues(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 + (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 + (new Epetra_FEVector(map,true)); } - std::pair - 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 (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 &indices, - const std::vector &values) - { - Assert (indices.size() == values.size(), - ExcDimensionMismatch(indices.size(),values.size())); - - set (indices.size(), &indices[0], &values[0]); } void - Vector::set (const std::vector &indices, - const ::dealii::Vector &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(indices)), - const_cast(values)); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - } + last_action = Zero; - - - void - Vector::add (const std::vector &indices, - const std::vector &values) - { - Assert (indices.size() == values.size(), - ExcDimensionMismatch(indices.size(),values.size())); - - add (indices.size(), &indices[0], &values[0]); + vector = std::auto_ptr (new Epetra_FEVector (map,true)); } void - Vector::add (const std::vector &indices, - const ::dealii::Vector &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(indices)), - const_cast(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 (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(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 (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 (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; iMap().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 (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; jExtractView (&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(val[i]) << ' '; - else - for (unsigned int i=0; i(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 - (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 - (new Epetra_MultiVector(map,1,!fast))) - {} - - - - void - LocalizedVector::reinit (unsigned int n) + template + Vector & + Vector::operator = (const ::dealii::Vector &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 (new Epetra_FEVector(map,false)); } - vector = std::auto_ptr - (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 - (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 indices (v.size()); + for (unsigned int i=0; iReplaceGlobalValues(v.size(), 0, v.begin(), + &indices[0]); + AssertThrow (ierr == 0, ExcTrilinosError()); - - - LocalizedVector & - LocalizedVector::operator = (const LocalizedVector &v) - { - reinit (v); return *this; } diff --git a/deal.II/lac/source/trilinos_vector_base.cc b/deal.II/lac/source/trilinos_vector_base.cc new file mode 100644 index 0000000000..4d841242dd --- /dev/null +++ b/deal.II/lac/source/trilinos_vector_base.cc @@ -0,0 +1,938 @@ +//--------------------------------------------------------------------------- +// $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 + +#include + +#ifdef DEAL_II_USE_TRILINOS + +#include + + +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 + (new Epetra_FEVector( + Epetra_Map(0,0,Epetra_MpiComm(MPI_COMM_WORLD))))) +#else + vector(std::auto_ptr + (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 + (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 (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; iMap().NumGlobalElements(); + } + + + + unsigned int + VectorBase::local_size () const + { + return (unsigned int) vector->Map().NumMyElements(); + } + + + + std::pair + 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 &indices, + const std::vector &values) + { + Assert (indices.size() == values.size(), + ExcDimensionMismatch(indices.size(),values.size())); + + set (indices.size(), &indices[0], &values[0]); + } + + + + void + VectorBase::set (const std::vector &indices, + const ::dealii::Vector &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; iReplaceGlobalValues (n_elements, + (int*)(const_cast(indices)), + const_cast(values)); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + VectorBase::add (const std::vector &indices, + const std::vector &values) + { + Assert (indices.size() == values.size(), + ExcDimensionMismatch(indices.size(),values.size())); + + add (indices.size(), &indices[0], &values[0]); + } + + + + void + VectorBase::add (const std::vector &indices, + const ::dealii::Vector &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; iSumIntoGlobalValues (n_elements, + (int*)(const_cast(indices)), + const_cast(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(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; iSumIntoMyValue(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; jExtractView (&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(val[i]) << ' '; + else + for (unsigned int i=0; i(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 -- 2.39.5