]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Changed implementation of Trilinos vectors in a way similar to the implementation...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Sep 2008 18:11:22 +0000 (18:11 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Sep 2008 18:11:22 +0000 (18:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@16829 0785d39b-7218-0410-832d-ea1e28bc413d

15 files changed:
deal.II/lac/include/lac/trilinos_block_sparse_matrix.h
deal.II/lac/include/lac/trilinos_block_vector.h
deal.II/lac/include/lac/trilinos_precondition.h
deal.II/lac/include/lac/trilinos_precondition_amg.h
deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/include/lac/trilinos_vector.h
deal.II/lac/include/lac/trilinos_vector_base.h [new file with mode: 0644]
deal.II/lac/include/lac/vector.h
deal.II/lac/include/lac/vector.templates.h
deal.II/lac/source/trilinos_block_vector.cc [new file with mode: 0644]
deal.II/lac/source/trilinos_precondition.cc
deal.II/lac/source/trilinos_precondition_amg.cc
deal.II/lac/source/trilinos_sparse_matrix.cc
deal.II/lac/source/trilinos_vector.cc
deal.II/lac/source/trilinos_vector_base.cc [new file with mode: 0644]

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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.