]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Trilinos vector and sparse matrix classes now only store a pointer to a user-provided...
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 8 Sep 2008 14:47:42 +0000 (14:47 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 8 Sep 2008 14:47:42 +0000 (14:47 +0000)
git-svn-id: https://svn.dealii.org/trunk@16764 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_precondition_amg.h
deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/include/lac/trilinos_vector.h
deal.II/lac/source/trilinos_precondition_amg.cc
deal.II/lac/source/trilinos_sparse_matrix.cc
deal.II/lac/source/trilinos_vector.cc

index 17f38b9b599b5cc660e6c491f363810e8dd4da3c..0b953b9011fdfcb71531da751e851d7598d99e04 100755 (executable)
@@ -54,10 +54,10 @@ namespace TrilinosWrappers
  * This class implements an algebraic multigrid (AMG) preconditioner
  * based on the Trilinos ML implementation.  What this class does is
  * twofold.  When the initialize() function is invoked, a ML
- * preconditioner object is created based on the matrix
- * that we want the preconditioner to be based on. A call of the
- * respective <code>vmult</code> function does call the respective
- * operation in the Trilinos package, where it is called
+ * preconditioner object is created based on the matrix that we want
+ * the preconditioner to be based on. A call of the respective
+ * <code>vmult</code> function does call the respective operation in
+ * the Trilinos package, where it is called
  * <code>ApplyInverse</code>. Use of this class is explained in the
  * @ref step_31 "step-31" tutorial program.
  *
@@ -71,11 +71,12 @@ namespace TrilinosWrappers
  * fact that some of the entries in the preconditioner matrix are zero
  * and hence can be neglected.
  *
- * The implementation is able to distinguish between matrices from 
- * elliptic problems and convection dominated problems. We use the standard
- * options provided by Trilinos ML for elliptic problems, except that we use a 
- * Chebyshev smoother instead of a symmetric Gauss-Seidel smoother. 
- * For most elliptic problems, Chebyshev is better than Gauss-Seidel (SSOR).
+ * The implementation is able to distinguish between matrices from
+ * elliptic problems and convection dominated problems. We use the
+ * standard options provided by Trilinos ML for elliptic problems,
+ * except that we use a Chebyshev smoother instead of a symmetric
+ * Gauss-Seidel smoother.  For most elliptic problems, Chebyshev is
+ * better than Gauss-Seidel (SSOR).
  *
  * @ingroup TrilinosWrappers
  * @ingroup Preconditioners
@@ -99,9 +100,9 @@ namespace TrilinosWrappers
                                        * multilevel hierarchy for the
                                        * solution of a linear system
                                        * with the given matrix. The
-                                       * function uses the matrix 
-                                       * format specified in
-                                       * TrilinosWrappers::SparseMatrix.
+                                       * function uses the matrix
+                                       * format specified in
+                                       * TrilinosWrappers::SparseMatrix.
                                        */
       void initialize (const SparseMatrix                    &matrix,
                       const bool                             elliptic = true,
@@ -115,11 +116,12 @@ namespace TrilinosWrappers
                                        * multilevel hierarchy for the
                                        * solution of a linear system
                                        * with the given matrix. This
-                                       * function takes a deal.ii matrix
-                                       * and copies the content into a
-                                       * Trilinos matrix, so the function
-                                       * can be considered rather 
-                                       * inefficient.
+                                       * function takes a deal.ii
+                                       * matrix and copies the
+                                       * content into a Trilinos
+                                       * matrix, so the function can
+                                       * be considered rather
+                                       * inefficient.
                                        */
       void initialize (const ::dealii::SparseMatrix<double>  &deal_ii_sparse_matrix,
                       const bool                             elliptic = true,
@@ -129,26 +131,29 @@ namespace TrilinosWrappers
                       const bool                            output_details = false);
 
                                       /**
-                                       * This function can be used 
-                                       * for a faster recalculation of
-                                       * the preconditioner construction 
-                                       * when the matrix entries 
-                                       * underlying the preconditioner 
-                                       * have changed, 
-                                       * but the matrix sparsity pattern
-                                       * has remained the same. What this
-                                       * function does is to take the 
-                                       * already generated coarsening 
-                                       * structure, compute the AMG 
-                                       * prolongation and restriction 
-                                       * according to a smoothed aggregation
-                                       * strategy and then builds the whole
-                                       * multilevel hiearchy. This function
-                                       * can be considerably faster in that
-                                       * case, since the coarsening pattern
-                                       * is usually the most difficult thing
-                                       * to do when setting up the
-                                       * AMG ML preconditioner.
+                                       * This function can be used
+                                       * for a faster recalculation
+                                       * of the preconditioner
+                                       * construction when the matrix
+                                       * entries underlying the
+                                       * preconditioner have changed,
+                                       * but the matrix sparsity
+                                       * pattern has remained the
+                                       * same. What this function
+                                       * does is to take the already
+                                       * generated coarsening
+                                       * structure, compute the AMG
+                                       * prolongation and restriction
+                                       * according to a smoothed
+                                       * aggregation strategy and
+                                       * then builds the whole
+                                       * multilevel hiearchy. This
+                                       * function can be considerably
+                                       * faster in that case, since
+                                       * the coarsening pattern is
+                                       * usually the most difficult
+                                       * thing to do when setting up
+                                       * the AMG ML preconditioner.
                                        */
       void reinit ();
 
@@ -178,11 +183,23 @@ namespace TrilinosWrappers
                                        */
       boost::shared_ptr<ML_Epetra::MultiLevelPreconditioner> multigrid_operator;
 
+                                       /**
+                                       * Internal communication
+                                       * pattern in case the matrix
+                                       * needs to be copied from
+                                       * deal.II format.
+                                       */
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
       Epetra_MpiComm     communicator;
 #else
       Epetra_SerialComm  communicator;
 #endif
+
+                                       /**
+                                       * Internal Trilinos map in
+                                       * case the matrix needs to be
+                                       * copied from deal.II format.
+                                       */
       boost::shared_ptr<Epetra_Map>       Map;
 
                                       /**
index 5cda675d0e5916fd82c26523e3d0a44e2a6ca493..2399203e18ec9bcf69d2fa57889abf5551c4b7d5 100755 (executable)
@@ -51,14 +51,14 @@ namespace TrilinosWrappers
   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.
+ * 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.
+ * 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
@@ -66,55 +66,53 @@ namespace TrilinosWrappers
     class const_iterator
     {
       private:
-                                         /**
-                                          * Accessor class for iterators
-                                          */
+                                       /**
+                                       * Accessor class for iterators
+                                       */
         class Accessor
         {
           public:
-                                             /**
-                                              * Constructor. Since we use
-                                              * accessors only for read
-                                              * access, a const matrix
-                                              * pointer is sufficient.
-                                              */
+                                       /**
+                                       * 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.
-                                              */
+                                       /**
+                                       * Row number of the element
+                                       * represented by this object.
+                                       */
             unsigned int row() const;
 
-                                             /**
-                                              * Index in row of the element
-                                              * represented by this
-                                              * object.
-                                              */
+                                       /**
+                                       * Index in row of the element
+                                       * represented by this object.
+                                       */
             unsigned int index() const;
 
-                                             /**
-                                              * Column number of the
-                                              * element represented by
-                                              * this object.
-                                              */
+                                       /**
+                                       * Column number of the element
+                                       * represented by this object.
+                                       */
             unsigned int column() const;
 
-                                             /**
-                                              * Value of this matrix entry.
-                                              */
+                                      /**
+                                       * Value of this matrix entry.
+                                       */
             TrilinosScalar value() const;
 
-                                             /**
-                                              * Exception
-                                              */
+                                       /**
+                                       * Exception
+                                       */
             DeclException0 (ExcBeyondEndOfMatrix);
-                                             /**
-                                              * Exception
-                                              */
+
+                                      /**
+                                       * Exception
+                                       */
             DeclException3 (ExcAccessToNonlocalRow,
                             int, int, int,
                             << "You tried to access row " << arg1
@@ -123,126 +121,128 @@ namespace TrilinosWrappers
                             << " are stored locally and can be accessed.");
             
           private:
-                                             /**
-                                              * The matrix accessed.
-                                              */
+                                       /**
+                                       * The matrix accessed.
+                                       */
             mutable SparseMatrix *matrix;
 
-                                             /**
-                                              * Current row number.
-                                              */
+                                       /**
+                                       * Current row number.
+                                       */
             unsigned int a_row;
 
-                                             /**
-                                              * Current index in 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.
-                                              */
+                                       /**
+                                       * 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.
-                                              */
+                                       /**
+                                       * 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.
-                                              */
+                                      /**
+                                       * 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.
-                                              */
+                                       /**
+                                       * 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.
-                                          */ 
+                                       /**
+                                       * 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.
-                                          */
+                                       /**
+                                       * Prefix increment.
+                                       */
         const_iterator& operator++ ();
 
-                                         /**
-                                          * Postfix increment.
-                                          */
+                                       /**
+                                       * Postfix increment.
+                                       */
         const_iterator operator++ (int);
 
-                                         /**
-                                          * Dereferencing operator.
-                                          */
+                                       /**
+                                       * Dereferencing operator.
+                                       */
         const Accessor& operator* () const;
 
-                                         /**
-                                          * Dereferencing operator.
-                                          */
+                                       /**
+                                       * Dereferencing operator.
+                                       */
         const Accessor* operator-> () const;
 
-                                         /**
-                                          * Comparison. True, if
-                                          * both iterators point to
-                                          * the same matrix
-                                          * position.
-                                          */
+                                       /**
+                                       * Comparison. True, if both
+                                       * iterators point to the same
+                                       * matrix position.
+                                       */
         bool operator == (const const_iterator&) const;
-                                         /**
-                                          * Inverse of <tt>==</tt>.
-                                          */
+
+                                       /**
+                                       * 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.
-                                          */
+                                       /**
+                                       * 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
-                                          */
+                                      /**
+                                        * Exception
+                                       */
         DeclException2 (ExcInvalidIndexWithinRow,
                         int, int,
                         << "Attempt to access element " << arg2
@@ -250,10 +250,10 @@ namespace TrilinosWrappers
                         << " which doesn't have that many elements.");
         
       private:
-                                         /**
-                                          * Store an object of the
-                                          * accessor class.
-                                          */
+                                       /**
+                                        * Store an object of the
+                                        * accessor class.
+                                       */
         Accessor accessor;
     };
     
@@ -261,33 +261,34 @@ namespace TrilinosWrappers
   
   
 /**
- * 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 be written and hence need to be
- * forwarded to the owner process. 
- * 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 has been installed with the same compiler as is used
- * for generating deal.II.
+ * 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 be written and hence need to be forwarded to the owner
+ * process.  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 has been installed with the same compiler
+ * as is used for generating deal.II.
  *
  * 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.
+ * 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.
+ * 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
@@ -297,23 +298,25 @@ namespace TrilinosWrappers
   {
     public:
                                        /**
-                                        * A structure that describes some of
-                                        * the traits of this class in terms of
-                                        * its run-time behavior. Some other
-                                        * classes (such as the block matrix
-                                        * classes) that take one or other of
-                                        * the matrix classes as its template
-                                        * parameters can tune their behavior
-                                        * based on the variables in this
-                                        * class.
+                                        * A structure that describes
+                                        * some of the traits of this
+                                        * class in terms of its
+                                        * run-time behavior. Some
+                                        * other classes (such as the
+                                        * block matrix classes) that
+                                        * take one or other of the
+                                        * matrix classes as its
+                                        * template parameters can tune
+                                        * their behavior based on the
+                                        * variables in this class.
                                         */
       struct Traits
       {
-                                           /**
-                                            * It is safe to elide additions of
-                                            * zeros to individual elements of
-                                            * this matrix.
-                                            */
+                                       /**
+                                        * It is safe to elide additions of
+                                        * zeros to individual elements of
+                                       * this matrix.
+                                       */
           static const bool zero_addition_can_be_elided = true;
       };
 
@@ -324,8 +327,9 @@ namespace TrilinosWrappers
       typedef MatrixIterators::const_iterator const_iterator;
 
                                        /**
-                                        * Declare a typedef in analogy to all
-                                        * the other container classes.
+                                        * Declare a typedef in analogy
+                                        * to all the other container
+                                        * classes.
                                         */
       typedef TrilinosScalar value_type;
       
@@ -335,134 +339,152 @@ namespace TrilinosWrappers
       SparseMatrix ();
 
                                        /**
-                                        * Constructor using an Epetra_Map
-                                       * and a maximum number of nonzero
-                                       * matrix entries. Note that this
-                                       * number does not need to be exact,
-                                       * and it is even allowed that 
-                                       * the actual matrix structure 
-                                       * has more nonzero entries than
-                                       * specified in the constructor.
-                                       * However it is still advantageous to
-                                       * provide good estimates here since
-                                       * this will considerably increase
-                                       * the performance of the matrix.
+                                        * Constructor using an
+                                       * Epetra_Map and a maximum
+                                       * number of nonzero matrix
+                                       * entries. Note that this
+                                       * number does not need to be
+                                       * exact, and it is even
+                                       * allowed that the actual
+                                       * matrix structure has more
+                                       * nonzero entries than
+                                       * specified in the
+                                       * constructor.  However it is
+                                       * still advantageous to
+                                       * provide good estimates here
+                                       * since this will considerably
+                                       * increase the performance of
+                                       * the matrix.
                                         */
       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 
+                                        * 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);
 
                                        /**
-                                        * This function is similar to the
-                                       * one above, but it now takes two
-                                       * different Epetra maps for rows
-                                       * and columns. This interface is
-                                       * meant to be used for generating
-                                       * rectangular matrices, where one
-                                       * map specifies the parallel 
-                                       * distribution of rows and the other
-                                       * the one of the columns. This is
-                                       * in contrast to the first 
-                                       * constructor, where the same map
-                                       * is used for both the number of 
-                                       * rows and the number of columns.
-                                       * The number
-                                       * of columns per row is specified
-                                       * by the maximum number of entries.
+                                        * This function is similar to
+                                       * the one above, but it now
+                                       * takes two different Epetra
+                                       * maps for rows and
+                                       * columns. This interface is
+                                       * meant to be used for
+                                       * generating rectangular
+                                       * matrices, where one map
+                                       * specifies the parallel
+                                       * distribution of rows and the
+                                       * other the one of the
+                                       * columns. This is in contrast
+                                       * to the first constructor,
+                                       * where the same map is used
+                                       * for both the number of rows
+                                       * and the number of columns.
+                                       * The number of columns per
+                                       * row is specified by the
+                                       * maximum number of entries.
                                         */
       SparseMatrix (const Epetra_Map   &InputRowMap,
                    const Epetra_Map   &InputColMap,
                    const unsigned int  n_max_entries_per_row);
 
                                        /**
-                                        * This function is similar to the
-                                       * one above, but it now takes two
-                                       * different Epetra maps for rows
-                                       * and columns. This interface is
-                                       * meant to be used for generating
-                                       * rectangular matrices, where one
-                                       * map specifies the parallel 
-                                       * distribution of rows and the other
-                                       * the one of the columns. The 
-                                       * vector n_entries_per_row specifies
-                                       * the number of entries in each 
-                                       * row of the newly generated matrix.
+                                        * This function is similar to
+                                       * the one above, but it now
+                                       * takes two different Epetra
+                                       * maps for rows and
+                                       * columns. This interface is
+                                       * meant to be used for
+                                       * generating rectangular
+                                       * matrices, where one map
+                                       * specifies the parallel
+                                       * distribution of rows and the
+                                       * other the one of the
+                                       * columns. The vector
+                                       * n_entries_per_row specifies
+                                       * the number of entries in
+                                       * each row of the newly
+                                       * generated matrix.
                                         */
       SparseMatrix (const Epetra_Map                &InputRowMap,
                    const Epetra_Map                &InputColMap,
                    const std::vector<unsigned int> &n_entries_per_row);
 
                                        /**
-                                        * Destructor. Made virtual so that one
-                                        * can use pointers to this class.
+                                        * Destructor. Made virtual so
+                                        * that one can use pointers to
+                                        * this class.
                                         */
       virtual ~SparseMatrix ();
 
                                       /**
-                                       * Copy the given matrix to this
-                                       * one.
+                                       * Copy the given matrix to
+                                       * this one.
                                        *
-                                       * The function returns a reference to
-                                       * <tt>*this</tt>.
+                                       * The function returns a
+                                       * reference to <tt>*this</tt>.
                                        */
       SparseMatrix &
       copy_from (const SparseMatrix &source);
 
                                        /**
-                                        * This function initializes the
-                                       * Trilinos matrix with a deal.II
-                                       * sparsity pattern, i.e. it 
-                                       * makes the Trilinos Epetra
-                                       * matrix know the position of
-                                       * nonzero entries according to
-                                       * the sparsity pattern. Note that,
-                                       * when using this function, the
-                                       * matrix must already be initialized
-                                       * with a suitable Epetra_Map that
-                                       * describes the distribution of 
-                                       * the matrix among the MPI 
+                                        * This function initializes
+                                       * the Trilinos matrix with a
+                                       * deal.II sparsity pattern,
+                                       * i.e. it makes the Trilinos
+                                       * Epetra matrix know the
+                                       * position of nonzero entries
+                                       * according to the sparsity
+                                       * pattern. Note that, when
+                                       * using this function, the
+                                       * matrix must already be
+                                       * initialized with a suitable
+                                       * Epetra_Map that describes
+                                       * the distribution of the
+                                       * matrix among the MPI
                                        * processes. Otherwise, an
                                        * error will be thrown.
                                         */
       void reinit (const SparsityPattern &sparsity_pattern);
 
                                       /**
-                                        * This function is initializes the
-                                       * Trilinos Epetra matrix according
-                                       * to the specified sparsity_pattern,
-                                       * and also reassigns the matrix 
-                                       * rows to different processes 
+                                        * This function is initializes
+                                       * the Trilinos Epetra matrix
+                                       * according to the specified
+                                       * sparsity_pattern, and also
+                                       * reassigns the matrix rows to
+                                       * different processes
                                        * according to a user-supplied
-                                       * Epetra map. This might be useful
-                                       * when the matrix structure changes,
-                                       * e.g. when the grid is refined.
+                                       * Epetra map. This might be
+                                       * useful when the matrix
+                                       * structure changes, e.g. when
+                                       * the grid is refined.
                                         */
       void reinit (const Epetra_Map       &input_map,
                   const SparsityPattern  &sparsity_pattern);
 
                                       /**
-                                        * This function is similar to the
-                                       * other initialization function above,
-                                       * but now also reassigns the matrix 
-                                       * rows and columns according to 
-                                       * two user-supplied Epetra maps.
-                                       * To be used for rectangular 
+                                        * This function is similar to
+                                       * the other initialization
+                                       * function above, but now also
+                                       * reassigns the matrix rows
+                                       * and columns according to two
+                                       * user-supplied Epetra maps.
+                                       * To be used for rectangular
                                        * matrices.
                                         */
       void reinit (const Epetra_Map       &input_row_map,
@@ -470,36 +492,38 @@ namespace TrilinosWrappers
                   const SparsityPattern  &sparsity_pattern);
 
                                       /** 
-                                       * This function copies the 
-                                       * content in <tt>sparse_matrix</tt>
-                                       * to the current matrix.
+                                       * This function copies the
+                                       * content in
+                                       * <tt>sparse_matrix</tt> to
+                                       * the current matrix.
                                        */
       void reinit (const SparseMatrix &sparse_matrix);
 
                                       /**
-                                        * This function initializes the
-                                       * Trilinos matrix using the deal.II
-                                       * sparse matrix and the entries stored
-                                       * therein. It uses a threshold 
-                                       * to copy only elements with 
-                                       * modulus larger than the 
-                                       * threshold (so zeros in the 
-                                       * deal.II matrix can be filtered
-                                       * away).
+                                        * This function initializes
+                                       * the Trilinos matrix using
+                                       * the deal.II sparse matrix
+                                       * and the entries stored
+                                       * therein. It uses a threshold
+                                       * to copy only elements with
+                                       * modulus larger than the
+                                       * threshold (so zeros in the
+                                       * deal.II matrix can be
+                                       * filtered away).
                                         */
       void reinit (const Epetra_Map                     &input_map,
                   const ::dealii::SparseMatrix<double> &dealii_sparse_matrix,
                   const double                          drop_tolerance=1e-13);
 
                                       /**
-                                        * This function is similar to the
-                                       * other initialization function with 
-                                       * deal.II sparse matrix input above,
-                                       * but now takes Epetra maps for both
-                                       * the rows and the columns of the
-                                       * matrix.
-                                       * To be used for rectangular 
-                                       * matrices.
+                                        * This function is similar to
+                                       * the other initialization
+                                       * function with deal.II sparse
+                                       * matrix input above, but now
+                                       * takes Epetra maps for both
+                                       * the rows and the columns of
+                                       * the matrix. To be used for
+                                       * rectangular matrices.
                                         */
       void reinit (const Epetra_Map                      &input_row_map,
                   const Epetra_Map                      &input_col_map,
@@ -507,41 +531,49 @@ namespace TrilinosWrappers
                   const double                           drop_tolerance=1e-13);
 
                                        /**
-                                        * Release all memory and return
-                                        * to a state just like after
-                                        * having called the default
-                                        * constructor.
+                                        * Release all memory and
+                                        * return to a state just like
+                                        * after having called the
+                                        * default constructor.
                                         */
       void clear ();
 
                                        /**
-                                        * 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
+                                        * 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 ();
 
                                        /**
-                                        * 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
+                                        * 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 &
@@ -551,14 +583,14 @@ namespace TrilinosWrappers
                                         * 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
+                                       * 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.
                                        */
@@ -567,17 +599,17 @@ namespace TrilinosWrappers
                 const TrilinosScalar value);
 
                                        /**
-                                        * Add @p value to the
-                                        * element (<i>i,j</i>).
+                                        * 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
+                                       * 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.
                                         */
