]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added a rudimentary implementation of Trilinos sparse matrix and vector wrappers...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 21 Aug 2008 01:07:49 +0000 (01:07 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 21 Aug 2008 01:07:49 +0000 (01:07 +0000)
git-svn-id: https://svn.dealii.org/trunk@16617 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_sparse_matrix.h [new file with mode: 0755]
deal.II/lac/include/lac/trilinos_vector.h [new file with mode: 0755]
deal.II/lac/source/trilinos_sparse_matrix.cc [new file with mode: 0755]
deal.II/lac/source/trilinos_vector.cc [new file with mode: 0755]

diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h
new file mode 100755 (executable)
index 0000000..dee47ba
--- /dev/null
@@ -0,0 +1,1220 @@
+//---------------------------------------------------------------------------
+//    $Id: trilinos_sparse_matrix.h 16014 2008-04-30 16:18:03Z young $
+//    Version: $Name$
+//
+//    Copyright (C) 2004, 2005, 2006, 2007, 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_sparse_matrix_h
+#define __deal2__trilinos_sparse_matrix_h
+
+
+#include <base/config.h>
+#include <base/subscriptor.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/exceptions.h>
+#include <lac/trilinos_vector.h>
+
+#include <boost/shared_ptr.hpp>
+#include <vector>
+#include <cmath>
+#include <memory>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+#  define TrilinosScalar double
+#  include <Epetra_FECrsMatrix.h>
+#  include <Epetra_Map.h>
+#  include <Epetra_CrsGraph.h>
+#  ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+#    include <Epetra_MpiComm.h>
+#    include "mpi.h"
+#  else
+#    include "Epetra_SerialComm.h"
+#  endif
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace TrilinosWrappers
+{
+                                   // forward declarations
+  class Vector;
+  class SparseMatrix;
+
+  namespace MatrixIterators
+  {
+/**
+ * STL conforming iterator. This class acts as an iterator walking over the
+ * elements of Trilinos matrices. The implementation of this class is similar
+ * to the one for PETSc matrices.
+ *
+ * Note that Trilinos does not give any guarantees as to the order of elements
+ * within each row. Note also that accessing the elements of a full matrix
+ * surprisingly only shows the nonzero elements of the matrix, not all
+ * elements.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ */
+    class const_iterator
+    {
+      private:
+                                         /**
+                                          * Accessor class for iterators
+                                          */
+        class Accessor
+        {
+          public:
+                                             /**
+                                              * Constructor. Since we use
+                                              * accessors only for read
+                                              * access, a const matrix
+                                              * pointer is sufficient.
+                                              */
+            Accessor (const SparseMatrix  *matrix,
+                      const unsigned int   row,
+                      const unsigned int   index);
+
+                                             /**
+                                              * Row number of the element
+                                              * represented by this
+                                              * object.
+                                              */
+            unsigned int row() const;
+
+                                             /**
+                                              * Index in row of the element
+                                              * represented by this
+                                              * object.
+                                              */
+            unsigned int index() const;
+
+                                             /**
+                                              * Column number of the
+                                              * element represented by
+                                              * this object.
+                                              */
+            unsigned int column() const;
+
+                                             /**
+                                              * Value of this matrix entry.
+                                              */
+            TrilinosScalar value() const;
+
+                                             /**
+                                              * Exception
+                                              */
+            DeclException0 (ExcBeyondEndOfMatrix);
+                                             /**
+                                              * Exception
+                                              */
+            DeclException3 (ExcAccessToNonlocalRow,
+                            int, int, int,
+                            << "You tried to access row " << arg1
+                            << " of a distributed matrix, but only rows "
+                            << arg2 << " through " << arg3
+                            << " are stored locally and can be accessed.");
+            
+          private:
+                                             /**
+                                              * The matrix accessed.
+                                              */
+            mutable SparseMatrix *matrix;
+
+                                             /**
+                                              * Current row number.
+                                              */
+            unsigned int a_row;
+
+                                             /**
+                                              * Current index in row.
+                                              */
+            unsigned int a_index;
+
+                                             /**
+                                              * Cache where we store the
+                                              * column indices of the present
+                                              * row. This is necessary, since
+                                              * Trilinos makes access to the
+                                              * elements of its matrices
+                                              * rather hard, and it is much
+                                              * more efficient to copy all
+                                              * column entries of a row once
+                                              * when we enter it than
+                                              * repeatedly asking Trilinos for
+                                              * individual ones. This also
+                                              * makes some sense since it is
+                                              * likely that we will access
+                                              * them sequentially anyway.
+                                              *
+                                              * In order to make copying of
+                                              * iterators/accessor of
+                                              * acceptable performance, we
+                                              * keep a shared pointer to these
+                                              * entries so that more than one
+                                              * accessor can access this data
+                                              * if necessary.
+                                              */
+            boost::shared_ptr<const std::vector<unsigned int> > colnum_cache;
+
+                                             /**
+                                              * Similar cache for the values
+                                              * of this row.
+                                              */
+            boost::shared_ptr<const std::vector<TrilinosScalar> > value_cache;
+            
+                                             /**
+                                              * Discard the old row caches
+                                              * (they may still be used by
+                                              * other accessors) and generate
+                                              * new ones for the row pointed
+                                              * to presently by this accessor.
+                                              */
+            void visit_present_row ();
+
+                                             /**
+                                              * Make enclosing class a
+                                              * friend.
+                                              */
+            friend class const_iterator;
+        };
+        
+      public:
+          
+                                         /**
+                                          * Constructor. Create an iterator
+                                          * into the matrix @p matrix for the
+                                          * given row and the index within it.
+                                          */ 
+        const_iterator (const SparseMatrix   *matrix,
+                        const unsigned int  row,
+                        const unsigned int  index);
+          
+                                         /**
+                                          * Prefix increment.
+                                          */
+        const_iterator& operator++ ();
+
+                                         /**
+                                          * Postfix increment.
+                                          */
+        const_iterator operator++ (int);
+
+                                         /**
+                                          * Dereferencing operator.
+                                          */
+        const Accessor& operator* () const;
+
+                                         /**
+                                          * Dereferencing operator.
+                                          */
+        const Accessor* operator-> () const;
+
+                                         /**
+                                          * Comparison. True, if
+                                          * both iterators point to
+                                          * the same matrix
+                                          * position.
+                                          */
+        bool operator == (const const_iterator&) const;
+                                         /**
+                                          * Inverse of <tt>==</tt>.
+                                          */
+        bool operator != (const const_iterator&) const;
+
+                                         /**
+                                          * Comparison
+                                          * operator. Result is true
+                                          * if either the first row
+                                          * number is smaller or if
+                                          * the row numbers are
+                                          * equal and the first
+                                          * index is smaller.
+                                          */
+        bool operator < (const const_iterator&) const;
+
+                                         /**
+                                          * Exception
+                                          */
+        DeclException2 (ExcInvalidIndexWithinRow,
+                        int, int,
+                        << "Attempt to access element " << arg2
+                        << " of row " << arg1
+                        << " which doesn't have that many elements.");
+        
+      private:
+                                         /**
+                                          * Store an object of the
+                                          * accessor class.
+                                          */
+        Accessor accessor;
+    };
+    
+  }
+  
+  
+/**
+ * This class implements a wrapper to use the Trilinos distributed sparse matrix
+ * class Epetra_FECrsMatrix. This is precisely the kind of matrix we deal with
+ * all the time - we most likely 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
+ * SparseMatrix 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 doubles.
+ *
+ * Note that Trilinos only guarantees that operations do what you expect if the
+ * functions @p GlobalAssemble  has been called
+ * after matrix assembly. Therefore, you need to call
+ * SparseMatrix::compress() before you actually use the matrix. This also
+ * calls @p FillComplete that compresses the storage format for sparse
+ * matrices by discarding unused elements. Trilinos allows to continue with
+ * assembling the matrix after calls to these functions, but since there are
+ * no more free entries available after that any more, it is better to only
+ * call SparseMatrix::compress() once at the end of the assembly stage and
+ * before the matrix is actively used.
+ * 
+ * @ingroup TrilinosWrappers
+ * @ingroup Matrix1
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ */
+  class SparseMatrix : public Subscriptor
+  {
+    public:
+                                       /**
+                                        * Declare a typedef for the iterator
+                                        * class.
+                                        */
+      typedef MatrixIterators::const_iterator const_iterator;
+
+                                       /**
+                                        * Declare a typedef in analogy to all
+                                        * the other container classes.
+                                        */
+      typedef TrilinosScalar value_type;
+      
+                                       /**
+                                        * Default constructor.
+                                        */
+      SparseMatrix ();
+
+                                       /**
+                                        * Constructor using an Epetra_Map
+                                       * and a maximum number of nonzero
+                                       * matrix entries.
+                                        */
+      SparseMatrix (const Epetra_Map   &InputMap,
+                   const unsigned int  n_max_entries_per_row);
+
+                                       /**
+                                        * Same as before, but now use the
+                                       * exact number of nonzeros in each
+                                       * matrix row. Since we know the
+                                       * number of elements in the matrix
+                                       * exactly in this case, we can 
+                                       * already allocate the right amount
+                                       * of memory, which makes the 
+                                       * creation process including the
+                                       * insertion of nonzero elements
+                                       * by the respective 
+                                       * SparseMatrix::reinit call 
+                                       * considerably faster.
+                                        */
+      SparseMatrix (const Epetra_Map                &InputMap,
+                   const std::vector<unsigned int> &n_entries_per_row);
+
+                                       /**
+                                        * Destructor. Made virtual so that one
+                                        * can use pointers to this class.
+                                        */
+      virtual ~SparseMatrix ();
+                                       /**
+                                        * This operator assigns a scalar to a
+                                        * matrix. Since this does usually not
+                                        * make much sense (should we set all
+                                        * matrix entries to this value? Only
+                                        * the nonzero entries of the sparsity
+                                        * pattern?), this operation is only
+                                        * allowed if the actual value to be
+                                        * assigned is zero. This operator only
+                                        * exists to allow for the obvious
+                                        * notation <tt>matrix=0</tt>, which
+                                        * sets all elements of the matrix to
+                                        * zero, but keeps the sparsity pattern
+                                        * previously used.
+                                        */
+      SparseMatrix &
+      operator = (const double d);
+
+                                       /**
+                                        * This function initializes the
+                                       * Trilinos matrix by attaching all
+                                       * the elements to the sparsity
+                                       * pattern provided as deal argument.
+                                       * This function uses a user-
+                                       * provided maximum number of 
+                                       * elements per row. If that is
+                                       * not directly available, use one of the
+                                       * other reinit functions.
+                                        */
+      void reinit (const CompressedSparsityPattern &sparsity_pattern,
+                  const unsigned int               n_max_entries_per_row);
+
+                                       /**
+                                        * This function initializes the
+                                       * Trilinos matrix by attaching all
+                                       * the elements to the sparsity
+                                       * pattern provided as deal argument,
+                                       * now calculating the maximum number
+                                       * of nonzeros from the sparsity
+                                       * pattern internally.
+                                        */
+      void reinit (const CompressedSparsityPattern &sparsity_pattern);
+
+                                      /**
+                                        * This function is similar to the
+                                       * other initialization function above,
+                                       * but now also reassigns the matrix 
+                                       * rows according to a user-supplied
+                                       * Epetra map. This might be used
+                                       * when the matrix structure changes,
+                                       * e.g. when the grid is refined.
+                                        */
+      void reinit (const Epetra_Map                &input_map,
+                  const CompressedSparsityPattern &sparsity_pattern);
+
+                                       /**
+                                        * Release all memory and return
+                                        * to a state just like after
+                                        * having called the default
+                                        * constructor.
+                                        */
+      void clear ();
+
+                                       /**
+                                        * Set the element (<i>i,j</i>)
+                                        * to @p value.
+                                       *
+                                       * This function
+                                       * adds a new entry to the
+                                       * matrix if it didn't exist
+                                       * before, very much in
+                                       * contrast to the SparseMatrix
+                                       * class which throws an error
+                                       * if the entry does not exist.
+                                       * If <tt>value</tt> is not a
+                                       * finite number an exception
+                                       * is thrown.
+                                       */
+      void set (const unsigned int i,
+                const unsigned int j,
+                const TrilinosScalar value);
+
+                                       /**
+                                        * Add @p value to the
+                                        * element (<i>i,j</i>).
+                                       *
+                                       * This function
+                                       * adds a new entry to the
+                                       * matrix if it didn't exist
+                                       * before, very much in
+                                       * contrast to the SparseMatrix
+                                       * class which throws an error
+                                       * if the entry does not exist.
+                                       * If <tt>value</tt> is not a
+                                       * finite number an exception
+                                       * is thrown.
+                                        */
+      void add (const unsigned int i,
+                const unsigned int j,
+                const TrilinosScalar value);
+
+                                       /**
+                                        * Remove all elements from
+                                        * this <tt>row</tt> by setting
+                                        * them to zero. The function
+                                        * does not modify the number
+                                        * of allocated nonzero
+                                        * entries, it only sets some
+                                        * entries to zero. It may drop
+                                        * them from the sparsity
+                                        * pattern, though (but retains
+                                        * the allocated memory in case
+                                        * new entries are again added
+                                        * later).
+                                        *
+                                        * This operation is used in
+                                        * eliminating constraints (e.g. due to
+                                        * hanging nodes) and makes sure that
+                                        * we can write this modification to
+                                        * the matrix without having to read
+                                        * entries (such as the locations of
+                                        * non-zero elements) from it --
+                                        * without this operation, removing
+                                        * constraints on parallel matrices is
+                                        * a rather complicated procedure.
+                                        *
+                                        * The second parameter can be used to
+                                        * set the diagonal entry of this row
+                                        * to a value different from zero. The
+                                        * default is to set it to zero.
+                                        */
+     // void clear_row (const unsigned int row,
+     //                 const TrilinosScalar  new_diag_value = 0);
+
+                                       /**
+                                        * Same as clear_row(), except that it
+                                        * works on a number of rows at once.
+                                        *
+                                        * The second parameter can be used to
+                                        * set the diagonal entries of all
+                                        * cleared rows to something different
+                                        * from zero. Note that all of these
+                                        * diagonal entries get the same value
+                                        * -- if you want different values for
+                                        * the diagonal entries, you have to
+                                        * set them by hand.
+                                        */
+     // void clear_rows (const std::vector<unsigned int> &rows,
+     //                  const TrilinosScalar                new_diag_value = 0);
+      
+                                       /**
+                                        * Trilinos matrices store their own
+                                        * sparsity patterns. So, in analogy to
+                                        * our own SparsityPattern class,
+                                        * this function compresses the
+                                        * sparsity pattern and allows the
+                                        * resulting matrix to be used in all
+                                        * other operations where before only
+                                        * assembly functions were
+                                        * allowed. This function must
+                                        * therefore be called once you have
+                                        * assembled the matrix.
+                                        */
+      void compress ();
+      
+                                       /**
+                                        * Return the value of the entry
+                                        * (<i>i,j</i>).  This may be an
+                                        * expensive operation and you should
+                                        * always take care where to call this
+                                        * function. In contrast to the
+                                        * respective function in the
+                                        * @p SparseMatrix class, we don't
+                                        * throw an exception if the respective
+                                        * entry doesn't exist in the sparsity
+                                        * pattern of this class, since Trilinos
+                                        * does not transmit this information.
+                                       * On the other hand, an exception
+                                       * will be thrown when the requested
+                                       * element is not saved on the calling
+                                       * process.
+                                        *
+                                        * This function is therefore exactly
+                                        * equivalent to the <tt>el()</tt> function.
+                                        */
+      TrilinosScalar operator () (const unsigned int i,
+                                 const unsigned int j) const;
+
+                                       /**
+                                        * Return the value of the matrix entry
+                                        * (<i>i,j</i>). If this entry does not
+                                        * exist in the sparsity pattern, then
+                                        * zero is returned. While this may be
+                                        * convenient in some cases, note that
+                                        * it is simple to write algorithms
+                                        * that are slow compared to an optimal
+                                        * solution, since the sparsity of the
+                                        * matrix is not used.
+                                        */
+      TrilinosScalar el (const unsigned int i,
+                        const unsigned int j) const;
+
+                                       /**
+                                        * Return the main diagonal
+                                        * element in the <i>i</i>th
+                                        * row. This function throws an
+                                        * error if the matrix is not
+                                        * quadratic.
+                                        *
+                                        * TODO: Trilinos can access the 
+                                       * diagonal faster. Implement this!
+                                        */
+      TrilinosScalar diag_element (const unsigned int i) const;
+      
+                                       /**
+                                        * Return the number of rows in this
+                                        * matrix.
+                                        */
+      unsigned int m () const;
+
+                                       /**
+                                        * Return the number of columns in this
+                                        * matrix.
+                                        */
+      unsigned int n () const;
+
+                                       /**
+                                        * Return the local dimension of the
+                                        * matrix, i.e. the number of rows
+                                        * stored on the present MPI
+                                        * process. For sequential matrices,
+                                        * this number is the same as m(),
+                                        * but for parallel matrices 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 rows of
+                                       * this matrix are stored
+                                       * locally. The first number is
+                                       * the index of the first
+                                       * row stored, the second
+                                       * the index of the one past
+                                       * the last one that is stored
+                                       * locally. If this is a
+                                       * sequential matrix, then the
+                                       * result will be the pair
+                                       * (0,m()), 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;
+
+                                       /**
+                                        * Return the number of nonzero
+                                        * elements of this
+                                        * matrix. 
+                                        */
+      unsigned int n_nonzero_elements () const;
+
+                                       /**
+                                        * Number of entries in a specific row.
+                                        */
+      unsigned int row_length (const unsigned int row) const;
+      
+                                       /**
+                                        * Return the l1-norm of the matrix, that is
+                                        * $|M|_1=max_{all columns j}\sum_{all 
+                                        * rows i} |M_ij|$,
+                                        * (max. sum of columns).
+                                        * This is the
+                                        * natural matrix norm that is compatible
+                                        * to the l1-norm for vectors, i.e.
+                                        * $|Mv|_1\leq |M|_1 |v|_1$.
+                                        * (cf. Haemmerlin-Hoffmann:
+                                        * Numerische Mathematik)
+                                        */
+      TrilinosScalar l1_norm () const;
+
+                                       /**
+                                        * Return the linfty-norm of the
+                                        * matrix, that is
+                                        * $|M|_infty=max_{all rows i}\sum_{all 
+                                        * columns j} |M_ij|$,
+                                        * (max. sum of rows).
+                                        * This is the
+                                        * natural matrix norm that is compatible
+                                        * to the linfty-norm of vectors, i.e.
+                                        * $|Mv|_infty \leq |M|_infty |v|_infty$.
+                                        * (cf. Haemmerlin-Hoffmann:
+                                        * Numerische Mathematik)
+                                        */
+      TrilinosScalar linfty_norm () const;
+
+                                       /**
+                                        * Return the frobenius norm of the
+                                        * matrix, i.e. the square root of the
+                                        * sum of squares of all entries in the
+                                        * matrix.
+                                        */
+      TrilinosScalar frobenius_norm () const;
+      
+                                       /**
+                                        * Multiply the entire matrix by a
+                                        * fixed factor.
+                                        */
+      SparseMatrix & operator *= (const TrilinosScalar factor);
+    
+                                       /**
+                                        * Divide the entire matrix by a
+                                        * fixed factor.
+                                        */
+      SparseMatrix & operator /= (const TrilinosScalar factor);
+
+                                       /**
+                                        * Matrix-vector multiplication:
+                                        * let <i>dst = M*src</i> with
+                                        * <i>M</i> being this matrix.
+                                        *
+                                        * Source and destination must
+                                        * not be the same vector.
+                                       *
+                                       * 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.
+                                        */
+      void vmult (Vector       &dst,
+                  const Vector &src) const;
+
+                                       /**
+                                        * Matrix-vector multiplication: let
+                                        * <i>dst = M<sup>T</sup>*src</i> with
+                                        * <i>M</i> being this matrix. This
+                                        * function does the same as vmult()
+                                        * but takes the transposed matrix.
+                                        *
+                                        * Source and destination must
+                                        * not be the same vector.
+                                       *
+                                       * 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.
+                                        */
+      void Tvmult (Vector       &dst,
+                   const Vector &src) const;
+
+                                       /**
+                                        * Adding Matrix-vector
+                                        * multiplication. Add
+                                        * <i>M*src</i> on <i>dst</i>
+                                        * with <i>M</i> being this
+                                        * matrix.
+                                        *
+                                        * Source and destination must
+                                        * not be the same vector.
+                                       *
+                                       * 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.
+                                        */
+      void vmult_add (Vector       &dst,
+                      const Vector &src) const;
+
+                                       /**
+                                        * Adding Matrix-vector
+                                        * multiplication. Add
+                                        * <i>M<sup>T</sup>*src</i> to
+                                        * <i>dst</i> with <i>M</i> being
+                                        * this matrix. This function
+                                        * does the same as vmult_add()
+                                        * but takes the transposed
+                                        * matrix.
+                                        *
+                                        * Source and destination must
+                                        * not be the same vector.
+                                       *
+                                       * 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.
+                                        */
+      void Tvmult_add (Vector       &dst,
+                       const Vector &src) const;
+
+                                       /**
+                                        * Return the square of the norm
+                                        * of the vector $v$ with respect
+                                        * to the norm induced by this
+                                        * matrix,
+                                        * i.e. $\left(v,Mv\right)$. This
+                                        * is useful, e.g. in the finite
+                                        * element context, where the
+                                        * $L_2$ norm of a function
+                                        * equals the matrix norm with
+                                        * respect to the mass matrix of
+                                        * the vector representing the
+                                        * nodal values of the finite
+                                        * element function.
+                                        *
+                                        * Obviously, the matrix needs to
+                                        * be quadratic for this operation.
+                                        *
+                                        * The implementation of this function
+                                        * is not as efficient as the one in
+                                        * the @p SparseMatrix class used in
+                                        * deal.II (i.e. the original one, not
+                                        * the Trilinos wrapper class) since Trilinos
+                                        * doesn't support this operation and
+                                        * needs a temporary vector.
+                                       *
+                                       * 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.
+                                        */
+      TrilinosScalar matrix_norm_square (const Vector &v) const;
+
+                                       /**
+                                        * Compute the matrix scalar
+                                        * product $\left(u,Mv\right)$.
+                                        *
+                                        * The implementation of this function
+                                        * is not as efficient as the one in
+                                        * the @p SparseMatrix class used in
+                                        * deal.II (i.e. the original one, not
+                                        * the Trilinos wrapper class) since Trilinos
+                                        * doesn't support this operation and
+                                        * needs a temporary vector.
+                                       *
+                                       * 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.
+                                        */
+      TrilinosScalar matrix_scalar_product (const Vector &u,
+                                           const Vector &v) const;
+
+                                       /**
+                                        * Compute the residual of an
+                                        * equation <i>Mx=b</i>, where
+                                        * the residual is defined to be
+                                        * <i>r=b-Mx</i>. Write the
+                                        * residual into
+                                        * @p dst. The
+                                        * <i>l<sub>2</sub></i> norm of
+                                        * the residual vector is
+                                        * returned.
+                                        *
+                                        * Source <i>x</i> and destination
+                                        * <i>dst</i> must not be the same
+                                        * vector.
+                                       *
+                                       * 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.
+                                        */
+      TrilinosScalar residual (Vector       &dst,
+                              const Vector &x,
+                              const Vector &b) const;
+
+                                       /**
+                                        * STL-like iterator with the
+                                        * first entry.
+                                        */
+      const_iterator begin () const;
+
+                                       /**
+                                        * Final iterator.
+                                        */
+      const_iterator end () const;
+
+                                       /**
+                                        * STL-like iterator with the
+                                        * first entry of row @p r.
+                                        *
+                                        * Note that if the given row is empty,
+                                        * i.e. does not contain any nonzero
+                                        * entries, then the iterator returned by
+                                        * this function equals
+                                        * <tt>end(r)</tt>. Note also that the
+                                        * iterator may not be dereferencable in
+                                        * that case.
+                                        */
+      const_iterator begin (const unsigned int r) const;
+
+                                       /**
+                                        * Final iterator of row <tt>r</tt>. It
+                                        * points to the first element past the
+                                        * end of line @p r, or past the end of
+                                        * the entire sparsity pattern.
+                                        *
+                                        * Note that the end iterator is not
+                                        * necessarily dereferencable. This is in
+                                        * particular the case if it is the end
+                                        * iterator for the last row of a matrix.
+                                        */
+      const_iterator end (const unsigned int r) const;
+
+                                        /**  
+                                        * Make an in-place transpose of a 
+                                        * matrix.
+                                        */
+      void transpose ();
+
+                                        /**  
+                                        * Test whether a matrix is symmetric.
+                                        * Default tolerance is zero. 
+                                        * TODO: Not implemented.
+                                        */
+      bool is_symmetric (const double tol = 0.0);
+
+                                        /** 
+                                        * Test whether a matrix is Hermitian, 
+                                        * i.e. it is the complex conjugate 
+                                        * of its transpose. 
+                                        * TODO: Not implemented.
+                                        */
+      bool is_hermitian ();
+
+                                        /** 
+                                        * Abstract Trilinos object that helps view 
+                                        * in ASCII other Trilinos objects. Currently
+                                        * this function is not implemented.
+                                        * TODO: Not implemented.
+                                        */
+      void write_ascii ();
+      
+                                        // TODO: Write an overloading
+                                        // of the operator << for output.
+                                        // Since the underlying Trilinos 
+                                        // object supports it, this should 
+                                        // be very easy.
+
+                                       /**
+                                        * Exception
+                                        */
+      DeclException1 (ExcTrilinosError,
+                      int,
+                      << "An error with error number " << arg1
+                      << " occured while calling a Trilinos function");
+
+                                       /**
+                                        * Exception
+                                        */
+      DeclException0 (ExcSourceEqualsDestination);
+
+                                       /**
+                                        * Exception
+                                        */
+      DeclException4 (ExcAccessToNonlocalElement,
+                     int, int, int, int,
+                     << "You tried to access element (" << arg1
+                     << "/" << arg2 << ")"
+                     << " of a distributed matrix, but only rows "
+                     << arg3 << " through " << arg4
+                     << " are stored locally and can be accessed.");
+
+                                       /**
+                                        * The Epetra Trilinos mapping that
+                                       * assigns parts of the matrix to
+                                       * the individual processes.
+                                       */
+      Epetra_Map map;
+
+                                       /**
+                                        * A sparse matrix object in
+                                        * Trilinos to be used for 
+                                       * finite element based problems
+                                       * which allows for assembling into
+                                       * non-local elements. 
+                                       * The actual type, a sparse
+                                        * matrix, is set in the constructor.
+                                        */
+      std::auto_ptr<Epetra_FECrsMatrix> matrix;
+
+    protected:
+                                       /**
+                                        * 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;
+  };
+
+
+
+#ifndef DOXYGEN
+// -------------------------- inline and template functions ----------------------
+
+
+  namespace MatrixIterators
+  {
+
+    inline
+    const_iterator::Accessor::
+    Accessor (const SparseMatrix   *matrix,
+              const unsigned int  row,
+              const unsigned int  index)
+                    :
+                    matrix(const_cast<SparseMatrix*>(matrix)),
+                    a_row(row),
+                    a_index(index)
+    {
+      visit_present_row ();
+    }
+
+
+    inline
+    unsigned int
+    const_iterator::Accessor::row() const
+    {
+      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
+      return a_row;
+    }
+
+
+    inline
+    unsigned int
+    const_iterator::Accessor::column() const
+    {
+      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
+      return (*colnum_cache)[a_index];
+    }
+
+
+    inline
+    unsigned int
+    const_iterator::Accessor::index() const
+    {
+      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
+      return a_index;
+    }
+
+
+    inline
+    TrilinosScalar
+    const_iterator::Accessor::value() const
+    {
+      Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
+      return (*value_cache)[a_index];
+    }
+
+
+    inline
+    const_iterator::
+    const_iterator(const SparseMatrix   *matrix,
+                   const unsigned int  row,
+                   const unsigned int  index)
+                    :
+                    accessor(matrix, row, index)
+    {}
+
+
+
+    inline
+    const_iterator &
+    const_iterator::operator++ ()
+    {
+      Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
+
+      ++accessor.a_index;
+
+                                       // if at end of line: do one step, then
+                                       // cycle until we find a row with a
+                                       // nonzero number of entries
+      if (accessor.a_index >= accessor.colnum_cache->size())
+        {
+          accessor.a_index = 0;
+          ++accessor.a_row;
+      
+          while ((accessor.a_row < accessor.matrix->m())
+                 &&
+                 (accessor.matrix->row_length(accessor.a_row) == 0))
+            ++accessor.a_row;
+
+          accessor.visit_present_row();
+        }
+      return *this;
+    }
+
+
+    inline
+    const_iterator
+    const_iterator::operator++ (int)
+    {
+      const const_iterator old_state = *this;
+      ++(*this);
+      return old_state;
+    }
+
+
+    inline
+    const const_iterator::Accessor &
+    const_iterator::operator* () const
+    {
+      return accessor;
+    }
+
+
+    inline
+    const const_iterator::Accessor *
+    const_iterator::operator-> () const
+    {
+      return &accessor;
+    }
+
+
+    inline
+    bool
+    const_iterator::
+    operator == (const const_iterator& other) const
+    {
+      return (accessor.a_row == other.accessor.a_row &&
+              accessor.a_index == other.accessor.a_index);
+    }
+
+
+    inline
+    bool
+    const_iterator::
+    operator != (const const_iterator& other) const
+    {
+      return ! (*this == other);
+    }
+
+
+    inline
+    bool
+    const_iterator::
+    operator < (const const_iterator& other) const
+    {
+      return (accessor.row() < other.accessor.row() ||
+              (accessor.row() == other.accessor.row() &&
+               accessor.index() < other.accessor.index()));
+    }
+    
+  }
+  
+  
+  inline
+  TrilinosScalar
+  SparseMatrix::operator() (const unsigned int i,
+                          const unsigned int j) const
+  {
+    return el(i,j);
+  }
+
+  
+
+  inline
+  SparseMatrix::const_iterator
+  SparseMatrix::begin() const
+  {
+    return const_iterator(this, 0, 0);
+  }
+
+
+  inline
+  SparseMatrix::const_iterator
+  SparseMatrix::end() const
+  {
+    return const_iterator(this, m(), 0);
+  }
+
+
+  inline
+  SparseMatrix::const_iterator
+  SparseMatrix::begin(const unsigned int r) const
+  {
+    Assert (r < m(), ExcIndexRange(r, 0, m()));
+    if (row_length(r) > 0)
+      return const_iterator(this, r, 0);
+    else
+      return end (r);
+  }
+
+
+  inline
+  SparseMatrix::const_iterator
+  SparseMatrix::end(const unsigned int r) const
+  {
+    Assert (r < m(), ExcIndexRange(r, 0, m()));
+
+                                     // place the iterator on the first entry
+                                     // past this line, or at the end of the
+                                     // matrix
+    for (unsigned int i=r+1; i<m(); ++i)
+      if (row_length(i) > 0)
+        return const_iterator(this, i, 0);
+    
+                                     // if there is no such line, then take the
+                                     // end iterator of the matrix
+    return end();
+  }
+
+
+  
+  inline
+  bool
+  SparseMatrix::in_local_range (const unsigned int index) const
+  {
+    int begin, end;
+    begin = matrix->RowMap().MinMyGID();
+    end = matrix->RowMap().MaxMyGID();
+    
+    return ((index >= static_cast<unsigned int>(begin)) &&
+            (index < static_cast<unsigned int>(end)));
+  }
+
+#endif // DOXYGEN      
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif // DEAL_II_USE_TRILINOS
+
+
+/*----------------------------   trilinos_matrix_base.h     ---------------------------*/
+
+#endif
+/*----------------------------   trilinos_matrix_base.h     ---------------------------*/
diff --git a/deal.II/lac/include/lac/trilinos_vector.h b/deal.II/lac/include/lac/trilinos_vector.h
new file mode 100755 (executable)
index 0000000..7998b13
--- /dev/null
@@ -0,0 +1,1013 @@
+//---------------------------------------------------------------------------
+//    $Id: Trilinos_vector_base.h 15453 2008-01-29 1:00:00Z marquez $
+//    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_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>
+
+#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"
+
+DEAL_II_NAMESPACE_OPEN
+
+                                 // forward declaration
+template <typename number> class Vector;
+
+
+/**
+ * A namespace in which wrapper classes for Trilinos objects reside.
+ *
+ * @ingroup TrilinosWrappers
+ * @ingroup Vectors
+ * @see @ref SoftwareTrilinos
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ */
+namespace TrilinosWrappers
+{
+                                   // forward declaration
+  class Vector;
+
+                                  /**
+                                   * @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 (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.
+ *
+ * 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.
+ *
+  * <h3>Parallel communication model</h3>
+ *
+ * The parallel functionality of Trilinos is built on top of the Message Passing
+ * Interface (MPI). MPI's communication model is built on collective
+ * communications: if one process wants something from another, that other
+ * process has to be willing to accept this communication. A process cannot
+ * query data from another process by calling a remote function, without that
+ * other process expecting such a transaction. The consequence is that most of
+ * the operations in the base class of this class have to be called
+ * collectively. For example, if you want to compute the l2 norm of a parallel
+ * vector, @em all processes across which this vector is shared have to call
+ * the @p l2_norm function. If you don't do this, but instead only call the @p
+ * l2_norm function on one process, then the following happens: This one
+ * process will call one of the collective MPI functions and wait for all the
+ * other processes to join in on this. Since the other processes don't call
+ * this function, you will either get a time-out on the first process, or,
+ * worse, by the time the next a callto a Trilinos function generates an MPI
+ * message on the other processes , you will get a cryptic message that only a
+ * subset of processes attempted a communication. These bugs can be very hard
+ * to figure out, unless you are well-acquainted with the communication model
+ * of MPI, and know which functions may generate MPI messages.
+ *
+ * One particular case, where an MPI message may be generated unexpectedly is
+ * discussed below.
+ *
+ * <h3>Accessing individual elements of a vector</h3>
+ *
+ * Trilinos does allow read access to individual elements of a vector, but in the
+ * distributed case only to elements that are stored locally. We implement
+ * this through calls like <tt>d=vec(i)</tt>. However, if you access an
+ * element outside the locally stored range, an exception is generated.
+ *
+ * In contrast to read access, Trilinos (and the respective deal.II wrapper
+ * classes) allow to write (or add) to individual elements of vectors, even if
+ * they are stored on a different process. You can do this writing, for
+ * example, <tt>vec(i)=d</tt> or <tt>vec(i)+=d</tt>, or similar
+ * operations. There is one catch, however, that may lead to very confusing
+ * error messages: Trilinos requires application programs to call the compress()
+ * function when they switch from adding, to elements to writing to
+ * elements. The reasoning is that all processes might accumulate addition
+ * operations to elements, even if multiple processes write to the same
+ * elements. By the time we call compress() the next time, all these additions
+ * are executed. However, if one process adds to an element, and another
+ * overwrites to it, the order of execution would yield non-deterministic
+ * behavior if we don't make sure that a synchronisation with compress()
+ * happens in between.
+ *
+ * In order to make sure these calls to compress() happen at the appropriate
+ * time, the deal.II wrappers keep a state variable that store which is the
+ * presently allowed operation: additions or writes. If it encounters an
+ * operation of the opposite kind, it calls compress() and flips the
+ * state. This can sometimes lead to very confusing behavior, in code that may
+ * for example look like this:
+ * @verbatim
+ *   TrilinosWrappers::Vector vector;
+ *   ...
+ *                   // do some write operations on the vector
+ *   for (unsigned int i=0; i<vector->size(); ++i)
+ *     vector(i) = i;
+ *
+ *                   // do some additions to vector elements, but
+ *                   // only for some elements
+ *   for (unsigned int i=0; i<vector->size(); ++i)
+ *     if (some_condition(i) == true)
+ *       vector(i) += 1;
+ *
+ *                   // do another collective operation
+ *   const double norm = vector->l2_norm();
+ * @endverbatim
+ *
+ * This code can run into trouble: by the time we see the first addition
+ * operation, we need to flush the overwrite buffers for the vector, and the
+ * deal.II library will do so by calling compress(). However, it will only do
+ * so for all processes that actually do an addition -- if the condition is
+ * never true for one of the processes, then this one will not get to the
+ * actual compress() call, whereas all the other ones do. This gets us into
+ * trouble, since all the other processes hang in the call to flush the write
+ * buffers, while the one other process advances to the call to compute the l2
+ * norm. At this time, you will get an error that some operation was attempted
+ * by only a subset of processes. This behavior may seem surprising, unless
+ * you know that write/addition operations on single elements may trigger this
+ * behavior.
+ *
+ * The problem described here may be avoided by placing additional calls to
+ * compress(), or making sure that all processes do the same type of
+ * operations at the same time, for example by placing zero additions if
+ * necessary.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ */
+  class Vector
+  {
+    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. It doesn't do
+                                        * anything, derived classes will have
+                                        * to initialize the data.
+                                        */
+      Vector ();
+                                       /**
+                                        * One of the constructors that
+                                       * actually build a vector. This
+                                       * one requires prior knowledge
+                                        * of the size of the vector and
+                                       * a communicator from 
+                                        * Epetra_CommSerial or Epetra_CommMpi,
+                                       * depending on whether we use a 
+                                       * serial or parallel MPI-based program.
+                                       * This command distributes the
+                                       * vector linearly among the processes,
+                                       * from the beginning to the end, 
+                                       * so you might want to use some
+                                       * more advanced mapping and the
+                                       * third constructor with argument
+                                       * Epetra_Map.
+                                        */
+      Vector (unsigned int GlobalSize, Epetra_Comm &Comm);
+
+                                       /**
+                                       * Third constructor. It takes an
+                                       * Epetra_Map that already knows how
+                                       * to distribute the individual 
+                                       * components among the MPI processors.
+                                        */
+      Vector (const Epetra_Map &InputMap);
+
+                                       /**
+                                        * Copy constructor. Sets the dimension
+                                        * to that of the given vector and uses
+                                       * the map of that vector, and
+                                        * copies all elements.
+                                        */
+      Vector (const Vector &v);
+
+                                       /**
+                                        * Destructor
+                                        */
+      virtual ~Vector ();
+
+                                      /** 
+                                       * Reinit functionality. This function
+                                       * destroys the old vector content 
+                                       * and generates a new one based on
+                                       * the input map.
+                                       */
+      void reinit (const Epetra_Map &input_map);
+
+                                       /**
+                                        * 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);
+
+                                       /**
+                                        * 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
+                                       * function 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);
+
+      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);
+
+                                       /**
+                                        * Scaling and multiple addition.
+                                        * <tt>*this = s*(*this) + a*V + b*W + c*X</tt>.
+                                        */
+      void sadd (const TrilinosScalar s,
+                 const TrilinosScalar a,
+                 const Vector        &V,
+                 const TrilinosScalar b,
+                 const Vector        &W,
+                 const TrilinosScalar c,
+                 const Vector        &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 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);
+
+                                       /**
+                                        * 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 (Vector &v);
+
+                                      /**
+                                       * 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.");
+
+
+                                       /**
+                                        * An Epetra map used to map vector data
+                                        * accross multiple processes. This is
+                                        * the communicator and data distribution
+                                       * object common to all
+                                        * Trilinos objects used by deal.II.
+                                       * TODO: we probably only need a pointer
+                                       * to the map, since the information
+                                       * is provided from outside and there
+                                       * is no need to copy the map. 
+                                       * Especially not when we have many
+                                       * vectors based on the same map.
+                                        */
+      Epetra_Map map;
+
+                                       /**
+                                        * An Epetra distibuted vector type.
+                                        * Requires an existing Epetra_Map for
+                                        * storing data.
+                                        */
+      std::auto_ptr<Epetra_FEVector> vector;
+
+
+                                       /**
+                                        * 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;
+
+                                       /**
+                                        * Make the reference class a friend.
+                                        */
+      friend class internal::VectorReference;
+
+  };
+
+
+
+// ------------------- inline and template functions --------------
+
+/**
+ * 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 Wolfgang Bangerth, 2004
+ */
+  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
+    {
+                                       // 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
+  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
+  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
+
+#endif // DEAL_II_USE_TRILINOS
+
+/*----------------------------   trilinos_vector_base.h     ---------------------------*/
+
+#endif
+/*----------------------------   trilinos_vector_base.h     ---------------------------*/
diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc
new file mode 100755 (executable)
index 0000000..5ce2e21
--- /dev/null
@@ -0,0 +1,687 @@
+//---------------------------------------------------------------------------
+//    $Id: trilinos_matrix_base.cc 16014 2008-04-30 16:18:03Z young $
+//    Version: $Name$
+//
+//    Copyright (C) 2004, 2005, 2006, 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.h>
+#include <lac/trilinos_sparse_matrix.h>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+  namespace MatrixIterators
+  {
+    void
+    SparseMatrix::const_iterator::Accessor::
+    visit_present_row ()
+    {
+                                       // if we are asked to visit the
+                                       // past-the-end line, then simply
+                                       // release all our caches and go on
+                                       // with life
+      if (this->a_row == matrix->m())
+        {
+          colnum_cache.reset ();
+          value_cache.reset ();
+
+          return;
+        }
+      
+                                       // otherwise first flush Trilinos caches
+      matrix->compress ();
+
+                                       // get a representation of the present
+                                       // row
+      int ncols;
+      int colnums = matrix->n();
+      TrilinosScalar *values = new TrilinosScalar(colnums);
+      
+      int ierr;
+      ierr = matrix->matrix->ExtractGlobalRowCopy((int)this->a_row, colnums,  
+                                                 ncols, &(values[0]));
+      AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+                                       // copy it into our caches if the line
+                                       // isn't empty. if it is, then we've
+                                       // done something wrong, since we
+                                       // shouldn't have initialized an
+                                       // iterator for an empty line (what
+                                       // would it point to?)
+      Assert (ncols != 0, ExcInternalError());
+      colnum_cache.reset (new std::vector<unsigned int> (colnums,
+                                                         colnums+ncols));
+      value_cache.reset (new std::vector<TrilinosScalar> (values, values+ncols));
+    }
+  }
+
+
+                                     // The constructor is actually the
+                                     // only point where we have to check
+                                      // whether we build a serial or
+                                      // a parallel Trilinos matrix.
+                                      // In the end, it does not even
+                                      // matter how many threads there
+                                      // are, but only if we use an
+                                      // MPI compiler or a standard 
+                                     // compiler. So, one thread on
+                                      // an MPI compiler will still get
+                                      // a parallel interface.
+  SparseMatrix::SparseMatrix ()
+                  :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+                 map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
+#else
+                 map (0,0,Epetra_SerialComm()),
+#endif
+                 matrix (std::auto_ptr<Epetra_FECrsMatrix>
+                                (new Epetra_FECrsMatrix(Copy, map, 0))),
+                  last_action (Insert)
+  {}
+
+
+  SparseMatrix::SparseMatrix (const Epetra_Map  &InputMap,
+                             const unsigned int n_max_entries_per_row)
+                  :
+                 map (InputMap),
+                 matrix (std::auto_ptr<Epetra_FECrsMatrix>
+                                (new Epetra_FECrsMatrix(Copy, map, 
+                                       int(n_max_entries_per_row), false))),
+                  last_action (Insert)
+  {}
+
+  SparseMatrix::SparseMatrix (const Epetra_Map  &InputMap,
+                             const std::vector<unsigned int> &n_entries_per_row)
+                  :
+                 map (InputMap),
+                 matrix (std::auto_ptr<Epetra_FECrsMatrix>
+                   (new Epetra_FECrsMatrix(Copy, map, 
+                     (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
+                     true))),
+                  last_action (Insert)
+  {}
+
+
+  SparseMatrix::~SparseMatrix ()
+  {
+  }
+
+
+  
+  void
+  SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern,
+                       const unsigned int               n_max_entries_per_row)
+  {
+
+    unsigned int n_rows = sparsity_pattern.n_rows();
+    
+    Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
+           ExcDimensionMismatch (matrix->NumGlobalRows(),
+                                 sparsity_pattern.n_rows()));
+           
+    std::vector<double> values(n_max_entries_per_row, 0.);
+    std::vector<int>    row_indices(n_max_entries_per_row);
+    
+    for (unsigned int row=0; row<n_rows; ++row)
+      {
+       const int row_length = sparsity_pattern.row_length(row);
+       row_indices.resize (row_length, 0);
+        values.resize (row_length, 0.);
+       
+       for (int col=0; col< row_length; ++col)
+         row_indices[col] = sparsity_pattern.column_number (row, col);
+       
+       matrix->InsertGlobalValues(row, row_length,
+                                   &values[0], &row_indices[0]);
+      }
+
+                                 // In the end, the matrix is to
+                                 // be compressed in order to be
+                                 // really ready. However, that is
+                                 // a collective operation, so it
+                                 // has to be called on all processes
+                                 // by the user, whereas this function
+                                 // should only be used on one processor
+                                 // since our sparsity pattern data
+                                 // types are all serial.
+  }
+
+  
+  
+  void
+  SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern)
+  {
+    unsigned int n_rows = sparsity_pattern.n_rows();
+    
+    Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
+           ExcDimensionMismatch (matrix->NumGlobalRows(),
+                                 sparsity_pattern.n_rows()));
+
+    std::vector<int> n_entries_per_row(n_rows);
+    
+    for (unsigned int row=0; row<n_rows; ++row)
+      n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
+    
+    const unsigned int n_max_entries_per_row = *std::max_element (
+                   &n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
+
+    reinit (sparsity_pattern, n_max_entries_per_row);
+  }
+  
+  
+  
+  void
+  SparseMatrix::reinit (const Epetra_Map                &input_map,
+                       const CompressedSparsityPattern &sparsity_pattern)
+  {
+
+    unsigned int n_rows = sparsity_pattern.n_rows();
+    matrix.reset();
+    map = input_map;
+    Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
+           ExcDimensionMismatch (input_map.NumGlobalElements(),
+                                 sparsity_pattern.n_rows()));
+    
+    std::vector<int> n_entries_per_row(n_rows);
+    
+    for (unsigned int row=0; row<n_rows; ++row)
+      n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
+    
+    matrix = std::auto_ptr<Epetra_FECrsMatrix>
+             (new Epetra_FECrsMatrix(Copy, map, &n_entries_per_row[0], true));
+    
+    const unsigned int n_max_entries_per_row = *std::max_element (
+                   &n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
+
+    reinit (sparsity_pattern, n_max_entries_per_row);
+  }
+
+
+
+  void
+  SparseMatrix::clear ()
+  {
+                                     // When we clear the matrix,
+                                    // reset the pointer and 
+                                    // generate an empty matrix.
+    matrix.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
+
+    matrix = std::auto_ptr<Epetra_FECrsMatrix> 
+             (new Epetra_FECrsMatrix(Copy, map, 0));
+  }
+
+
+
+  void
+  SparseMatrix::compress ()
+  {
+                                     // flush buffers
+    int ierr;
+    ierr = matrix->GlobalAssemble (true);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    ierr = matrix->OptimizeStorage ();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  SparseMatrix &
+  SparseMatrix::operator = (const double d)
+  {
+    Assert (d==0, ExcScalarAssignmentOnlyForZeroValue());
+
+    compress ();
+
+    const int ierr = matrix->PutScalar(d);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+
+  void
+  SparseMatrix::set (const unsigned int   i,
+                    const unsigned int   j,
+                    const TrilinosScalar value)
+  {
+
+    Assert (numbers::is_finite(value),
+           ExcMessage("The given value is not finite but either "
+                      "infinite or Not A Number (NaN)"));
+
+    if (last_action == Add)
+      {
+        int ierr;
+        ierr = matrix->GlobalAssemble(false);
+        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+       last_action = Insert;
+      }
+
+    int trilinos_i = i;
+    int trilinos_j = j;
+
+    const int ierr = matrix->ReplaceGlobalValues (trilinos_i, 1,
+                                                 const_cast<double*>(&value), 
+                                                 &trilinos_j);
+
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  void
+  SparseMatrix::add (const unsigned int   i,
+                    const unsigned int   j,
+                    const TrilinosScalar value)
+  {
+
+    Assert (numbers::is_finite(value), 
+           ExcMessage("The given value is not finite but either "
+                      "infinite or Not A Number (NaN)"));
+
+    if (last_action == Insert)
+      {
+        int ierr;
+        ierr = matrix->GlobalAssemble(false);
+        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+        last_action = Add;
+      }
+
+                                     // we have to do above actions in any
+                                     // case to be consistent with the MPI
+                                     // communication model (see the
+                                     // comments in the documentation of
+                                     // TrilinosWrappers::MPI::Vector), but we
+                                     // can save some work if the addend is
+                                     // zero
+    if (value == 0)
+      return;
+
+    int trilinos_i = i;
+    int trilinos_j = j;
+
+    const int ierr = matrix->SumIntoGlobalValues (trilinos_i, 1, 
+                                                 const_cast<double*>(&value), 
+                                                 &trilinos_j);
+
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  TrilinosScalar
+  SparseMatrix::el (const unsigned int i,
+                   const unsigned int j) const
+  {
+                                      // Extract local indices in
+                                      // the matrix.
+    int trilinos_i = matrix->LRID(i), trilinos_j = matrix->LRID(j);
+    TrilinosScalar value = 0.;
+    
+                                      // If the data is not on the
+                                      // present processor, we can't
+                                      // continue.
+    if ((trilinos_i == -1 ) || (trilinos_j == -1))
+      {
+        Assert (false, ExcAccessToNonlocalElement(i, j, local_range().first,
+                                                 local_range().second));
+      }
+    else
+    {
+                                     // Check whether the matrix 
+                                     // already is transformed to
+                                     // local indices.
+      if (!matrix->Filled())
+       matrix->FillComplete(true);
+
+                                     // Prepare pointers for extraction
+                                     // of a view of the row.
+      int nnz_present = matrix->NumMyEntries(trilinos_i);
+      int nnz_extracted;
+      int *col_indices;
+      TrilinosScalar *values;
+
+                                     // Generate the view and make
+                                     // sure that we have not generated
+                                     // an error.
+      int ierr = matrix->ExtractMyRowView(trilinos_i, nnz_extracted,
+                                         values, col_indices);
+      Assert (ierr==0, ExcTrilinosError(ierr));
+
+      Assert (nnz_present == nnz_extracted,
+             ExcDimensionMismatch(nnz_present, nnz_extracted));
+
+                                     // Search the index where we
+                                     // look for the value, and then
+                                     // finally get it.
+      int* index = std::find(&col_indices[0],&col_indices[0] + nnz_present,
+                            trilinos_j);
+
+      int position;
+      if (!index)
+       value = 0;
+      else
+       {
+         position = (int)(index - &(col_indices[0]));
+         value = values[position];
+       }
+    }
+
+    return value;
+  }
+
+
+
+  TrilinosScalar
+  SparseMatrix::diag_element (const unsigned int i) const
+  {
+    Assert (m() == n(), ExcNotQuadratic());
+
+                                     // this doesn't seem to work any
+                                     // different than any other element
+    return el(i,i);
+  }
+
+
+
+  unsigned int
+  SparseMatrix::m () const
+  {
+    int n_rows = matrix->NumGlobalRows();
+
+    return n_rows;
+  }
+
+
+
+  unsigned int
+  SparseMatrix::n () const
+  {
+    int n_cols = matrix -> NumGlobalCols();
+    return n_cols;
+  }
+
+
+
+  unsigned int
+  SparseMatrix::local_size () const
+  {
+    int n_rows = matrix -> NumMyRows();
+
+    return n_rows;
+  }
+
+
+
+  std::pair<unsigned int, unsigned int>
+  SparseMatrix::local_range () const
+  {
+    int begin, end;
+    begin = matrix->RowMap().MinMyGID();
+    end = matrix->RowMap().MaxMyGID();
+    
+    return std::make_pair (begin, end);
+  }
+
+
+
+  unsigned int
+  SparseMatrix::n_nonzero_elements () const
+  {
+    int nnz = matrix->NumGlobalNonzeros();
+
+    return static_cast<unsigned int>(nnz);
+  }
+
+
+
+  unsigned int
+  SparseMatrix::row_length (const unsigned int row) const
+  {
+    Assert (row < m(), ExcInternalError());
+
+                                // get a representation of the present
+                                // row
+    int ncols = -1;
+    int local_row = matrix->RowMap().LID(row);
+
+                                // on the processor who owns this
+                                // row, we'll have a non-negative
+                                // value.
+    if (local_row >= 0)
+      {
+       int ierr = matrix->NumMyRowEntries (local_row, ncols);
+       AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+      }
+
+    return ncols;
+  }
+
+
+  TrilinosScalar
+  SparseMatrix::l1_norm () const
+  {
+    if (!matrix->Filled())
+      matrix->FillComplete();
+
+    TrilinosScalar result = matrix->NormOne();
+
+    return result;
+  }
+  
+  
+
+  TrilinosScalar
+  SparseMatrix::linfty_norm () const
+  {
+    if (!matrix->Filled())
+      matrix->FillComplete();
+
+    TrilinosScalar result = matrix->NormInf();
+
+    return result;
+  }
+
+
+
+  TrilinosScalar
+  SparseMatrix::frobenius_norm () const
+  {
+    if (!matrix->Filled())
+      matrix->FillComplete();
+
+    TrilinosScalar result = matrix->NormFrobenius();
+
+    return result;
+  }
+
+
+
+  SparseMatrix &
+  SparseMatrix::operator *= (const TrilinosScalar a)
+  {
+    const int ierr = matrix->Scale (a);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+
+  SparseMatrix &
+  SparseMatrix::operator /= (const TrilinosScalar a)
+  {
+    const TrilinosScalar factor = 1./a;
+
+    const int ierr = matrix->Scale (factor);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+  void
+  SparseMatrix::vmult (Vector       &dst,
+                      const Vector &src) const
+  {
+    Assert (&src != &dst, ExcSourceEqualsDestination());
+
+    if (!matrix->Filled())
+      matrix->FillComplete();
+
+    const int ierr = matrix->Multiply (false, *(src.vector), *(dst.vector));
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  void
+  SparseMatrix::Tvmult (Vector       &dst,
+                       const Vector &src) const
+  {
+    Assert (&src != &dst, ExcSourceEqualsDestination());
+
+    if (!matrix->Filled())
+      matrix->FillComplete();
+
+    const int ierr = matrix->Multiply (true, *(src.vector), *(dst.vector));
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  void
+  SparseMatrix::vmult_add (Vector       &dst,
+                          const Vector &src) const
+  {
+    Assert (&src != &dst, ExcSourceEqualsDestination());
+
+    Vector tmp = dst;
+    vmult (dst, src);
+    dst += tmp;
+  }
+
+
+
+  void
+  SparseMatrix::Tvmult_add (Vector       &dst,
+                           const Vector &src) const
+  {
+    Assert (&src != &dst, ExcSourceEqualsDestination());
+
+    Vector tmp = dst;
+    vmult (dst, src);
+    dst += tmp;
+  }
+
+
+
+  TrilinosScalar
+  SparseMatrix::matrix_norm_square (const Vector &v) const
+  {
+    Vector tmp(v.map);
+    vmult (tmp, v);
+    return tmp*v;
+  }
+
+
+
+  TrilinosScalar
+  SparseMatrix::matrix_scalar_product (const Vector &u,
+                                      const Vector &v) const
+  {
+    Vector tmp(v.map);
+    vmult (tmp, v);
+    return u*tmp;
+  }
+
+
+
+  TrilinosScalar
+  SparseMatrix::residual (Vector       &dst,
+                         const Vector &x,
+                         const Vector &b) const
+  {
+    vmult (dst, x);
+    dst -= b;
+    dst *= -1;
+
+    return dst.l2_norm();
+  }
+
+
+
+                                 // TODO: Currently this only flips
+                                  // a flag that tells Trilinos that
+                                  // any application should be done with
+                                  // the transpose. However, the matrix
+                                  // structure is not reset.
+  void
+  SparseMatrix::transpose () 
+  {
+    int ierr;
+
+    if (!matrix->UseTranspose())
+      {
+        ierr = matrix->SetUseTranspose (true);
+        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+      }
+    else
+      {
+        ierr = matrix->SetUseTranspose (false);
+        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+      }
+  }
+
+
+
+  bool
+  SparseMatrix::is_symmetric (const double tolerance) 
+  {
+    //bool truth;
+    if (tolerance == 0)
+      Assert (false, ExcNotImplemented());
+
+    return false;
+  }  
+
+
+
+  bool
+  SparseMatrix::is_hermitian () 
+  {
+    //bool truth;
+
+    Assert (false, ExcNotImplemented());
+    return false;
+  }  
+
+  void
+  SparseMatrix::write_ascii ()
+  {
+    Assert (false, ExcNotImplemented());
+  }
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS
diff --git a/deal.II/lac/source/trilinos_vector.cc b/deal.II/lac/source/trilinos_vector.cc
new file mode 100755 (executable)
index 0000000..18f8afc
--- /dev/null
@@ -0,0 +1,786 @@
+//---------------------------------------------------------------------------
+//    $Id: trilinos_vector_base.cc 15631 2008-01-29 01:00:00Z marquez $
+//    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.h>
+
+#include <cmath>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+  namespace internal
+  {
+    VectorReference::operator TrilinosScalar () const // i believe useless with
+                                                      // trilinos
+    {
+      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()-1));
+      return *(*(vector.vector))[index];
+    }
+  }
+
+  Vector::Vector ()
+                  :
+#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_FEVector> 
+                        (new Epetra_FEVector(map))),
+                  last_action (Insert)
+  {}
+
+  Vector::Vector (unsigned int GlobalSize, Epetra_Comm &Comm)
+                  :
+                 map (GlobalSize, 0, Comm),
+                 vector (std::auto_ptr<Epetra_FEVector> 
+                         (new Epetra_FEVector(map))),
+                  last_action (Insert)
+  {}
+
+  
+  Vector::Vector (const Epetra_Map &InputMap)
+                  :
+                 map (InputMap),
+                 vector (std::auto_ptr<Epetra_FEVector> 
+                         (new Epetra_FEVector(map))),
+                  last_action (Insert)
+  {}
+  
+  
+  Vector::Vector (const Vector &v)
+                  :
+                 map (v.map),
+                 vector(std::auto_ptr<Epetra_FEVector> 
+                        (new Epetra_FEVector(*(v.vector)))),
+                  last_action (Insert)
+  {}
+  
+
+
+  Vector::~Vector ()
+  {}
+
+
+
+  void
+  Vector::reinit (const Epetra_Map &input_map)
+  {
+    vector.reset();
+    map = input_map;
+    
+    vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(input_map));
+    last_action = Insert;
+  }
+
+
+
+  void
+  Vector::clear ()
+  {
+                                     // When we clear the matrix,
+                                    // reset the pointer and 
+                                    // generate an empty matrix.
+    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));
+  }
+
+
+
+  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));
+  }
+
+
+
+  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)"));
+
+    const int ierr = vector->PutScalar(s);
+
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    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))[i]!=(*vector)[i]) return false;
+    return true;
+  }
+
+
+
+  bool
+  Vector::operator != (const Vector &v) const
+  {
+    Assert (size() == v.size(),
+            ExcDimensionMismatch(size(), v.size()));
+
+    return (!(*this==v));
+  }
+
+
+
+  unsigned int
+  Vector::size () const
+  {
+    return (unsigned int) vector->Map().NumGlobalElements();
+  }
+
+
+
+  unsigned int
+  Vector::local_size () const
+  {
+    return (unsigned int) vector->Map().NumMyElements();
+  }
+
+
+
+  std::pair<unsigned int, unsigned int>
+  Vector::local_range () const
+  {
+    int begin, end;
+    begin = vector->Map().MinMyGID();
+    end = vector->Map().MaxMyGID();
+    return std::make_pair (begin, end);
+  }
+
+
+
+  TrilinosScalar
+  Vector::el (const unsigned int index) const
+  {
+                                      // Extract local indices in
+                                      // the vector.
+    int trilinos_i = map.LID(index);
+    TrilinosScalar value = 0.;
+    if (trilinos_i == -1 )
+      {
+        Assert (false, ExcAccessToNonlocalElement(index, local_range().first,
+                                                 local_range().second));
+      }
+    else
+      value = (*vector)[0][trilinos_i];
+    
+    return value;
+  }
+
+
+  void
+  Vector::set (const std::vector<unsigned int>    &indices,
+              const std::vector<TrilinosScalar>  &values)
+  {
+    Assert (indices.size() == values.size(),
+           ExcMessage ("Function called with arguments of different sizes"));
+    
+    set (indices.size(), &indices[0], &values[0]);
+  }
+
+  
+  void
+  Vector::set (const unsigned int    n_elements,
+              const unsigned int   *indices,
+              const TrilinosScalar *values)
+  {
+    if (last_action == Add)
+      {
+       vector->GlobalAssemble(Add);
+       last_action = Insert;
+      }
+
+    const int ierr= vector->ReplaceGlobalValues (n_elements, 
+                                     (int*)(const_cast<unsigned int*>(indices)), 
+                                     const_cast<TrilinosScalar*>(values));
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+  
+  void
+  Vector::add (const std::vector<unsigned int>    &indices,
+              const std::vector<TrilinosScalar>  &values)
+  {
+    Assert (indices.size() == values.size(),
+           ExcMessage ("Function called with arguments of different sizes"));
+    
+    add (indices.size(), &indices[0], &values[0]);
+  }
+
+  
+  
+  void
+  Vector::add (const unsigned int    n_elements,
+              const unsigned int   *indices,
+              const TrilinosScalar *values)
+  {
+    if (last_action == Insert)
+      {
+       vector->GlobalAssemble(Insert);
+       last_action = Add;
+      }
+
+    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;
+    int ierr;
+
+    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));
+
+    return d;
+  }
+
+
+
+  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)
+      {
+        if (*ptr != 0)
+          {
+            flag = false;
+            break;
+          }
+        ++ptr;
+      }
+
+    return flag;
+  }
+
+
+
+  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;
+      }
+
+    return flag;
+  }
+
+
+
+  Vector &
+  Vector::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);
+    Assert (ierr == 0, ExcTrilinosError(ierr));
+
+    return *this;
+  }
+
+
+
+  Vector &
+  Vector::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);
+    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));
+    return *this;
+  }
+
+
+
+  Vector &
+  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)"));
+
+    std::vector<TrilinosScalar> list (size(), s);
+
+    int* index = new int[size()];
+    for (unsigned int i=0; i<size(); i++)
+      index[i]=i;
+    
+    const int ierr = vector->SumIntoGlobalValues(size(), index, &list[0]);
+
+    delete[] index;
+    
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  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.0);
+
+    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.0, *(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.0);
+    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)"));
+
+    Assert (size() == v.size(),
+            ExcDimensionMismatch (size(), v.size()));
+
+    const int ierr = vector->Update(a, *(v.vector), 0.0);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
+
+  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 (size() == v.size(),
+            ExcDimensionMismatch (size(), v.size()));
+
+    Assert (size() == w.size(),
+            ExcDimensionMismatch (size(), w.size()));
+
+    const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 0.0);
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+  }
+
+
+
+  void
+  Vector::ratio (const Vector &v,
+                const Vector &w)
+  {
+    Assert (size() == v.size(),
+            ExcDimensionMismatch (size(), v.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));
+  }
+
+
+
+  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!
+    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.
+    std::auto_ptr<Epetra_FEVector> tmp;
+    tmp = v.vector;
+    v.vector = vector;
+    vector = tmp;
+  }
+
+}
+
+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.