@@ -597,79 +629,96 @@ namespace TrilinosWrappers
                                         * pattern, though (but retains
                                         * the allocated memory in case
                                         * new entries are again added
-                                        * later). Note that this
-                                       * is a global operation, so this
-                                       * needs to be done on all 
-                                       * MPI processes.
+                                        * later). Note that this is a
+                                        * global operation, so this
+                                        * needs to be done on all MPI
+                                        * processes.
                                         *
                                         * 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 o peration, removing
-                                        * constraints on parallel matrices is
-                                        * a rather complicated procedure.
+                                        * 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 o peration, 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.
+                                        * 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.
+                                        * 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.
+                                        * 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);
 
                                        /**
-                                        * 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.
+                                        * 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.
+                                        * 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.
+                                        * 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;
@@ -680,32 +729,30 @@ namespace TrilinosWrappers
                                         * 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.
+                                        * Return the number of rows in
+                                        * this matrix.
                                         */
       unsigned int m () const;
 
                                        /**
-                                        * Return the number of columns in this
-                                        * matrix.
+                                        * 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.
+                                        * 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,
@@ -718,15 +765,15 @@ namespace TrilinosWrappers
                                        * 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
+                                       * 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>
@@ -741,24 +788,25 @@ namespace TrilinosWrappers
 
                                        /**
                                         * Return the number of nonzero
-                                        * elements of this
-                                        * matrix. 
+                                        * elements of this matrix.
                                         */
       unsigned int n_nonzero_elements () const;
 
                                        /**
-                                        * Number of entries in a specific row.
+                                        * 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.
+                                        * 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)
@@ -766,76 +814,87 @@ namespace TrilinosWrappers
       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$.
+                                        * 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
+                                        * 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.
+                                        * Multiply the entire matrix
+                                        * by a fixed factor.
                                         */
       SparseMatrix & operator *= (const TrilinosScalar factor);
     
                                        /**
-                                        * Divide the entire matrix by a
-                                        * fixed factor.
+                                        * Divide the entire matrix by
+                                        * fixed factor.
                                         */
       SparseMatrix & operator /= (const TrilinosScalar factor);
 
                                        /**
-                                        * Matrix-vector multiplication:
-                                        * let <i>dst = M*src</i> with
-                                        * <i>M</i> being this matrix.
+                                        * 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.
+                                       * 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.
+                                        * 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.
+                                       * 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;
@@ -850,13 +909,16 @@ namespace TrilinosWrappers
                                         * 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.
+                                       * 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;
@@ -865,59 +927,70 @@ namespace TrilinosWrappers
                                         * 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.
+                                        * <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.
+                                       * 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,
+                                        * 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.
+                                        * 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.
+                                        * 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.
+                                        * 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.
+                                       * 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;
 
@@ -925,21 +998,27 @@ namespace TrilinosWrappers
                                         * 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.
+                                        * 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.
+                                       * 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;
@@ -947,35 +1026,38 @@ namespace TrilinosWrappers
                                        /**
                                         * 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
+                                        * 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.
+                                        * 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.
+                                       * 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;
 
                                       /**
-                                       * Add <tt>matrix</tt> scaled by
-                                       * <tt>factor</tt> to this matrix,
-                                       * i.e. the matrix <tt>factor*matrix</tt>
-                                       * is added to <tt>this</tt>.
+                                       * Add <tt>matrix</tt> scaled
+                                       * by <tt>factor</tt> to this
+                                       * matrix, i.e. the matrix
+                                       * <tt>factor*matrix</tt> is
+                                       * added to <tt>this</tt>.
                                        */
       void add (const TrilinosScalar  factor,
                const SparseMatrix   &matrix);
@@ -995,64 +1077,75 @@ namespace TrilinosWrappers
                                         * 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.
+                                        * 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.
+                                        * 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.
+                                        * 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.
-                                        */
+                                       /**  
+                                       * Make an in-place transpose
+                                       * of a matrix.
+                                       */
       void transpose ();
 
-                                        /**  
-                                        * Test whether a matrix is symmetric.
-                                        * Default tolerance is zero. 
-                                        * TODO: Not implemented.
-                                        */
+                                       /**  
+                                       * 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.
-                                        */
+                                       /** 
+                                       * 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.
-                                        */
+                                       /** 
+                                       * Abstract Trilinos object
+                                       * that helps view in ASCII
+                                       * other Trilinos
+                                       * objects. Currently this
+                                       * function is not
+                                       * implemented.  TODO: Not
+                                       * implemented.
+                                       */
       void write_ascii ();
       
                                       /**
-                                       * Print the matrix to the given
-                                       * stream, using the format
-                                       * <tt>(line,col) value</tt>,
-                                       * i.e. one nonzero entry of the
-                                       * matrix per line.
+                                       * Print the matrix to the
+                                       * given stream, using the
+                                       * format <tt>(line,col)
+                                       * value</tt>, i.e. one nonzero
+                                       * entry of the matrix per
+                                       * line.
                                        */
       void print (std::ostream &out) const;
     
@@ -1098,55 +1191,71 @@ namespace TrilinosWrappers
 
     private:
                                        /**
-                                       * Epetra Trilinos mapping of the 
-                                       * matrix rows that
-                                       * assigns parts of the matrix to
-                                       * the individual processes.
-                                       * TODO: is it possible to only use
-                                       * a pointer?
-                                       */
-      Epetra_Map row_map;
+                                       * Dummy pointer that is used
+                                       * internally in the
+                                       * constructor for empty
+                                       * matrices.
+                                       */
+      const Epetra_Map *dummy_map;
 
                                        /**
-                                        * Pointer to the user-supplied 
-                                       * Epetra Trilinos mapping of the 
-                                       * matrix columns that
-                                       * assigns parts of the matrix to
-                                       * the individual processes.
+                                       * Pointer to Epetra Trilinos
+                                       * mapping of the matrix rows
+                                       * that assigns parts of the
+                                       * matrix to the individual
+                                       * processes. This map is
+                                       * provided either via the
+                                       * constructor or in a reinit
+                                       * function.
                                        */
-      Epetra_Map col_map;
+      Epetra_Map *row_map;
 
                                        /**
-                                        * Trilinos doesn't allow to mix additions
-                                        * to matrix entries and overwriting
-                                        * them (to make synchronisation of
+                                        * Pointer to the user-supplied
+                                       * Epetra Trilinos mapping of
+                                       * the matrix columns that
+                                       * assigns parts of the matrix
+                                       * to the individual processes.
+                                       */
+      Epetra_Map *col_map;
+  
+                                       /**
+                                        * Trilinos doesn't allow to
+                                        * mix additions to matrix
+                                        * entries and overwriting them
+                                        * (to make synchronisation of
                                         * parallel computations
-                                        * simpler). The way we do it is to,
-                                        * for each access operation, store
-                                        * whether it is an insertion or an
-                                        * addition. If the previous one was of
-                                        * different type, then we first have
-                                        * to flush the Trilinos buffers;
-                                        * otherwise, we can simply go on.
-                                       * Luckily, Trilinos has an object
-                                       * for this which does already all
-                                       * the parallel communications in
-                                       * such a case, so we simply use their
-                                       * model, which stores whether the
-                                       * last operation was an addition
-                                       * or an insertion.
+                                        * 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;
 
     public:
                                        /**
                                         * 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.
+                                        * 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;
 
@@ -1163,7 +1272,7 @@ namespace TrilinosWrappers
 
     inline
     const_iterator::Accessor::
-    Accessor (const SparseMatrix   *matrix,
+    Accessor (const SparseMatrix *matrix,
               const unsigned int  row,
               const unsigned int  index)
                     :
index 4e3fe402508b2f41e62a64dcc2366a998f031761..7871e248db5302df0e6fc06f50312780d271d814 100755 (executable)
@@ -51,113 +51,124 @@ namespace TrilinosWrappers
                                    // forward declaration
   class Vector;
 
-                                  /**
-                                   * @cond internal
-                                   */
+                                      /**
+                                       * @cond internal
+                                       */
 
 /**
- * A namespace for internal implementation details of the TrilinosWrapper
- * members.
- * @ingroup TrilinosWrappers
+ * 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
-                                      */
+                                       /**
+                                       * 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.
-                                          */
+                                       /**
+                                       * 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.
-                                          */
+                                       /**
+                                       * 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>.
-                                          */
+                                       /**
+                                       * 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->
-                                          */
+                                       /**
+                                       * 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->
-                                          */
+                                       /**
+                                       * 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>.
-                                          */
+                                      /**
+                                       * 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>.
-                                          */
+                                       /**
+                                       * 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.
-                                          */
+                                       /**
+                                       * Convert the reference to an
+                                       * actual value, i.e. return
+                                       * the value of the referenced
+                                       * element of the vector.
+                                       */
         operator TrilinosScalar () const;
 
-                                         /**
-                                          * Exception
-                                          */
+                                       /**
+                                       * Exception
+                                       */
         DeclException1 (ExcTrilinosError,
                         int,
                         << "An error with error number " << arg1
                         << " occured while calling a Trilinos function");
 
-                                         /**
-                                          * Exception
-                                          */
+                                       /**
+                                       * Exception
+                                       */
         DeclException3 (ExcAccessToNonLocalElement,
                         int, int, int,
                         << "You tried to access element " << arg1
@@ -166,111 +177,116 @@ namespace TrilinosWrappers
                         << " are stored locally and can be accessed.");
 
       private:
-                                         /**
-                                          * Point to the vector we are
-                                          * referencing.
-                                          */
+                                       /**
+                                       * Point to the vector we are
+                                       * referencing.
+                                       */
         Vector   &vector;
 
-                                         /**
-                                          * Index of the referenced element
-                                          * of the 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.
-                                          */
+                                       /**
+                                       * Make the vector class a
+                                       * friend, so that it can
+                                       * create objects of the
+                                       * present type.
+                                       */
         friend class ::dealii::TrilinosWrappers::Vector;
     };
   }
-                                   /**
-                                    * @endcond
-                                    */
+                                       /**
+                                       * @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.
+ * 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.
+ * 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.
+ * 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.
+ * 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.
+ * 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.
+ * 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 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;
+ * 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
@@ -282,23 +298,24 @@ namespace TrilinosWrappers
  *   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.
+ * 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.
+ * 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
  * @ingroup Vectors
@@ -309,11 +326,13 @@ namespace TrilinosWrappers
   {
     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.
+                                        * 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;
@@ -322,47 +341,35 @@ namespace TrilinosWrappers
       typedef const internal::VectorReference const_reference;
 
                                        /**
-                                        * Default constructor. It doesn't do
-                                        * anything, derived classes will have
-                                        * to initialize the data.
+                                        * Default constructor. It
+                                        * doesn't do anything, derived
+                                        * classes will have to
+                                        * initialize the data.
                                         */
       Vector ();
-                                       /**
-                                        * One of the constructors that
-                                       * actually builds 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
-                                       * constructor with argument
-                                       * Epetra_Map.
-                                        */
-      Vector (unsigned int GlobalSize, Epetra_Comm &Comm);
 
                                        /**
                                        * This constructor takes an
-                                       * Epetra_Map that already knows how
-                                       * to distribute the individual 
-                                       * components among the MPI processors,
-                                       * including the size of the vector.
+                                       * Epetra_Map that already
+                                       * knows how to distribute the
+                                       * individual components among
+                                       * the MPI processors,
+                                       * including the size of the
+                                       * vector.
                                         */
       Vector (const Epetra_Map &InputMap);
 
                                        /**
-                                        * Copy constructor. Sets the dimension
-                                        * to that of the given vector and uses
-                                       * the map of that vector, but
-                                        * does not copy any element. Instead,
-                                       * the memory will remain untouched
-                                       * in case <tt>fast</tt> is false and
-                                       * initialized with zero otherwise.
+                                        * Copy constructor. Sets the
+                                        * dimension to that of the
+                                        * given vector and uses the
+                                        * map of that vector, but does
+                                        * not copy any
+                                        * element. Instead, the memory
+                                        * will remain untouched in
+                                        * case <tt>fast</tt> is false
+                                        * and initialized with zero
+                                        * otherwise.
                                         */
       Vector (const Vector &v, 
              const bool    fast = false);
@@ -673,55 +680,61 @@ namespace TrilinosWrappers
       Vector & operator += (const Vector &V);
 
                                        /**
-                                        * Subtract the given vector from the
-                                        * present one.
+                                        * 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.
+                                        * 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>.
+                                        * 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>.
+                                        * 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>.
+                                        * 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>.
+                                        * 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>.
+                                        * 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.
+                                        * Scaling and multiple
+                                        * addition.
                                         */
       void sadd (const TrilinosScalar s,
                  const TrilinosScalar a,
@@ -730,8 +743,10 @@ namespace TrilinosWrappers
                  const Vector        &W);
 
                                        /**
-                                        * Scaling and multiple addition.
-                                        * <tt>*this = s*(*this) + a*V + b*W + c*X</tt>.
+                                        * Scaling and multiple
+                                        * addition.  <tt>*this =
+                                        * s*(*this) + a*V + b*W +
+                                        * c*X</tt>.
                                         */
       void sadd (const TrilinosScalar s,
                  const TrilinosScalar a,
@@ -744,64 +759,72 @@ namespace TrilinosWrappers
                                        /**
                                         * 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.
+                                        * 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>.
+                                        * Assignment <tt>*this =
+                                        * a*V</tt>.
                                         */
       void equ (const TrilinosScalar a, const Vector &V);
 
                                        /**
-                                        * Assignment <tt>*this = a*V + b*W</tt>.
+                                        * 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.
+                                        * 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.
+                                        * 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);
 
                                     /**
-                                     *  Output of vector in user-defined
-                                     *  format in analogy to the 
+                                     *  Output of vector in
+                                     *  user-defined format in
+                                     *  analogy to the
                                      *  dealii::Vector<number> class.
                                      */
     void print (const char* format = 0) const;
 
                                        /**
-                                        * Print to a
-                                        * stream. @p precision denotes
-                                        * the desired precision with
-                                        * which values shall be printed,
-                                        * @p scientific whether
+                                        * 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
+                                        * 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.
+                                        * are printed on a separate
+                                        * line each.
                                         */
       void print (std::ostream       &out,
                   const unsigned int  precision  = 3,
@@ -813,18 +836,18 @@ namespace TrilinosWrappers
                                         * 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.
+                                        * 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
@@ -864,49 +887,65 @@ namespace TrilinosWrappers
 
     private:
                                        /**
-                                        * 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.
+                                       * A dummy pointer to an
+                                       * Epetra_Map to be used for
+                                       * the default constructor.
+                                       */
+      Epetra_Map *dummy_map;
+
+                                       /**
+                                        * A pointer to 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.
                                         */
-      Epetra_Map map;
+      Epetra_Map *map;
 
                                       /**
-                                        * Trilinos doesn't allow to mix additions
-                                        * to matrix entries and overwriting
-                                        * them (to make synchronisation of
+                                        * 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.
+                                        * 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;
 
     public:
                                        /**
-                                        * An Epetra distibuted vector type.
-                                        * Requires an existing Epetra_Map for
-                                        * storing data.
-                                       * TODO: Should become private at some
-                                       * point.
+                                        * An Epetra distibuted vector
+                                        * type.  Requires an existing
+                                        * Epetra_Map for storing data.
+                                        * TODO: Should become private
+                                        * at some point.
                                         */
       std::auto_ptr<Epetra_FEVector> vector;
 
 
                                        /**
-                                        * Make the reference class a friend.
+                                        * Make the reference class a
+                                        * friend.
                                         */
       friend class internal::VectorReference;
 
@@ -917,9 +956,9 @@ namespace TrilinosWrappers
 // ------------------- inline and template functions --------------
 
 /**
- * Global function swap which overloads the default implementation
- * of the C standard library which uses a temporary object. The
- * function simply exchanges the data of the two vectors.
+ * Global function 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
@@ -947,10 +986,10 @@ namespace TrilinosWrappers
     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
+                                        // 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;
@@ -1029,10 +1068,10 @@ namespace TrilinosWrappers
   {
                                    // if the vectors have different sizes,
                                    // then first resize the present one
-    if (!map.SameAs(v.map))
+    if (! map->SameAs(*v.map))
       {
        map = v.map;
-         vector = std::auto_ptr<Epetra_FEVector> 
+       vector = std::auto_ptr<Epetra_FEVector> 
                        (new Epetra_FEVector(*v.vector));
       }
     else
index d9c44149a1e63b64e4886cca41cd1d72c3ae15ce..c265a76b1290792b69c33dc40ab40ea2d23d3c57 100755 (executable)
@@ -53,7 +53,7 @@ namespace TrilinosWrappers
     const unsigned int n_rows = matrix.m();
     const unsigned int null_space_dimension = null_space.size();
 
-                                    // Build the AMG preconditioner.
+                                       // Build the AMG preconditioner.
     Teuchos::ParameterList parameter_list;
   
     if (elliptic)
@@ -87,9 +87,10 @@ namespace TrilinosWrappers
                ExcDimensionMismatch(n_rows,
                                     null_space[0].size()));
        
-                                // Reshape null space as a contiguous
-                                // vector of doubles so that Trilinos
-                                // can read from it.
+                                       // Reshape null space as a
+                                       // contiguous vector of
+                                       // doubles so that Trilinos
+                                       // can read from it.
        null_space_modes.resize (n_rows * null_space_dimension, 0.);
        for (unsigned int d=0; d<null_space_dimension; ++d)
          for (unsigned int row=0; row<n_rows; ++row)
@@ -121,8 +122,9 @@ namespace TrilinosWrappers
   {
     const unsigned int n_rows = deal_ii_sparse_matrix.m();
   
-                                    // Init Epetra Matrix, avoid 
-                                    // storing the nonzero elements.
+                                       // Init Epetra Matrix, avoid
+                                       // storing the nonzero
+                                       // elements.
 
     Map.reset (new Epetra_Map(n_rows, 0, communicator));
 
@@ -155,23 +157,27 @@ namespace TrilinosWrappers
   }
 
 
-                                  // For the implementation of the
-                                  // <code>vmult</code> function we
-                                  // note that invoking a call of 
-                                  // the Trilinos preconditioner 
-                                  // requires us to use Epetra vectors
-                                  // as well. It is faster
-                                  // to provide a view, i.e., feed 
-                                  // Trilinos with a pointer to the
-                                  // data, so we avoid copying the
-                                  // content of the vectors during
-                                  // the iteration. In the declaration
-                                  // of the right hand side, we need
-                                  // to cast the source vector (that
-                                  // is <code>const</code> in all deal.II 
-                                  // calls) to non-constant value, as
-                                  // this is the way Trilinos wants to
-                                  // have them.
+                                       // For the implementation of
+                                       // the <code>vmult</code>
+                                       // function we note that
+                                       // invoking a call of the
+                                       // Trilinos preconditioner
+                                       // requires us to use Epetra
+                                       // vectors as well. It is
+                                       // faster to provide a view,
+                                       // i.e., feed Trilinos with a
+                                       // pointer to the data, so we
+                                       // avoid copying the content
+                                       // of the vectors during the
+                                       // iteration. In the
+                                       // declaration of the right
+                                       // hand side, we need to cast
+                                       // the source vector (that is
+                                       // <code>const</code> in all
+                                       // deal.II calls) to
+                                       // non-constant value, as this
+                                       // is the way Trilinos wants
+                                       // to have them.
   void PreconditionAMG::vmult (dealii::Vector<double>       &dst,
                               const dealii::Vector<double> &src) const
   {
index 9ed483cdfb4e5954331fab1113df40b04c95f385..d224a273641c1b5a6d11c248fe549ce68bf81294 100755 (executable)
@@ -60,7 +60,7 @@ namespace TrilinosWrappers
                                  // would it point to?)
       Assert (ncols != 0, ExcInternalError());
       colnum_cache.reset (new std::vector<unsigned int> (colnums,
-                                                       colnums+ncols));
+                                                        colnums+ncols));
       value_cache.reset (new std::vector<TrilinosScalar> (values, values+ncols));
     }
   }
@@ -80,35 +80,46 @@ namespace TrilinosWrappers
   SparseMatrix::SparseMatrix ()
                  :
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-                 row_map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+                  dummy_map (new Epetra_Map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD))),
 #else
-                 row_map (0, 0, Epetra_SerialComm()),
+                  dummy_map (new Epetra_Map (0, 0, Epetra_SerialComm())),
 #endif
+                 row_map (const_cast<Epetra_Map*>(dummy_map)),
                  col_map (row_map),
                  last_action (Insert),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                               (new Epetra_FECrsMatrix(Copy, row_map, 0)))
+                               (new Epetra_FECrsMatrix(Copy, *row_map, 0)))
   {}
 
   SparseMatrix::SparseMatrix (const Epetra_Map  &InputMap,
                              const unsigned int n_max_entries_per_row)
                  :
-                 row_map (InputMap),
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+                  dummy_map (new Epetra_Map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD))),
+#else
+                  dummy_map (new Epetra_Map (0, 0, Epetra_SerialComm())),
+#endif
+                  row_map (const_cast<Epetra_Map*>(&InputMap)),
                  col_map (row_map),
                  last_action (Insert),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                               (new Epetra_FECrsMatrix(Copy, row_map, 
+                               (new Epetra_FECrsMatrix(Copy, *row_map, 
                                        int(n_max_entries_per_row), false)))
   {}
 
-  SparseMatrix::SparseMatrix (const Epetra_Map  &InputMap,
+  SparseMatrix::SparseMatrix (const Epetra_Map                &InputMap,
                              const std::vector<unsigned int> &n_entries_per_row)
                  :
-                 row_map (InputMap),
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+                  dummy_map (new Epetra_Map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD))),
+#else
+                  dummy_map (new Epetra_Map (0, 0, Epetra_SerialComm())),
+#endif
+                  row_map (const_cast<Epetra_Map*>(&InputMap)),
                  col_map (row_map),
                  last_action (Insert),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                   (new Epetra_FECrsMatrix(Copy, row_map, 
+                   (new Epetra_FECrsMatrix(Copy, *row_map, 
                      (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
                      true)))
   {}
@@ -117,23 +128,33 @@ namespace TrilinosWrappers
                              const Epetra_Map  &InputColMap,
                              const unsigned int n_max_entries_per_row)
                  :
-                 row_map (InputRowMap),
-                 col_map (InputColMap),
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+                  dummy_map (new Epetra_Map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD))),
+#else
+                  dummy_map (new Epetra_Map (0, 0, Epetra_SerialComm())),
+#endif
+                  row_map (const_cast<Epetra_Map*>(&InputRowMap)),
+                  col_map (const_cast<Epetra_Map*>(&InputColMap)),
                  last_action (Insert),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                               (new Epetra_FECrsMatrix(Copy, row_map, col_map, 
+                               (new Epetra_FECrsMatrix(Copy, *row_map, *col_map, 
                                        int(n_max_entries_per_row), false)))
   {}
 
-  SparseMatrix::SparseMatrix (const Epetra_Map  &InputRowMap,
-                             const Epetra_Map  &InputColMap,
+  SparseMatrix::SparseMatrix (const Epetra_Map                &InputRowMap,
+                             const Epetra_Map                &InputColMap,
                              const std::vector<unsigned int> &n_entries_per_row)
                  :
-                 row_map (InputRowMap),
-                 col_map (InputColMap),
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+                  dummy_map (new Epetra_Map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD))),
+#else
+                  dummy_map (new Epetra_Map (0, 0, Epetra_SerialComm())),
+#endif
+                  row_map (const_cast<Epetra_Map*>(&InputRowMap)),
+                  col_map (const_cast<Epetra_Map*>(&InputColMap)),
                  last_action (Insert),
                  matrix (std::auto_ptr<Epetra_FECrsMatrix>
-                   (new Epetra_FECrsMatrix(Copy, row_map, col_map, 
+                   (new Epetra_FECrsMatrix(Copy, *row_map, *col_map, 
                      (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
                      true)))
   {}
@@ -142,6 +163,7 @@ namespace TrilinosWrappers
 
   SparseMatrix::~SparseMatrix ()
   {
+    delete dummy_map;
   }
 
 
@@ -167,7 +189,7 @@ namespace TrilinosWrappers
                                  // sparsity pattern on that
                                  // processor, and only broadcast the
                                  // pattern afterwards.
-    if (row_map.Comm().MyPID() == 0)
+    if (row_map->Comm().MyPID() == 0)
       {
        Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
                ExcDimensionMismatch (matrix->NumGlobalRows(),
@@ -228,11 +250,12 @@ namespace TrilinosWrappers
                                  // column map. Maybe find something
                                  // more out about this...
     //reinit (input_map, input_map, sparsity_pattern);
-    matrix.reset();
+    
+    matrix.reset();   
 
     unsigned int n_rows = sparsity_pattern.n_rows();
 
-    if (row_map.Comm().MyPID() == 0)
+    if (row_map->Comm().MyPID() == 0)
       {
        Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
                ExcDimensionMismatch (input_map.NumGlobalElements(),
@@ -242,8 +265,8 @@ namespace TrilinosWrappers
                                      sparsity_pattern.n_cols()));
       }
 
-    row_map = input_map;
-    col_map = input_map;
+    row_map = const_cast<Epetra_Map*>(&input_map);
+    col_map = row_map;
 
     std::vector<int> n_entries_per_row(n_rows);
 
@@ -251,7 +274,7 @@ namespace TrilinosWrappers
       n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
 
     matrix = std::auto_ptr<Epetra_FECrsMatrix>
-             (new Epetra_FECrsMatrix(Copy, row_map, &n_entries_per_row[0], 
+             (new Epetra_FECrsMatrix(Copy, *row_map, &n_entries_per_row[0], 
                                      false));
 
     reinit (sparsity_pattern);
@@ -268,7 +291,7 @@ namespace TrilinosWrappers
 
     unsigned int n_rows = sparsity_pattern.n_rows();
 
-    if (row_map.Comm().MyPID() == 0)
+    if (input_row_map.Comm().MyPID() == 0)
       {
        Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
                ExcDimensionMismatch (input_row_map.NumGlobalElements(),
@@ -278,8 +301,8 @@ namespace TrilinosWrappers
                                      sparsity_pattern.n_cols()));
       }
 
-    row_map = input_row_map;
-    col_map = input_col_map;
+    row_map = const_cast<Epetra_Map*>(&input_row_map);
+    col_map = const_cast<Epetra_Map*>(&input_col_map);
 
     std::vector<int> n_entries_per_row(n_rows);
 
@@ -287,7 +310,7 @@ namespace TrilinosWrappers
       n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
 
     matrix = std::auto_ptr<Epetra_FECrsMatrix>
-             (new Epetra_FECrsMatrix(Copy, row_map, col_map,
+             (new Epetra_FECrsMatrix(Copy, *row_map, *col_map,
                                      &n_entries_per_row[0], false));
 
     reinit (sparsity_pattern);
@@ -334,8 +357,8 @@ namespace TrilinosWrappers
            ExcDimensionMismatch (input_col_map.NumGlobalElements(),
                                  dealii_sparse_matrix.n()));
 
-    row_map = input_row_map;
-    col_map = input_col_map;
+    row_map = const_cast<Epetra_Map*>(&input_row_map);
+    col_map = const_cast<Epetra_Map*>(&input_col_map);
 
     std::vector<int> n_entries_per_row(n_rows);
 
@@ -344,7 +367,7 @@ namespace TrilinosWrappers
          dealii_sparse_matrix.get_sparsity_pattern().row_length(row);
 
     matrix = std::auto_ptr<Epetra_FECrsMatrix>
-             (new Epetra_FECrsMatrix(Copy, row_map, col_map,
+             (new Epetra_FECrsMatrix(Copy, *row_map, *col_map,
                                      &n_entries_per_row[0], true));
 
     std::vector<double> values;
@@ -382,14 +405,12 @@ namespace TrilinosWrappers
                                  // the pointer and generate an
                                  // empty matrix.
     matrix.reset();
-#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-    row_map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
-#else
-    row_map = Epetra_Map (0,0,Epetra_SerialComm()),
-#endif
+    
+    row_map = const_cast<Epetra_Map*>(dummy_map);
+    col_map = row_map;
 
     matrix = std::auto_ptr<Epetra_FECrsMatrix> 
-             (new Epetra_FECrsMatrix(Copy, row_map, 0));
+             (new Epetra_FECrsMatrix(Copy, *row_map, 0));
   }
 
 
@@ -399,10 +420,10 @@ namespace TrilinosWrappers
   {
                                  // flush buffers
     int ierr;
-    if (row_map.SameAs(col_map))
+    if (row_map->SameAs(*col_map))
       ierr = matrix->GlobalAssemble ();
     else
-      ierr = matrix->GlobalAssemble (col_map, row_map);
+      ierr = matrix->GlobalAssemble (*col_map, *row_map);
     
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
@@ -440,10 +461,10 @@ namespace TrilinosWrappers
     if (last_action == Add)
       {
        int ierr;
-       if (row_map.SameAs(col_map))
+       if (row_map->SameAs(*col_map))
          ierr = matrix->GlobalAssemble (false);
        else
-         ierr = matrix->GlobalAssemble(col_map, row_map, false);
+         ierr = matrix->GlobalAssemble(*col_map, *row_map, false);
        
        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
        last_action = Insert;
@@ -475,10 +496,10 @@ namespace TrilinosWrappers
     if (last_action == Insert)
       {
        int ierr;
-       if (row_map.SameAs(col_map))
+       if (row_map->SameAs(*col_map))
          ierr = matrix->GlobalAssemble (false);
        else
-         ierr = matrix->GlobalAssemble(col_map, row_map, false);
+         ierr = matrix->GlobalAssemble(*col_map, *row_map, false);
        
        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
 
@@ -623,6 +644,10 @@ namespace TrilinosWrappers
   {
     Assert (m() == n(), ExcNotQuadratic());
 
+                                 // Trilinos doesn't seem to have a
+                                 // more efficient way to access the
+                                 // diagonal than by just using the
+                                 // standard el(i,j) function.
     return el(i,i);
   }
 
@@ -684,14 +709,14 @@ namespace TrilinosWrappers
   {
     Assert (row < m(), ExcInternalError());
 
-                               // get a representation of the present
-                               // row
+                                 // 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.
+                                 // 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);
@@ -772,9 +797,9 @@ namespace TrilinosWrappers
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
     
-    Assert (col_map.SameAs(src.vector->Map()),
+    Assert (col_map->SameAs(src.vector->Map()),
            ExcMessage ("Column map of matrix does not fit with vector map!"));
-    Assert (row_map.SameAs(dst.vector->Map()),
+    Assert (row_map->SameAs(dst.vector->Map()),
            ExcMessage ("Row map of matrix does not fit with vector map!"));
 
     if (!matrix->Filled())
@@ -792,9 +817,9 @@ namespace TrilinosWrappers
   {
     Assert (&src != &dst, ExcSourceEqualsDestination());
 
-    Assert (row_map.SameAs(src.vector->Map()),
+    Assert (row_map->SameAs(src.vector->Map()),
            ExcMessage ("Row map of matrix does not fit with vector map!"));
-    Assert (col_map.SameAs(dst.vector->Map()),
+    Assert (col_map->SameAs(dst.vector->Map()),
            ExcMessage ("Column map of matrix does not fit with vector map!"));
 
     if (!matrix->Filled())
@@ -876,13 +901,15 @@ namespace TrilinosWrappers
     Assert (rhs.m() == m(), ExcDimensionMismatch (rhs.m(), m()));
     Assert (rhs.n() == n(), ExcDimensionMismatch (rhs.n(), n()));
     
-                                    // I bet that there must be a better way
-                                    // to do this but it has not been found:
-                                    // currently, we simply go through each
-                                    // row of the argument matrix, copy it,
-                                    // scale it, and add it to the current
-                                    // matrix. that's probably not the most
-                                    // efficient way to do things.
+                                 // I bet that there must be a
+                                  // better way to do this but it
+                                  // has not been found: currently,
+                                 // we simply go through each row
+                                 // of the argument matrix, copy
+                                 // it, scale it, and add it to
+                                 // the current matrix. that's
+                                 // probably not the most
+                                 // efficient way to do things.
     const std::pair<unsigned int, unsigned int>
       local_range = rhs.local_range();
 
@@ -917,11 +944,12 @@ namespace TrilinosWrappers
   
 
 
-                                 // 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.
+                                 // 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. Can we leave it like this?
   void
   SparseMatrix::transpose () 
   {
@@ -973,7 +1001,7 @@ namespace TrilinosWrappers
 
 
                                  // As of now, no particularly neat
-                                 // ouput is generated in case of 
+                                 // ouput is generated in case of
                                  // multiple processors.
   void SparseMatrix::print (std::ostream &out) const
   {
index d57c9afd14366c28c80945f7ecaf90fd1785f2c7..4d7ae7d37f1d202b30024563b333ef15cceeab92 100755 (executable)
@@ -29,63 +29,74 @@ namespace TrilinosWrappers
       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
+                                        // 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));
+      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))[0][index];
     }
   }
 
+
+
   Vector::Vector ()
                   :
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-                 map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
+                  dummy_map (new Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD))),
 #else
-                 map (0,0,Epetra_SerialComm()),
+                 dummy_map (new Epetra_Map (0,0,Epetra_SerialComm())),
 #endif
+                 map (const_cast<Epetra_Map*>(dummy_map)),
                   last_action (Insert),
                  vector(std::auto_ptr<Epetra_FEVector> 
-                        (new Epetra_FEVector(map)))
+                        (new Epetra_FEVector(*map)))
   {}
 
-  Vector::Vector (unsigned int GlobalSize, Epetra_Comm &Comm)
-                  :
-                 map (GlobalSize, 0, Comm),
-                  last_action (Insert),
-                 vector (std::auto_ptr<Epetra_FEVector> 
-                         (new Epetra_FEVector(map)))
-  {}
 
   
   Vector::Vector (const Epetra_Map &InputMap)
                   :
-                 map (InputMap),
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+                  dummy_map (new Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD))),
+#else
+                 dummy_map (new Epetra_Map (0,0,Epetra_SerialComm())),
+#endif
+                 map (const_cast<Epetra_Map*>(&InputMap)),
                   last_action (Insert),
                  vector (std::auto_ptr<Epetra_FEVector> 
-                         (new Epetra_FEVector(map)))
+                         (new Epetra_FEVector(*map)))
   {}
   
+
   
   Vector::Vector (const Vector &v,
                  const bool    fast)
                   :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+                  dummy_map (new Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD))),
+#else
+                 dummy_map (new Epetra_Map (0,0,Epetra_SerialComm())),
+#endif
                  map (v.map),
                   last_action (Insert),
                  vector(std::auto_ptr<Epetra_FEVector> 
-                        (new Epetra_FEVector(v.map,!fast)))
+                        (new Epetra_FEVector(*map,!fast)))
   {}
   
 
 
   Vector::~Vector ()
-  {}
+  {
+    vector.reset();
+    delete dummy_map;
+  }
 
 
 
@@ -93,9 +104,9 @@ namespace TrilinosWrappers
   Vector::reinit (const Epetra_Map &input_map)
   {
     vector.reset();
-    map = input_map;
+    map = const_cast<Epetra_Map*> (&input_map);
 
-    vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(input_map));
+    vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(*map));
     last_action = Insert;
   }
 
@@ -107,10 +118,9 @@ namespace TrilinosWrappers
   {
     vector.reset();
 
-    if (!map.SameAs(v.map))
-      map = v.map;
+    map = v.map;
 
-    vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(v.map,!fast));
+    vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(*map,!fast));
     last_action = Insert;
   }
 
@@ -119,17 +129,13 @@ namespace TrilinosWrappers
   void
   Vector::clear ()
   {
-                                     // When we clear the matrix,
-                                    // reset the pointer and 
-                                    // generate an empty matrix.
+                                     // When we clear the vector,
+                                    // reset the pointer and generate
+                                    // an empty vector.
     vector.reset();
-#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
-    map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)),
-#else
-    map = Epetra_Map (0,0,Epetra_SerialComm()),
-#endif
+    map = dummy_map;
 
-    vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(map));
+    vector = std::auto_ptr<Epetra_FEVector> (new Epetra_FEVector(*map));
     last_action = Insert;
   }
 
@@ -218,9 +224,9 @@ namespace TrilinosWrappers
   TrilinosScalar
   Vector::el (const unsigned int index) const
   {
-                                      // Extract local indices in
-                                      // the vector.
-    int trilinos_i = map.LID(index);
+                                        // Extract local indices in
+                                        // the vector.
+    int trilinos_i = map->LID(index);
     TrilinosScalar value = 0.;
     if (trilinos_i == -1 )
       {
@@ -388,8 +394,9 @@ namespace TrilinosWrappers
   Vector::real_type
   Vector::lp_norm (const TrilinosScalar p) const
   {
-                                     // get a representation of the vector and
-                                     // loop over all the elements
+                                        // 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);
@@ -399,6 +406,7 @@ namespace TrilinosWrappers
     TrilinosScalar sum=0;
 
     const TrilinosScalar * ptr  = start_ptr;
+
                                        // add up elements
     while (ptr != start_ptr+size())
       sum += std::pow(std::fabs(*ptr++), p);
@@ -686,9 +694,9 @@ namespace TrilinosWrappers
            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
+                                        // 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));
 
@@ -766,9 +774,10 @@ namespace TrilinosWrappers
 
 
   
-                                    // TODO: up to now only local data
-                                     // printed out! Find a way to neatly
-                                    // output distributed data...
+                                       // TODO: up to now only local
+                                        // data printed out! Find a
+                                        // way to neatly output
+                                        // distributed data...
   void
   Vector::print (const char *format) const
   {
@@ -796,11 +805,12 @@ namespace TrilinosWrappers
   {
     AssertThrow (out, ExcIO());
 
-                                     // get a representation of the vector and
-                                     // loop over all the elements 
-                                     // TODO: up to now only local data
-                                     // printed out! Find a way to neatly
-                                    // output distributed data...
+                                        // get a representation of the
+                                        // vector and loop over all
+                                        // the elements TODO: up to
+                                        // now only local data printed
+                                        // out! Find a way to neatly
+                                        // output distributed data...
     TrilinosScalar *val;
     int leading_dimension;
     int ierr = vector->ExtractView (&val, &leading_dimension);
@@ -820,8 +830,8 @@ namespace TrilinosWrappers
         out << static_cast<double>(val[i]) << std::endl;
     out << std::endl;
 
-                                     // restore the representation of the
-                                     // vector
+                                        // restore the representation
+                                        // of the vector
     AssertThrow (out, ExcIO());
   }
 
@@ -830,9 +840,9 @@ namespace TrilinosWrappers
   void
   Vector::swap (Vector &v)
   {
-                                    // Just swap the pointers to the 
-                                    // two Epetra vectors that hold all
-                                    // the data.
+                                        // Just swap the pointers to
+                                        // the two Epetra vectors that
+                                        // hold all the data.
     Vector *p_v = &v, *p_this = this;
     Vector* tmp = p_v;
     p_v = p_this;
@@ -840,6 +850,7 @@ namespace TrilinosWrappers
   }
 
 
+
   unsigned int
   Vector::memory_consumption () const
   {

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.