]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added a function TrilinosWrappers::SparsityPattern that can generate distributed...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 27 Dec 2008 20:42:06 +0000 (20:42 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 27 Dec 2008 20:42:06 +0000 (20:42 +0000)
git-svn-id: https://svn.dealii.org/trunk@18021 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/dofs/dof_constraints.cc
deal.II/deal.II/source/dofs/dof_tools.cc
deal.II/lac/include/lac/block_sparsity_pattern.h
deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/include/lac/trilinos_sparsity_pattern.h [new file with mode: 0755]
deal.II/lac/source/block_sparsity_pattern.cc
deal.II/lac/source/trilinos_block_sparse_matrix.cc
deal.II/lac/source/trilinos_solver_block.cc
deal.II/lac/source/trilinos_sparse_matrix.cc
deal.II/lac/source/trilinos_sparsity_pattern.cc [new file with mode: 0755]

index 242130442447add018e661739af43f984ec24783..80c251d5eec7ca9f6a8203790a12551f567ef146 100644 (file)
@@ -2090,16 +2090,23 @@ add_entries_local_to_global<CompressedSetSparsityPattern> (const std::vector<uns
                                              CompressedSetSparsityPattern       &,
                                              const bool,
                                              const Table<2,bool> &) const;
-
 template void ConstraintMatrix::
 add_entries_local_to_global<CompressedSimpleSparsityPattern> (const std::vector<unsigned int> &,
                                              CompressedSimpleSparsityPattern       &,
                                              const bool,
                                              const Table<2,bool> &) const;
+#ifdef DEAL_II_USE_TRILINOS
+template void ConstraintMatrix::
+add_entries_local_to_global<TrilinosWrappers::SparsityPattern>
+                                             (const std::vector<unsigned int> &,
+                                             TrilinosWrappers::SparsityPattern &,
+                                             const bool,
+                                             const Table<2,bool> &) const;
+#endif
 
 template void ConstraintMatrix::
 add_entries_local_to_global<BlockSparsityPattern> (const std::vector<unsigned int> &,
-                                             BlockSparsityPattern       &,
+                                             BlockSparsityPattern &,
                                              const bool,
                                              const Table<2,bool> &) const;
 template void ConstraintMatrix::
@@ -2109,14 +2116,21 @@ add_entries_local_to_global<BlockCompressedSparsityPattern> (const std::vector<u
                                              const Table<2,bool> &) const;
 template void ConstraintMatrix::
 add_entries_local_to_global<BlockCompressedSetSparsityPattern> (const std::vector<unsigned int> &,
-                                             BlockCompressedSetSparsityPattern       &,
+                                             BlockCompressedSetSparsityPattern &,
                                              const bool,
                                              const Table<2,bool> &) const;
-
 template void ConstraintMatrix::
 add_entries_local_to_global<BlockCompressedSimpleSparsityPattern> (const std::vector<unsigned int> &,
-                                              BlockCompressedSimpleSparsityPattern       &,
+                                              BlockCompressedSimpleSparsityPattern &,
                                               const bool,
                                              const Table<2,bool> &) const;
+#ifdef DEAL_II_USE_TRILINOS
+template void ConstraintMatrix::
+add_entries_local_to_global<TrilinosWrappers::BlockSparsityPattern> 
+                                             (const std::vector<unsigned int> &,
+                                             TrilinosWrappers::BlockSparsityPattern &,
+                                             const bool,
+                                             const Table<2,bool> &) const;
+#endif
 
 DEAL_II_NAMESPACE_CLOSE
index 9bba2f39f0ffa7df5ae705c17cfe7de57c113a9b..7fae1afa9f97c2076da5c9c7dac9f6e39960d0b3 100644 (file)
@@ -17,6 +17,7 @@
 #include <base/quadrature_lib.h>
 #include <base/table.h>
 #include <base/template_constraints.h>
+#include <base/utilities.h>
 #include <grid/tria.h>
 #include <grid/tria_iterator.h>
 #include <grid/intergrid_map.h>
@@ -33,6 +34,7 @@
 #include <lac/compressed_sparsity_pattern.h>
 #include <lac/compressed_set_sparsity_pattern.h>
 #include <lac/compressed_simple_sparsity_pattern.h>
+#include <lac/trilinos_sparsity_pattern.h>
 #include <lac/block_sparsity_pattern.h>
 #include <lac/vector.h>
 #include <numerics/vectors.h>
@@ -66,7 +68,24 @@ DoFTools::make_sparsity_pattern (const DH               &dof,
   dofs_on_this_cell.reserve (max_dofs_per_cell(dof));
   typename DH::active_cell_iterator cell = dof.begin_active(),
                                    endc = dof.end();
+
+                                  // In case we work with a distributed
+                                  // sparsity pattern of Trilinos type, we
+                                  // only have to do the work if the
+                                  // current cell is owned by the calling
+                                  // processor. Otherwise, just continue.
   for (; cell!=endc; ++cell) 
+#ifdef DEAL_II_USE_TRILINOS
+    if ((types_are_equal<SparsityPattern,TrilinosWrappers::SparsityPattern>::value
+        ||
+        types_are_equal<SparsityPattern,TrilinosWrappers::BlockSparsityPattern>::value)
+       &&
+       cell->subdomain_id() != 
+       Utilities::Trilinos::get_this_mpi_process(Utilities::Trilinos::comm_world()))
+      continue;
+    else
+        
+#endif
     {
       const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
       dofs_on_this_cell.resize (dofs_per_cell);
@@ -156,7 +175,24 @@ DoFTools::make_sparsity_pattern (
   std::vector<unsigned int> dofs_on_this_cell(fe_collection.max_dofs_per_cell());
   typename DH::active_cell_iterator cell = dof.begin_active(),
                                    endc = dof.end();
-  for (; cell!=endc; ++cell) 
+
+                                  // In case we work with a distributed
+                                  // sparsity pattern of Trilinos type, we
+                                  // only have to do the work if the
+                                  // current cell is owned by the calling
+                                  // processor. Otherwise, just continue.
+  for (; cell!=endc; ++cell)
+#ifdef DEAL_II_USE_TRILINOS
+    if ((types_are_equal<SparsityPattern,TrilinosWrappers::SparsityPattern>::value
+        ||
+        types_are_equal<SparsityPattern,TrilinosWrappers::BlockSparsityPattern>::value)
+       &&
+       cell->subdomain_id() != 
+       Utilities::Trilinos::get_this_mpi_process(Utilities::Trilinos::comm_world()))
+      continue;
+    else
+        
+#endif
     {
       const unsigned int fe_index
        = cell->active_fe_index();
@@ -5394,6 +5430,16 @@ DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
  CompressedSimpleSparsityPattern &sparsity,
  const ConstraintMatrix &,
  const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
+                               TrilinosWrappers::SparsityPattern>
+(const DoFHandler<deal_II_dimension> &dof,
+ TrilinosWrappers::SparsityPattern    &sparsity,
+ const ConstraintMatrix &,
+ const bool);
+#endif
+
 template void
 DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
                                BlockSparsityPattern>
@@ -5422,6 +5468,15 @@ DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
  BlockCompressedSimpleSparsityPattern      &sparsity,
  const ConstraintMatrix &,
  const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
+                               TrilinosWrappers::BlockSparsityPattern>
+(const DoFHandler<deal_II_dimension> &dof,
+ TrilinosWrappers::BlockSparsityPattern                &sparsity,
+ const ConstraintMatrix &,
+ const bool);
+#endif
 
 template void
 DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
@@ -5430,8 +5485,6 @@ DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
  SparsityPattern    &sparsity,
  const ConstraintMatrix &,
  const bool);
-
-
 template void
 DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
                                CompressedSparsityPattern>
@@ -5439,7 +5492,6 @@ DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
  CompressedSparsityPattern    &sparsity,
  const ConstraintMatrix &,
  const bool);
-
 template void
 DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
                                CompressedSetSparsityPattern>
@@ -5447,7 +5499,6 @@ DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
  CompressedSetSparsityPattern    &sparsity,
  const ConstraintMatrix &,
  const bool);
-
 template void
 DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
                                CompressedSimpleSparsityPattern>
@@ -5455,7 +5506,15 @@ DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
  CompressedSimpleSparsityPattern    &sparsity,
  const ConstraintMatrix &,
  const bool);
-
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
+                               TrilinosWrappers::SparsityPattern>
+(const hp::DoFHandler<deal_II_dimension> &dof,
+ TrilinosWrappers::SparsityPattern    &sparsity,
+ const ConstraintMatrix &,
+ const bool);
+#endif
 
 template void
 DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
@@ -5485,6 +5544,16 @@ DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
  BlockCompressedSimpleSparsityPattern      &sparsity,
  const ConstraintMatrix &,
  const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
+                               TrilinosWrappers::BlockSparsityPattern>
+(const hp::DoFHandler<deal_II_dimension> &dof,
+ TrilinosWrappers::BlockSparsityPattern  &sparsity,
+ const ConstraintMatrix &,
+ const bool);
+#endif
+
 
 
 template void 
@@ -5519,6 +5588,17 @@ DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
  CompressedSimpleSparsityPattern&,
  const ConstraintMatrix &,
  const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void 
+DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
+                               TrilinosWrappers::SparsityPattern>
+(const DoFHandler<deal_II_dimension>&,
+ const Table<2,Coupling>&,
+ TrilinosWrappers::SparsityPattern&,
+ const ConstraintMatrix &,
+ const bool);
+#endif
+
 template void
 DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
                                BlockSparsityPattern>
@@ -5551,6 +5631,16 @@ DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
  BlockCompressedSimpleSparsityPattern&,
  const ConstraintMatrix &,
  const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<DoFHandler<deal_II_dimension>,
+                               TrilinosWrappers::BlockSparsityPattern>
+(const DoFHandler<deal_II_dimension>&,
+ const Table<2,Coupling>&,
+ TrilinosWrappers::BlockSparsityPattern&,
+ const ConstraintMatrix &,
+ const bool);
+#endif
 
 template void 
 DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
@@ -5584,6 +5674,17 @@ DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
  CompressedSimpleSparsityPattern&,
  const ConstraintMatrix &,
  const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void 
+DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
+                               TrilinosWrappers::SparsityPattern>
+(const hp::DoFHandler<deal_II_dimension>&,
+ const Table<2,Coupling>&,
+ TrilinosWrappers::SparsityPattern&,
+ const ConstraintMatrix &,
+ const bool);
+#endif
+
 template void
 DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
                                BlockSparsityPattern>
@@ -5616,6 +5717,17 @@ DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
  BlockCompressedSimpleSparsityPattern&,
  const ConstraintMatrix &,
  const bool);
+#ifdef DEAL_II_USE_TRILINOS
+template void
+DoFTools::make_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,
+                               TrilinosWrappers::BlockSparsityPattern>
+(const hp::DoFHandler<deal_II_dimension>&,
+ const Table<2,Coupling>&,
+ TrilinosWrappers::BlockSparsityPattern&,
+ const ConstraintMatrix &,
+ const bool);
+#endif
+
 
 
 template void
index 2952413051b0bf17fed0b48cd7204ba991cfd612..45b914cba8e3af9b50b209fb28b04dde75914a71 100644 (file)
@@ -20,6 +20,7 @@
 #include <base/subscriptor.h>
 #include <base/smartpointer.h>
 #include <lac/sparsity_pattern.h>
+#include <lac/trilinos_sparsity_pattern.h>
 #include <lac/compressed_sparsity_pattern.h>
 #include <lac/compressed_set_sparsity_pattern.h>
 #include <lac/compressed_simple_sparsity_pattern.h>
@@ -33,6 +34,12 @@ class BlockSparsityPattern;
 class BlockCompressedSparsityPattern;
 class BlockCompressedSimpleSparsityPattern;
 class BlockCompressedSetSparsityPattern;
+#ifdef DEAL_II_USE_TRILINOS
+namespace TrilinosWrappers
+{
+  class BlockSparsityPattern;
+}
+#endif
 
 /*! @addtogroup Sparsity
  *@{
@@ -687,8 +694,6 @@ typedef BlockCompressedSparsityPattern CompressedBlockSparsityPattern;
  * There are several, exchangeable variations of this class, see @ref Sparsity,
  * section 'Dynamic block sparsity patterns' for more information.
  *
- * This class is used in @ref step_22 "step-22".
- *
  * @author Wolfgang Bangerth, 2007
  */
 class BlockCompressedSetSparsityPattern : public BlockSparsityPatternBase<CompressedSetSparsityPattern>
@@ -776,6 +781,8 @@ class BlockCompressedSetSparsityPattern : public BlockSparsityPatternBase<Compre
  * There are several, exchangeable variations of this class, see @ref Sparsity,
  * section 'Dynamic block sparsity patterns' for more information.
  *
+ * This class is used in @ref step_22 "step-22" and @ref step_31 "step-31".
+ *
  * @author Timo Heister, 2008
  */
 class BlockCompressedSimpleSparsityPattern : public BlockSparsityPatternBase<CompressedSimpleSparsityPattern>
@@ -803,7 +810,7 @@ class BlockCompressedSimpleSparsityPattern : public BlockSparsityPatternBase<Com
                                      * you assign them sizes.
                                      */
     BlockCompressedSimpleSparsityPattern (const unsigned int n_rows,
-                                   const unsigned int n_columns);
+                                         const unsigned int n_columns);
 
                                     /**
                                      * Initialize the pattern with
@@ -847,6 +854,129 @@ class BlockCompressedSimpleSparsityPattern : public BlockSparsityPatternBase<Com
 
 
 
+
+#ifdef DEAL_II_USE_TRILINOS
+
+
+/**
+ * This class extends the base class to implement an array of Trilinos
+ * sparsity patterns that can be used to initialize Trilinos block sparse
+ * matrices that can be distributed among different processors. It is used
+ * in the same way as the BlockSparsityPattern except that it builds upon
+ * the TrilinosWrappers::SparsityPattern instead of the
+ * dealii::SparsityPattern. See the documentation of the
+ * BlockSparsityPattern for examples.
+ *
+ * This class is has properties of the "dynamic" type of @ref Sparsity (in
+ * the sense that it can extend the memory if too little elements were
+ * allocated), but otherwise is more like the basic deal.II SparsityPattern
+ * (in the sense that the method compress() needs to be called before the
+ * pattern can be used).
+ *
+ * This class is used in @ref step_32 "step-32".
+ *
+ * @author Martin Kronbichler, 2008
+ */
+namespace TrilinosWrappers
+{
+  class BlockSparsityPattern : 
+    public dealii::BlockSparsityPatternBase<SparsityPattern>
+  {
+    public:
+
+                                    /**
+                                     * Initialize the matrix empty,
+                                     * that is with no memory
+                                     * allocated. This is useful if
+                                     * you want such objects as
+                                     * member variables in other
+                                     * classes. You can make the
+                                     * structure usable by calling
+                                     * the reinit() function.
+                                     */
+      BlockSparsityPattern ();
+
+                                    /**
+                                     * Initialize the matrix with the
+                                     * given number of block rows and
+                                     * columns. The blocks themselves
+                                     * are still empty, and you have
+                                     * to call collect_sizes() after
+                                     * you assign them sizes.
+                                     */
+      BlockSparsityPattern (const unsigned int n_rows,
+                           const unsigned int n_columns);
+
+                                    /**
+                                     * Initialize the pattern with
+                                     * two BlockIndices for the block
+                                     * structures of matrix rows and
+                                     * columns. This function is
+                                     * equivalent to calling the
+                                     * previous constructor with the
+                                     * length of the two index vector
+                                     * and then entering the index
+                                     * values.
+                                     */
+      BlockSparsityPattern (const std::vector<unsigned int>& row_block_sizes,
+                           const std::vector<unsigned int>& col_block_sizes);
+
+                                    /**
+                                     * Initialize the pattern with an array
+                                     * Epetra_Map that specifies both rows
+                                     * and columns of the matrix (so the
+                                     * final matrix will be a square
+                                     * matrix), where the Epetra_Map
+                                     * specifies the parallel distribution
+                                     * of the degrees of freedom on the
+                                     * individual block.  This function is
+                                     * equivalent to calling the second
+                                     * constructor with the length of the
+                                     * mapping vector and then entering the
+                                     * index values.
+                                     */
+      BlockSparsityPattern (const std::vector<Epetra_Map>& input_maps);
+
+                                    /**
+                                     * Resize the matrix to a tensor
+                                     * product of matrices with
+                                     * dimensions defined by the
+                                     * arguments.
+                                     *
+                                     * The matrix will have as many
+                                     * block rows and columns as
+                                     * there are entries in the two
+                                     * arguments. The block at
+                                     * position (<i>i,j</i>) will
+                                     * have the dimensions
+                                     * <tt>row_block_sizes[i]</tt>
+                                     * times <tt>col_block_sizes[j]</tt>.
+                                     */
+      void reinit (const std::vector< unsigned int > &row_block_sizes,
+                  const std::vector< unsigned int > &col_block_sizes);
+
+                                    /**
+                                     * Resize the matrix to a square tensor
+                                     * product of matrices with parallel
+                                     * distribution according to the
+                                     * specifications in the array of
+                                     * Epetra_Maps.
+                                     */
+      void reinit (const std::vector<Epetra_Map>& input_maps);
+
+
+                                    /**
+                                     * Allow the use of the reinit
+                                     * functions of the base class as
+                                     * well.
+                                     */
+    using BlockSparsityPatternBase<SparsityPattern>::reinit;
+  };
+}
+
+#endif
+
+
 /*@}*/
 /*---------------------- Template functions -----------------------------------*/
 
index 21e9cf65e22e72996cdd9137d5c161e8918fef44..0a646f832572730fb254267e835563954fec85ef 100755 (executable)
 
 DEAL_II_NAMESPACE_OPEN
 
-                                   // forward declarations
-class SparsityPattern;
-class CompressedSparsityPattern;
-class CompressedSetSparsityPattern;
-class CompressedSimpleSparsityPattern;
 
 namespace TrilinosWrappers
 {
                                    // forward declarations
   class VectorBase;
   class SparseMatrix;
+  class SparsityPattern;
 
   namespace MatrixIterators
   {
@@ -60,10 +56,11 @@ namespace TrilinosWrappers
  * 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 stores the elements within each row in ascending
+ * order. This is opposed to the deal.II sparse matrix style where the
+ * diagonal element (if it exists) is stored before all other values, and
+ * the PETSc sparse matrices, where one can't guarantee a certain order of
+ * the elements.
  *
  * @ingroup TrilinosWrappers
  * @author Martin Kronbichler, Wolfgang Bangerth, 2008
@@ -83,9 +80,9 @@ namespace TrilinosWrappers
                                        * access, a const matrix
                                        * pointer is sufficient.
                                        */
-            Accessor (const SparseMatrix  *matrix,
-                      const unsigned int   row,
-                      const unsigned int   index);
+            Accessor (const SparseMatrix *matrix,
+                      const unsigned int  row,
+                      const unsigned int  index);
 
                                        /**
                                        * Row number of the element
@@ -199,7 +196,7 @@ namespace TrilinosWrappers
                                        * matrix for the given row and
                                        * the index within it.
                                        */ 
-        const_iterator (const SparseMatrix   *matrix,
+        const_iterator (const SparseMatrix *matrix,
                         const unsigned int  row,
                         const unsigned int  index);
           
@@ -344,22 +341,23 @@ 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
+                                       * setup. However, there should be no
+                                       * effect in the performance of
+                                       * matrix-vector products, since
+                                       * Trilinos wants to reorganize the
+                                       * matrix memory prior to use.
                                         */
       SparseMatrix (const Epetra_Map   &InputMap,
                    const unsigned int  n_max_entries_per_row);
@@ -383,29 +381,28 @@ namespace TrilinosWrappers
                    const std::vector<unsigned int> &n_entries_per_row);
 
                                        /**
-                                        * This constructor 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 takes
-                                        * care of the columns and the
-                                        * other one of the rows. Note
-                                        * that there is no real
-                                        * parallelism along the columns
-                                        * &ndash; the processor that
-                                        * owns a certain row always owns
-                                        * all the column elements, no
-                                        * matter how far they might be
-                                        * spread out. The second
-                                        * Epetra_Map is only used to
-                                        * specify the number of columns
-                                        * and for internal arragements
-                                        * when doing matrix-vector
-                                        * products with vectors based on
-                                        * that column map.
+                                        * This constructor 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 describes the parallel
+                                        * partitioning of the dofs
+                                        * associated with the matrix rows
+                                        * and the other one the partitioning
+                                        * of dofs in the matrix
+                                        * columns. Note that there is no
+                                        * real parallelism along the columns
+                                        * &ndash; the processor that owns a
+                                        * certain row always owns all the
+                                        * column elements, no matter how far
+                                        * they might be spread out. The
+                                        * second Epetra_Map is only used to
+                                        * specify the number of columns and
+                                        * for internal arragements when
+                                        * doing matrix-vector products with
+                                        * vectors based on that column map.
                                        *
                                        * The number of columns entries
                                        * per row is specified as the
@@ -424,16 +421,22 @@ namespace TrilinosWrappers
                                        * to be used for generating
                                        * rectangular matrices, where one
                                        * map specifies the parallel
-                                       * distribution of rows and the
-                                       * second one specifies the number of
-                                       * columns in the total matrix. It
-                                       * also provides information for the
+                                       * distribution of degrees of freedom
+                                       * associated with matrix rows and
+                                       * the second one specifies the
+                                       * parallel distribution the dofs
+                                       * associated with columns in the
+                                       * matrix. The second map also
+                                       * provides information for the
                                        * internal arrangement in matrix
-                                       * vector products, but is not used
-                                       * for the distribution of the
-                                       * columns &ndash; rather, all column
-                                       * elements of a row are stored on
-                                       * the same processor. The vector
+                                       * vector products (i.e., the
+                                       * distribution of vector this matrix
+                                       * is to be multiplied with), but is
+                                       * not used for the distribution of
+                                       * the columns &ndash; rather, all
+                                       * column elements of a row are
+                                       * stored on the same processor in
+                                       * any case. The vector
                                        * <tt>n_entries_per_row</tt>
                                        * specifies the number of entries in
                                        * each row of the newly generated
@@ -472,6 +475,12 @@ namespace TrilinosWrappers
                    const unsigned int               n,
                    const std::vector<unsigned int> &n_entries_per_row);
 
+                                       /**
+                                        * Generate a matrix from a Trilinos
+                                        * sparsity pattern object.
+                                        */
+      SparseMatrix (const SparsityPattern &InputSparsityPattern);
+
                                        /**
                                         * Copy constructor. Sets the
                                         * calling matrix to be the same
@@ -489,75 +498,64 @@ namespace TrilinosWrappers
       virtual ~SparseMatrix ();
 
                                        /**
-                                        * 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
+                                        * 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. This function is
-                                       * meant for use in serial
-                                       * programs, where there is no
-                                       * need to specify how the
-                                       * matrix is going to be
-                                       * distributed among the
-                                       * processors. This function
-                                       * works in parallel, too, but
-                                       * it is recommended to
-                                       * manually specify the
-                                       * parallel partioning of the
-                                       * matrix using an
-                                       * Epetra_Map. When run in
-                                       * parallel, it is currently
-                                       * necessary that each
-                                       * processor holds the
-                                       * sparsity_pattern structure
-                                       * because each processor sets
-                                       * its rows.
+                                       * pattern. This function is meant
+                                       * for use in serial programs, where
+                                       * there is no need to specify how
+                                       * the matrix is going to be
+                                       * distributed among different
+                                       * processors. This function works in
+                                       * parallel, too, but it is
+                                       * recommended to manually specify
+                                       * the parallel partioning of the
+                                       * matrix using an Epetra_Map. When
+                                       * run in parallel, it is currently
+                                       * necessary that each processor
+                                       * holds the sparsity_pattern
+                                       * structure because each processor
+                                       * sets its rows.
                                        *
-                                       * This is a collective
-                                       * operation that needs to be
-                                       * called on all processors in
-                                       * order to avoid a dead lock.
+                                       * This is a collective operation
+                                       * that needs to be called on all
+                                       * processors in order to avoid a
+                                       * dead lock.
                                         */
       template<typename SparsityType>
       void reinit (const SparsityType &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
-                                       * according to a user-supplied
-                                       * Epetra map. In programs
-                                       * following the style of the
-                                       * tutorial programs, this
-                                       * function (and the respective
-                                       * call for a rectangular matrix)
-                                       * are the natural way to
-                                       * initialize the matrix size,
-                                       * its distribution among the MPI
-                                       * processes (if run in parallel)
-                                       * as well as the locatoin of
-                                       * non-zero elements. Trilinos
+                                        * 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. In
+                                       * programs following the style of
+                                       * the tutorial programs, this
+                                       * function (and the respective call
+                                       * for a rectangular matrix) are the
+                                       * natural way to initialize the
+                                       * matrix size, its distribution
+                                       * among the MPI processes (if run in
+                                       * parallel) as well as the locatoin
+                                       * of non-zero elements. Trilinos
                                        * stores the sparsity pattern
-                                       * internally, so it won't be
-                                       * needed any more after this
-                                       * call, in contrast to the
-                                       * deal.II own object. In a
-                                       * parallel run, it is currently
-                                       * necessary that each processor
-                                       * holds the sparsity_pattern
-                                       * structure because each
-                                       * processor sets its
-                                       * rows.
+                                       * internally, so it won't be needed
+                                       * any more after this call, in
+                                       * contrast to the deal.II own
+                                       * object. In a parallel run, it is
+                                       * currently necessary that each
+                                       * processor holds the
+                                       * sparsity_pattern structure because
+                                       * each processor sets its rows.
                                        *
-                                       * This is a
-                                       * collective operation that
-                                       * needs to be called on all
+                                       * This is a collective operation
+                                       * that needs to be called on all
                                        * processors in order to avoid a
                                        * dead lock.
                                         */
@@ -566,18 +564,16 @@ namespace TrilinosWrappers
                   const SparsityType  &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.
+                                        * 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.
                                        *
-                                       * This is a
-                                       * collective operation that
-                                       * needs to be called on all
+                                       * This is a collective operation
+                                       * that needs to be called on all
                                        * processors in order to avoid a
                                        * dead lock.
                                         */
@@ -587,34 +583,42 @@ namespace TrilinosWrappers
                   const SparsityType  &sparsity_pattern);
 
                                       /**
-                                       * This function copies the
-                                       * content in
-                                       * <tt>sparse_matrix</tt> to
-                                       * the calling matrix.
+                                       * This function reinitializes the
+                                       * Trilinos sparse matrix from a
+                                       * (possibly distributed) Trilinos
+                                       * sparsity pattern.
+                                       *
+                                       * This is a collective operation
+                                       * that needs to be called on all
+                                       * processors in order to avoid a
+                                       * dead lock.
+                                       */
+      void reinit (const SparsityPattern &sparsity_pattern);
+
+                                      /**
+                                       * This function copies the content
+                                       * in <tt>sparse_matrix</tt> to the
+                                       * calling matrix.
                                        *
-                                       * This is a
-                                       * collective operation that
-                                       * needs to be called on all
+                                       * This is a collective operation
+                                       * that needs to be called on all
                                        * processors in order to avoid a
                                        * dead lock.
                                        */
       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).
                                        *
-                                       * This is a
-                                       * collective operation that
-                                       * needs to be called on all
+                                       * This is a collective operation
+                                       * that needs to be called on all
                                        * processors in order to avoid a
                                        * dead lock.
                                         */
@@ -622,27 +626,23 @@ namespace TrilinosWrappers
                   const double                          drop_tolerance=1e-13);
 
                                       /**
-                                        * 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). In contrast
-                                       * to the other reinit function
-                                       * with deal.II sparse matrix
-                                       * argument, this function
-                                       * takes a parallel
-                                       * partitioning specified by
-                                       * the user instead of
+                                        * 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). In
+                                       * contrast to the other reinit
+                                       * function with deal.II sparse
+                                       * matrix argument, this function
+                                       * takes a parallel partitioning
+                                       * specified by the user instead of
                                        * internally generating one.
                                        *
-                                       * This is a
-                                       * collective operation that
-                                       * needs to be called on all
+                                       * This is a collective operation
+                                       * that needs to be called on all
                                        * processors in order to avoid a
                                        * dead lock.
                                         */
@@ -651,18 +651,16 @@ namespace TrilinosWrappers
                   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. Chosen 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. Chosen for rectangular
+                                       * matrices.
                                        *
-                                       * This is a
-                                       * collective operation that
-                                       * needs to be called on all
+                                       * This is a collective operation
+                                       * that needs to be called on all
                                        * processors in order to avoid a
                                        * dead lock.
                                         */
@@ -679,58 +677,51 @@ namespace TrilinosWrappers
       void reinit (const Epetra_CrsMatrix &input_matrix);
 
                                        /**
-                                        * 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.
+                                        * This operator assigns a scalar to
+                                        * a matrix. Since this does usually
+                                        * not make much sense (should we set
+                                        * all matrix entries to this value?
+                                        * Only the nonzero entries of the
+                                        * sparsity pattern?), this operation
+                                        * is only allowed if the actual
+                                        * value to be assigned is zero. This
+                                        * operator only exists to allow for
+                                        * the obvious notation
+                                        * <tt>matrix=0</tt>, which sets all
+                                        * elements of the matrix to zero,
+                                        * but keeps the sparsity pattern
+                                        * previously used.
                                         */
       SparseMatrix &
        operator = (const double d);
 
                                        /**
-                                        * 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.
                                        *
-                                       * This is a
-                                       * collective operation that
-                                       * needs to be called on all
+                                       * This is a collective operation
+                                       * that needs to be called on all
                                        * processors in order to avoid a
                                        * dead lock.
                                         */
       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 assembled
-                                        * the matrix. This is a
+                                        * 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. This is a
                                         * collective operation, i.e., it
-                                        * needs to be run on all
-                                        * processors when used in
-                                        * parallel.
+                                        * needs to be run on all processors
+                                        * when used in parallel.
                                         */
       void compress ();
 
@@ -739,9 +730,9 @@ namespace TrilinosWrappers
                                        * i.e., whether compress() needs to
                                        * be called after an operation
                                        * requiring data exchange. A call to
-                                       * compress() is also after the
-                                       * method set() is called (even when
-                                       * working in serial).
+                                       * compress() is also needed when the
+                                       * method set() has been called (even
+                                       * when working in serial).
                                        */
       bool is_compressed () const;
 //@}
@@ -1741,7 +1732,7 @@ namespace TrilinosWrappers
                                        * to the individual processes.
                                        */
       Epetra_Map col_map;
-  
+
                                        /**
                                         * Trilinos doesn't allow to
                                         * mix additions to matrix
@@ -2476,7 +2467,7 @@ DEAL_II_NAMESPACE_CLOSE
 #endif // DEAL_II_USE_TRILINOS
 
 
-/*----------------------------   trilinos_sparse_matrix.h     ---------------------------*/
+/*-----------------------   trilinos_sparse_matrix.h     --------------------*/
 
 #endif
-/*----------------------------   trilinos_sparse_matrix.h     ---------------------------*/
+/*-----------------------   trilinos_sparse_matrix.h     --------------------*/
diff --git a/deal.II/lac/include/lac/trilinos_sparsity_pattern.h b/deal.II/lac/include/lac/trilinos_sparsity_pattern.h
new file mode 100755 (executable)
index 0000000..11707a9
--- /dev/null
@@ -0,0 +1,1272 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2008 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__trilinos_sparsity_pattern_h
+#define __deal2__trilinos_sparsity_pattern_h
+
+
+#include <base/config.h>
+#include <base/subscriptor.h>
+#include <lac/exceptions.h>
+
+#include <boost/shared_ptr.hpp>
+#include <vector>
+#include <cmath>
+#include <memory>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+#  include <Epetra_FECrsGraph.h>
+#  include <Epetra_Map.h>
+#  ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+#    include <Epetra_MpiComm.h>
+#    include "mpi.h"
+#  else
+#    include "Epetra_SerialComm.h"
+#  endif
+
+DEAL_II_NAMESPACE_OPEN
+
+                                   // forward declarations
+class SparsityPattern;
+class CompressedSparsityPattern;
+class CompressedSetSparsityPattern;
+class CompressedSimpleSparsityPattern;
+
+namespace TrilinosWrappers
+{
+                                   // forward declarations
+  class SparsityPattern;
+
+  namespace SparsityPatternIterators
+  {
+/**
+ * STL conforming iterator. This class acts as an iterator walking
+ * over the elements of Trilinos sparsity pattern.
+ *
+ * @ingroup TrilinosWrappers
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
+ */
+    class const_iterator
+    {
+      private:
+                                       /**
+                                       * Accessor class for iterators
+                                       */
+        class Accessor
+        {
+          public:
+                                       /**
+                                       * Constructor. Since we use
+                                       * accessors only for read
+                                       * access, a const matrix
+                                       * pointer is sufficient.
+                                       */
+            Accessor (const SparsityPattern *sparsity_pattern,
+                      const unsigned int     row,
+                      const unsigned int     index);
+
+                                       /**
+                                       * Row number of the element
+                                       * represented by this object.
+                                       */
+            unsigned int row() const;
+
+                                       /**
+                                       * Index in row of the element
+                                       * represented by this object.
+                                       */
+            unsigned int index() const;
+
+                                       /**
+                                       * Column number of the element
+                                       * represented by this object.
+                                       */
+            unsigned int column() const;
+
+                                       /**
+                                       * Exception
+                                       */
+            DeclException0 (ExcBeyondEndOfSparsityPattern);
+
+                                      /**
+                                       * Exception
+                                       */
+            DeclException3 (ExcAccessToNonlocalRow,
+                            int, int, int,
+                            << "You tried to access row " << arg1
+                            << " of a distributed sparsity pattern, "
+                           << " but only rows " << arg2 << " through " << arg3
+                            << " are stored locally and can be accessed.");
+            
+          private:
+                                       /**
+                                       * The matrix accessed.
+                                       */
+            mutable SparsityPattern *sparsity_pattern;
+
+                                       /**
+                                       * Current row number.
+                                       */
+            unsigned int a_row;
+
+                                       /**
+                                       * Current index in row.
+                                       */
+            unsigned int a_index;
+
+                                       /**
+                                       * Cache where we store the
+                                       * column indices of the
+                                       * present row. This is
+                                       * necessary, since Trilinos
+                                       * makes access to the elements
+                                       * of its matrices rather hard,
+                                       * and it is much more
+                                       * efficient to copy all column
+                                       * entries of a row once when
+                                       * we enter it than repeatedly
+                                       * asking Trilinos for
+                                       * individual ones. This also
+                                       * makes some sense since it is
+                                       * likely that we will access
+                                       * them sequentially anyway.
+                                       *
+                                       * In order to make copying of
+                                       * iterators/accessor of
+                                       * acceptable performance, we
+                                       * keep a shared pointer to
+                                       * these entries so that more
+                                       * than one accessor can access
+                                       * this data if necessary.
+                                       */
+            boost::shared_ptr<const std::vector<unsigned int> > colnum_cache;
+            
+                                      /**
+                                       * Discard the old row caches
+                                       * (they may still be used by
+                                       * other accessors) and
+                                       * generate new ones for the
+                                       * row pointed to presently by
+                                       * this accessor.
+                                       */
+            void visit_present_row ();
+
+                                       /**
+                                       * Make enclosing class a
+                                       * friend.
+                                       */
+            friend class const_iterator;
+        };
+        
+      public:
+          
+                                       /**
+                                       * Constructor. Create an
+                                       * iterator into the matrix @p
+                                       * matrix for the given row and
+                                       * the index within it.
+                                       */ 
+        const_iterator (const SparsityPattern *sparsity_pattern,
+                        const unsigned int     row,
+                        const unsigned int     index);
+          
+                                       /**
+                                       * Prefix increment.
+                                       */
+        const_iterator& operator++ ();
+
+                                       /**
+                                       * Postfix increment.
+                                       */
+        const_iterator operator++ (int);
+
+                                       /**
+                                       * Dereferencing operator.
+                                       */
+        const Accessor& operator* () const;
+
+                                       /**
+                                       * Dereferencing operator.
+                                       */
+        const Accessor* operator-> () const;
+
+                                       /**
+                                       * Comparison. True, if both
+                                       * iterators point to the same
+                                       * matrix position.
+                                       */
+        bool operator == (const const_iterator&) const;
+
+                                       /**
+                                       * Inverse of <tt>==</tt>.
+                                       */
+        bool operator != (const const_iterator&) const;
+
+                                       /**
+                                       * Comparison operator. Result
+                                       * is true if either the first
+                                       * row number is smaller or if
+                                       * the row numbers are equal
+                                       * and the first index is
+                                       * smaller.
+                                       */
+        bool operator < (const const_iterator&) const;
+
+                                      /**
+                                        * Exception
+                                       */
+        DeclException2 (ExcInvalidIndexWithinRow,
+                        int, int,
+                        << "Attempt to access element " << arg2
+                        << " of row " << arg1
+                        << " which doesn't have that many elements.");
+        
+      private:
+                                       /**
+                                        * Store an object of the
+                                        * accessor class.
+                                       */
+        Accessor accessor;
+    };
+    
+  }
+  
+  
+/**
+ * This class implements a wrapper class to use the Trilinos distributed
+ * sparsity pattern class Epetra_FECrsGraph. This class is designed to be
+ * used for construction of parallel Trilinos matrices. The functionality of
+ * this class is modeled after the existing sparsity pattern classes, with
+ * the difference that this class can work fully in parallel according to a
+ * partitioning of the sparsity pattern rows.
+ *
+ * This class has many similarities to the compressed sparsity pattern
+ * classes of deal.II (i.e., the classes CompressedSparsityPattern,
+ * CompressedSetSparsityPattern, and CompressedSimpleSparsityPattern), since
+ * it can dynamically add elements to the pattern without any memory being
+ * previously reserved for it. However, it also has a method
+ * SparsityPattern::compress(), that finalizes the pattern and enables its
+ * use with Trilinos sparse matrices.
+ * 
+ * @ingroup TrilinosWrappers
+ * @ingroup Sparsity
+ * @author Martin Kronbichler, 2008
+ */
+  class SparsityPattern : public Subscriptor
+  {
+    public:
+
+                                       /**
+                                        * Declare a typedef for the
+                                        * iterator class.
+                                        */
+      typedef SparsityPatternIterators::const_iterator const_iterator;
+      
+/**
+ * @name Constructors and initalization.
+ */
+//@{
+                                       /**
+                                        * Default constructor. Generates an
+                                        * empty (zero-size) sparsity
+                                        * pattern.
+                                        */
+      SparsityPattern ();
+
+                                       /**
+                                        * Constructor for a square sparsity
+                                       * pattern using an Epetra_Map and
+                                       * the number of nonzero entries in
+                                       * the rows of the sparsity
+                                       * pattern. Note that this number
+                                       * does not need to be exact, and it
+                                       * is even allowed that the actual
+                                       * sparsity 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 when creating the
+                                       * sparsity pattern.
+                                        */
+      SparsityPattern (const Epetra_Map   &InputMap,
+                      const unsigned int  n_entries_per_row = 1);
+
+                                       /**
+                                        * Same as before, but now use the
+                                       * exact number of nonzeros in each m
+                                       * row. Since we know the number of
+                                       * elements in the sparsity pattern
+                                       * exactly in this case, we can
+                                       * already allocate the right amount
+                                       * of memory, which makes the
+                                       * creation process by the respective
+                                       * SparsityPattern::reinit call
+                                       * considerably faster. However, this
+                                       * is a rather unusual situation,
+                                       * since knowing the number of
+                                       * entries in each row is usually
+                                       * connected to knowing the indices
+                                       * of nonzero entries, which the
+                                       * sparsity pattern is designed to
+                                       * describe.
+                                        */
+      SparsityPattern (const Epetra_Map                &InputMap,
+                      const std::vector<unsigned int> &n_entries_per_row);
+
+                                       /**
+                                        * This constructor 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 sparsity pattern,
+                                        * where one map describes the
+                                        * parallel partitioning of the dofs
+                                        * associated with the sparsity
+                                        * pattern rows and the other one of
+                                        * the sparsity pattern columns. Note
+                                        * that there is no real parallelism
+                                        * along the columns &ndash; the
+                                        * processor that owns a certain row
+                                        * always owns all the column
+                                        * elements, no matter how far they
+                                        * might be spread out. The second
+                                        * Epetra_Map is only used to specify
+                                        * the number of columns and for
+                                        * internal arragements when doing
+                                        * matrix-vector products with
+                                        * vectors based on that column map.
+                                       *
+                                       * The number of columns entries
+                                       * per row is specified as the
+                                       * maximum number of entries
+                                       * argument.
+                                        */
+      SparsityPattern (const Epetra_Map   &InputRowMap,
+                      const Epetra_Map   &InputColMap,
+                      const unsigned int  n_entries_per_row = 1);
+
+                                       /**
+                                        * This constructor 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
+                                       * second one specifies the
+                                       * distribution of degrees of freedom
+                                       * associated with matrix
+                                       * columns. This second map is
+                                       * however not used for the
+                                       * distribution of the columns
+                                       * themselves &ndash; rather, all
+                                       * column elements of a row are
+                                       * stored on the same processor. The
+                                       * vector <tt>n_entries_per_row</tt>
+                                       * specifies the number of entries in
+                                       * each row of the newly generated
+                                       * matrix.
+                                        */
+      SparsityPattern (const Epetra_Map                &InputRowMap,
+                      const Epetra_Map                &InputColMap,
+                      const std::vector<unsigned int> &n_entries_per_row);
+
+                                       /**
+                                        * Generate a sparsity pattern that
+                                        * is completely stored locally,
+                                        * having #m rows and #n columns. The
+                                        * resulting matrix will be
+                                        * completely stored locally.
+                                       *
+                                       * The number of columns entries per
+                                       * row is specified as the maximum
+                                       * number of entries argument. As
+                                       * above, this does not need to be an
+                                       * accurate number since the entries
+                                       * are allocated dynamically in a
+                                       * similar manner as for the deal.II
+                                       * CompressedSparsityPattern classes,
+                                       * but a good estimate will reduce
+                                       * the setup time of the sparsity
+                                       * pattern.
+                                        */
+      SparsityPattern (const unsigned int  m,
+                      const unsigned int  n,
+                      const unsigned int  n_entries_per_row = 1);
+
+                                       /**
+                                        * Generate a sparsity pattern that
+                                        * is completely stored locally,
+                                        * having #m rows and #n columns. The
+                                        * resulting matrix will be
+                                        * completely stored locally.
+                                       *
+                                       * The vector
+                                       * <tt>n_entries_per_row</tt>
+                                       * specifies the number of entries in
+                                       * each row.
+                                        */
+      SparsityPattern (const unsigned int               m,
+                      const unsigned int               n,
+                      const std::vector<unsigned int> &n_entries_per_row);
+
+                                       /**
+                                        * Copy constructor. Sets the calling
+                                        * sparsity pattern to be the same as
+                                        * the input sparsity pattern.
+                                        */
+      SparsityPattern (const SparsityPattern &SP);
+
+                                       /**
+                                        * Destructor. Made virtual so that
+                                        * one can use pointers to this
+                                        * class.
+                                        */
+      virtual ~SparsityPattern ();
+
+                                       /**
+                                        * Reinitialization function for
+                                       * generating a square sparsity
+                                       * pattern using an Epetra_Map and
+                                       * the number of nonzero entries in
+                                       * the rows of the sparsity
+                                       * pattern. Note that this number
+                                       * does not need to be exact, and it
+                                       * is even allowed that the actual
+                                       * sparsity 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 when creating the
+                                       * sparsity pattern.
+                                       *
+                                       * This function does not create any
+                                       * entries by itself, but provides
+                                       * the correct data structures that
+                                       * can be used by the respective
+                                       * add() function.
+                                        */
+      void  
+      reinit (const Epetra_Map   &InputMap,
+             const unsigned int  n_entries_per_row = 1);
+
+                                       /**
+                                        * Same as before, but now use the
+                                       * exact number of nonzeros in each m
+                                       * row. Since we know the number of
+                                       * elements in the sparsity pattern
+                                       * exactly in this case, we can
+                                       * already allocate the right amount
+                                       * of memory, which makes process of
+                                       * adding entries to the sparsity
+                                       * pattern considerably
+                                       * faster. However, this is a rather
+                                       * unusual situation, since knowing
+                                       * the number of entries in each row
+                                       * is usually connected to knowing
+                                       * the indices of nonzero entries,
+                                       * which the sparsity pattern is
+                                       * designed to describe.
+                                        */
+      void  
+      reinit (const Epetra_Map                &InputMap,
+             const std::vector<unsigned int> &n_entries_per_row);
+
+                                       /**
+                                        * This reinit 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 sparsity pattern,
+                                        * where one map describes the
+                                        * parallel partitioning of the dofs
+                                        * associated with the sparsity
+                                        * pattern rows and the other one of
+                                        * the sparsity pattern columns. Note
+                                        * that there is no real parallelism
+                                        * along the columns &ndash; the
+                                        * processor that owns a certain row
+                                        * always owns all the column
+                                        * elements, no matter how far they
+                                        * might be spread out. The second
+                                        * Epetra_Map is only used to specify
+                                        * the number of columns and for
+                                        * internal arragements when doing
+                                        * matrix-vector products with
+                                        * vectors based on that column map.
+                                       *
+                                       * The number of columns entries per
+                                       * row is specified by the argument
+                                       * <tt>n_entries_per_row</tt>.
+                                        */
+      void  
+      reinit (const Epetra_Map   &InputRowMap,
+             const Epetra_Map   &InputColMap,
+             const unsigned int  n_entries_per_row = 1);
+
+                                       /**
+                                        * This reinit 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
+                                       * second one specifies the
+                                       * distribution of degrees of freedom
+                                       * associated with matrix
+                                       * columns. This second map is
+                                       * however not used for the
+                                       * distribution of the columns
+                                       * themselves &ndash; rather, all
+                                       * column elements of a row are
+                                       * stored on the same processor. The
+                                       * vector <tt>n_entries_per_row</tt>
+                                       * specifies the number of entries in
+                                       * each row of the newly generated
+                                       * matrix.
+                                        */
+      void  
+      reinit (const Epetra_Map                &InputRowMap,
+             const Epetra_Map                &InputColMap,
+             const std::vector<unsigned int> &n_entries_per_row);
+
+                                       /**
+                                        * Initialize a sparsity pattern that
+                                        * is completely stored locally,
+                                        * having #m rows and #n columns. The
+                                        * resulting matrix will be
+                                        * completely stored locally.
+                                       *
+                                       * The number of columns entries per
+                                       * row is specified as the maximum
+                                       * number of entries argument. As
+                                       * above, this does not need to be an
+                                       * accurate number since the entries
+                                       * are allocated dynamically in a
+                                       * similar manner as for the deal.II
+                                       * CompressedSparsityPattern classes,
+                                       * but a good estimate will reduce
+                                       * the setup time of the sparsity
+                                       * pattern.
+                                        */
+      void  
+      reinit (const unsigned int  m,
+             const unsigned int  n,
+             const unsigned int  n_entries_per_row = 1);
+
+                                       /**
+                                        * Initialize a sparsity pattern that
+                                        * is completely stored locally,
+                                        * having #m rows and #n columns. The
+                                        * resulting matrix will be
+                                        * completely stored locally.
+                                       *
+                                       * The vector
+                                       * <tt>n_entries_per_row</tt>
+                                       * specifies the number of entries in
+                                       * each row.
+                                        */
+      void  
+      reinit (const unsigned int               m,
+             const unsigned int               n,
+             const std::vector<unsigned int> &n_entries_per_row);
+
+                                       /**
+                                        * Reinit function. Takes one of the
+                                        * deal.II sparsity patterns and a
+                                        * parallel partitioning of the rows
+                                        * and columns for initializing the
+                                        * current Trilinos sparsity pattern.
+                                        */
+      template<typename SparsityType>
+      void  
+      reinit (const Epetra_Map   &InputRowMap,
+             const Epetra_Map   &InputColMap,
+             const SparsityType &SP);
+
+                                       /**
+                                        * Reinit function. Takes one of the
+                                        * deal.II sparsity patterns and a
+                                        * parallel partitioning of the rows
+                                        * and columns for initializing the
+                                        * current Trilinos sparsity pattern.
+                                        */
+      template<typename SparsityType>
+      void  
+      reinit (const Epetra_Map   &InputMap,
+             const SparsityType &SP);
+
+                                       /**
+                                        * Copy function. Sets the calling
+                                        * sparsity pattern to be the same as
+                                        * the input sparsity pattern.
+                                        */
+      void  
+      copy_from (const SparsityPattern &SP);
+
+                                       /**
+                                        * Copy function from one of the
+                                        * deal.II sparsity patterns. If used
+                                        * in parallel, this function uses an
+                                        * ad-hoc partitioning of the rows
+                                        * and columns.
+                                        */
+      template<typename SparsityType>
+      void  
+      copy_from (const SparsityType &SP);
+
+                                       /**
+                                        * Release all memory and
+                                        * return to a state just like
+                                        * after having called the
+                                        * default constructor.
+                                       *
+                                       * This is a
+                                       * collective operation that
+                                       * needs to be called on all
+                                       * processors in order to avoid a
+                                       * dead lock.
+                                        */
+      void clear ();
+
+                                       /**
+                                        * In analogy to our own
+                                        * SparsityPattern class, this
+                                        * function compresses the sparsity
+                                        * pattern and allows the resulting
+                                        * pattern to be used for actually
+                                        * generating a matrix. This function
+                                        * also exchanges non-local data that
+                                        * might have accumulated during the
+                                        * addition of new elements. This
+                                        * function must therefore be called
+                                        * once the structure is fixed. This
+                                        * is a collective operation, i.e.,
+                                        * it needs to be run on all
+                                        * processors when used in parallel.
+                                        */
+      void compress ();
+//@}
+/**
+ * @name Information on the sparsity pattern
+ */
+//@{
+
+                                      /**
+                                       * Returns the state of the sparsity
+                                       * pattern, i.e., whether compress()
+                                       * needs to be called after an
+                                       * operation requiring data
+                                       * exchange.
+                                       */
+      bool is_compressed () const;
+
+                                       /**
+                                        * Gives the maximum number of
+                                        * entries per row on the current
+                                        * processor.
+                                        */
+      unsigned int max_entries_per_row () const;
+      
+                                       /**
+                                        * Return the number of rows in this
+                                        * sparsity pattern.
+                                        */
+      unsigned int n_rows () const;
+
+                                       /**
+                                        * Return the number of columns in
+                                        * this sparsity pattern.
+                                        */
+      unsigned int n_cols () const;
+
+                                       /**
+                                        * Return the local dimension of the
+                                        * sparsity pattern, i.e. the number
+                                        * of rows stored on the present MPI
+                                        * process. In the sequential case,
+                                        * this number is the same as
+                                        * n_rows(), but for parallel
+                                        * matrices it may be smaller.
+                                       *
+                                       * To figure out which elements
+                                       * exactly are stored locally,
+                                       * use local_range().
+                                        */
+      unsigned int local_size () const;
+
+                                       /**
+                                       * Return a pair of indices
+                                       * indicating which rows of this
+                                       * sparsity pattern 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,n_rows()), otherwise
+                                       * it will be a pair (i,i+n), where
+                                       * <tt>n=local_size()</tt>.
+                                       */
+      std::pair<unsigned int, unsigned int>
+       local_range () const;
+
+                                      /**
+                                       * Return whether @p index is
+                                       * in the local range or not,
+                                       * see also local_range().
+                                       */
+      bool in_local_range (const unsigned int index) const;
+
+                                       /**
+                                        * Return the number of nonzero
+                                        * elements of this sparsity pattern.
+                                        */
+      unsigned int n_nonzero_elements () const;
+
+                                       /**
+                                        * Number of entries in a
+                                        * specific row.
+                                        */
+      unsigned int row_length (const unsigned int row) const;
+
+                                       /**
+                                        * Return whether the object is
+                                        * empty. It is empty if no memory is
+                                        * allocated, which is the same as
+                                        * when both dimensions are zero.
+                                        */
+      bool empty () const;
+
+                                       /**
+                                        * Return whether the index
+                                        * (<i>i,j</i>) exists in the
+                                        * sparsity pattern (i.e., it may be
+                                        * non-zero) or not.
+                                        */
+      bool exists (const unsigned int i,
+                  const unsigned int j) const;
+
+                                      /**
+                                       * Determine an estimate for the
+                                       * memory consumption (in bytes)
+                                       * of this object. Currently not
+                                       * implemented for this class.
+                                       */
+      unsigned int memory_consumption () const;
+
+//@}
+/**
+ * @name Adding entries
+ */
+//@{
+                                       /**
+                                        * Add the element (<i>i,j</i>) to
+                                        * the sparsity pattern.
+                                        */
+      void add (const unsigned int i,
+                const unsigned int j);
+
+
+                                       /**
+                                        * Add several elements in one row to
+                                        * the sparsity pattern.
+                                        */
+      void add (const unsigned int    row,
+               const unsigned int    n_cols,
+               const unsigned int   *col_indices);
+//@}
+/**
+ * @name Iterators
+ */
+//@{
+
+                                       /**
+                                        * STL-like iterator with the
+                                        * first entry.
+                                        */
+      const_iterator begin () const;
+
+                                       /**
+                                        * Final iterator.
+                                        */
+      const_iterator end () const;
+
+                                       /**
+                                        * STL-like iterator with the
+                                        * first entry of row @p r.
+                                        *
+                                        * Note that if the given row
+                                        * is empty, i.e. does not
+                                        * contain any nonzero entries,
+                                        * then the iterator returned
+                                        * by this function equals
+                                        * <tt>end(r)</tt>. Note also
+                                        * that the iterator may not be
+                                        * dereferencable in that case.
+                                        */
+      const_iterator begin (const unsigned int r) const;
+
+                                       /**
+                                        * Final iterator of row
+                                        * <tt>r</tt>. It points to the
+                                        * first element past the end
+                                        * of line @p r, or past the
+                                        * end of the entire sparsity
+                                        * pattern.
+                                        *
+                                        * Note that the end iterator
+                                        * is not necessarily
+                                        * dereferencable. This is in
+                                        * particular the case if it is
+                                        * the end iterator for the
+                                        * last row of a matrix.
+                                        */
+      const_iterator end (const unsigned int r) const;
+
+//@}
+/**
+ * @name Input/Output
+ */
+//@{
+
+                                       /**
+                                       * 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 sparsity pattern to the
+                                       * given stream, using the format
+                                       * <tt>(line,col)</tt>.
+                                       */
+      void print (std::ostream &out) const;
+    
+                                        // TODO: Write an overloading
+                                        // of the operator << for output.
+                                        // Since the underlying Trilinos 
+                                        // object supports it, this should 
+                                        // be very easy.
+
+//@}
+                                    /** @addtogroup Exceptions
+                                     * @{ */
+                                       /**
+                                        * Exception
+                                        */
+      DeclException1 (ExcTrilinosError,
+                      int,
+                      << "An error with error number " << arg1
+                      << " occured while calling a Trilinos function");
+
+                                      /**
+                                       * Exception
+                                       */
+      DeclException2 (ExcInvalidIndex,
+                     int, int,
+                     << "The entry with index <" << arg1 << ',' << arg2
+                     << "> does not exist.");
+
+                                       /**
+                                        * Exception
+                                        */
+      DeclException0 (ExcSourceEqualsDestination);
+
+                                       /**
+                                        * Exception
+                                        */
+      DeclException4 (ExcAccessToNonLocalElement,
+                     int, int, int, int,
+                     << "You tried to access element (" << arg1
+                     << "/" << arg2 << ")"
+                     << " of a distributed matrix, but only rows "
+                     << arg3 << " through " << arg4
+                     << " are stored locally and can be accessed.");
+
+                                       /**
+                                        * Exception
+                                        */
+      DeclException2 (ExcAccessToNonPresentElement,
+                     int, int,
+                     << "You tried to access element (" << arg1
+                     << "/" << arg2 << ")"
+                     << " of a sparse matrix, but it appears to not"
+                     << " exist in the Trilinos sparsity pattern.");
+                                    //@}    
+    private:
+                                       /**
+                                       * 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 row_map;
+
+                                       /**
+                                        * 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;
+
+                                      /**
+                                       * A boolean variable to hold
+                                       * information on whether the
+                                       * vector is compressed or not.
+                                       */
+      bool compressed;
+
+                                       /**
+                                        * A sparsity pattern object in
+                                        * Trilinos to be used for finite
+                                        * element based problems which
+                                        * allows for adding non-local
+                                        * elements to the pattern.
+                                        */
+      std::auto_ptr<Epetra_FECrsGraph> graph;
+
+      friend class SparseMatrix;
+      friend class SparsityPatternIterators::const_iterator;
+  };
+
+
+
+// -------------------------- inline and template functions ----------------------
+
+
+#ifndef DOXYGEN
+
+  namespace SparsityPatternIterators
+  {
+
+    inline
+    const_iterator::Accessor::
+    Accessor (const SparsityPattern *sp,
+              const unsigned int     row,
+              const unsigned int     index)
+                    :
+                    sparsity_pattern(const_cast<SparsityPattern*>(sp)),
+                    a_row(row),
+                    a_index(index)
+    {
+      visit_present_row ();
+    }
+
+
+    inline
+    unsigned int
+    const_iterator::Accessor::row() const
+    {
+      Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern());
+      return a_row;
+    }
+
+
+
+    inline
+    unsigned int
+    const_iterator::Accessor::column() const
+    {
+      Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern());
+      return (*colnum_cache)[a_index];
+    }
+
+
+
+    inline
+    unsigned int
+    const_iterator::Accessor::index() const
+    {
+      Assert (a_row < sparsity_pattern->n_rows(), ExcBeyondEndOfSparsityPattern());
+      return a_index;
+    }
+
+
+
+    inline
+    const_iterator::
+    const_iterator(const SparsityPattern *sp,
+                   const unsigned int     row,
+                   const unsigned int     index)
+                    :
+                    accessor(sp, row, index)
+    {}
+
+
+
+    inline
+    const_iterator &
+    const_iterator::operator++ ()
+    {
+      Assert (accessor.a_row < accessor.sparsity_pattern->n_rows(), 
+             ExcIteratorPastEnd());
+
+      ++accessor.a_index;
+
+                                       // If at end of line: do one
+                                       // step, then cycle until we 
+                                       // find a row with a nonzero
+                                       // number of entries.
+      if (accessor.a_index >= accessor.colnum_cache->size())
+        {
+          accessor.a_index = 0;
+          ++accessor.a_row;
+      
+          while ((accessor.a_row < accessor.sparsity_pattern->n_rows())
+                 &&
+                 (accessor.sparsity_pattern->row_length(accessor.a_row) == 0))
+            ++accessor.a_row;
+
+          accessor.visit_present_row();
+        }
+      return *this;
+    }
+
+
+
+    inline
+    const_iterator
+    const_iterator::operator++ (int)
+    {
+      const const_iterator old_state = *this;
+      ++(*this);
+      return old_state;
+    }
+
+
+
+    inline
+    const const_iterator::Accessor &
+    const_iterator::operator* () const
+    {
+      return accessor;
+    }
+
+
+
+    inline
+    const const_iterator::Accessor *
+    const_iterator::operator-> () const
+    {
+      return &accessor;
+    }
+
+
+
+    inline
+    bool
+    const_iterator::
+    operator == (const const_iterator& other) const
+    {
+      return (accessor.a_row == other.accessor.a_row &&
+              accessor.a_index == other.accessor.a_index);
+    }
+
+
+
+    inline
+    bool
+    const_iterator::
+    operator != (const const_iterator& other) const
+    {
+      return ! (*this == other);
+    }
+
+
+
+    inline
+    bool
+    const_iterator::
+    operator < (const const_iterator& other) const
+    {
+      return (accessor.row() < other.accessor.row() ||
+              (accessor.row() == other.accessor.row() &&
+               accessor.index() < other.accessor.index()));
+    }
+    
+  }
+
+
+  inline
+  SparsityPattern::const_iterator
+  SparsityPattern::begin() const
+  {
+    return const_iterator(this, 0, 0);
+  }
+
+
+
+  inline
+  SparsityPattern::const_iterator
+  SparsityPattern::end() const
+  {
+    return const_iterator(this, n_rows(), 0);
+  }
+
+
+
+  inline
+  SparsityPattern::const_iterator
+  SparsityPattern::begin(const unsigned int r) const
+  {
+    Assert (r < n_rows(), ExcIndexRange(r, 0, n_rows()));
+    if (row_length(r) > 0)
+      return const_iterator(this, r, 0);
+    else
+      return end (r);
+  }
+
+
+
+  inline
+  SparsityPattern::const_iterator
+  SparsityPattern::end(const unsigned int r) const
+  {
+    Assert (r < n_rows(), ExcIndexRange(r, 0, n_rows()));
+
+                                     // place the iterator on the first entry
+                                     // past this line, or at the end of the
+                                     // matrix
+    for (unsigned int i=r+1; i<n_rows(); ++i)
+      if (row_length(i) > 0)
+        return const_iterator(this, i, 0);
+    
+                                     // if there is no such line, then take the
+                                     // end iterator of the matrix
+    return end();
+  }
+
+
+
+  inline
+  bool
+  SparsityPattern::in_local_range (const unsigned int index) const
+  {
+    int begin, end;
+    begin = graph->RowMap().MinMyGID();
+    end = graph->RowMap().MaxMyGID()+1;
+    
+    return ((index >= static_cast<unsigned int>(begin)) &&
+            (index < static_cast<unsigned int>(end)));
+  }
+
+
+
+  inline
+  bool
+  SparsityPattern::is_compressed () const
+  {
+    return compressed;
+  }
+
+
+
+  inline
+  bool
+  SparsityPattern::empty () const
+  {
+    return ((n_rows() == 0) && (n_cols() == 0));
+  }
+
+
+
+  inline
+  void
+  SparsityPattern::add (const unsigned int   i,
+                       const unsigned int   j)
+  {
+    add (i, 1, &j);
+  }
+
+
+
+  inline
+  void
+  SparsityPattern::add (const unsigned int    row,
+                       const unsigned int    n_cols,
+                       const unsigned int   *col_indices)
+  {
+    int * col_index_ptr = (int*)col_indices;
+    compressed = false;
+
+    int ierr;
+
+                                  // If the calling sparsity pattern owns
+                                  // the row to which we want to add
+                                  // values, we can directly call the
+                                  // Epetra_CrsGraph input function, which
+                                  // is much faster than the
+                                  // Epetra_FECrsGraph function.
+    if (row_map.MyGID(row) == true)
+      ierr = graph->Epetra_CrsGraph::InsertGlobalIndices(row, 
+                                                        n_cols,
+                                                        col_index_ptr);
+    else
+      {
+                                  // When we're at off-processor data, we
+                                  // have to stick with the standard
+                                  // SumIntoGlobalValues
+                                  // function. Nevertheless, the way we
+                                  // call it is the fastest one (any other
+                                  // will lead to repeated allocation and
+                                  // deallocation of memory in order to
+                                  // call the function we already use,
+                                  // which is very unefficient if writing
+                                  // one element at a time).
+
+       ierr = graph->InsertGlobalIndices (1, (int*)&row, n_cols, 
+                                          col_index_ptr);
+      }
+
+    //Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
+    AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
+  }
+
+
+#endif // DOXYGEN      
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif // DEAL_II_USE_TRILINOS
+
+
+/*--------------------   trilinos_sparsity_pattern.h     --------------------*/
+
+#endif
+/*--------------------   trilinos_sparsity_pattern.h     --------------------*/
index 6f6ea0be85a2b73fa6f3d88cab19d753d532f718..37481370fc7828aea57215020c870ec66737ac93 100644 (file)
@@ -309,7 +309,9 @@ template class BlockSparsityPatternBase<SparsityPattern>;
 template class BlockSparsityPatternBase<CompressedSparsityPattern>;
 template class BlockSparsityPatternBase<CompressedSimpleSparsityPattern>;
 template class BlockSparsityPatternBase<CompressedSetSparsityPattern>;
-
+#ifdef DEAL_II_USE_TRILINOS
+template class BlockSparsityPatternBase<TrilinosWrappers::SparsityPattern>;
+#endif
 
 
 
@@ -473,25 +475,27 @@ BlockCompressedSparsityPattern::reinit (
   this->collect_sizes();  
 }
 
-BlockCompressedSimpleSparsityPattern::BlockCompressedSimpleSparsityPattern ()
+
+
+BlockCompressedSetSparsityPattern::BlockCompressedSetSparsityPattern ()
 {}
 
 
 
-BlockCompressedSimpleSparsityPattern::
-BlockCompressedSimpleSparsityPattern (const unsigned int n_rows,
+BlockCompressedSetSparsityPattern::
+BlockCompressedSetSparsityPattern (const unsigned int n_rows,
                                const unsigned int n_columns)
                :
-               BlockSparsityPatternBase<CompressedSimpleSparsityPattern>(n_rows,
+               BlockSparsityPatternBase<CompressedSetSparsityPattern>(n_rows,
                                                                    n_columns)
 {}
 
 
-BlockCompressedSimpleSparsityPattern::
-BlockCompressedSimpleSparsityPattern (const std::vector<unsigned int>& row_indices,
+BlockCompressedSetSparsityPattern::
+BlockCompressedSetSparsityPattern (const std::vector<unsigned int>& row_indices,
                                const std::vector<unsigned int>& col_indices)
                :
-               BlockSparsityPatternBase<CompressedSimpleSparsityPattern>(row_indices.size(),
+               BlockSparsityPatternBase<CompressedSetSparsityPattern>(row_indices.size(),
                                                                    col_indices.size())
 {
   for (unsigned int i=0;i<row_indices.size();++i)
@@ -502,11 +506,11 @@ BlockCompressedSimpleSparsityPattern (const std::vector<unsigned int>& row_indic
 
 
 void
-BlockCompressedSimpleSparsityPattern::reinit (
+BlockCompressedSetSparsityPattern::reinit (
   const std::vector< unsigned int > &row_block_sizes,
   const std::vector< unsigned int > &col_block_sizes)
 {
-  BlockSparsityPatternBase<CompressedSimpleSparsityPattern>::reinit(row_block_sizes.size(), col_block_sizes.size());
+  BlockSparsityPatternBase<CompressedSetSparsityPattern>::reinit(row_block_sizes.size(), col_block_sizes.size());
   for (unsigned int i=0;i<row_block_sizes.size();++i)
     for (unsigned int j=0;j<col_block_sizes.size();++j)
       this->block(i,j).reinit(row_block_sizes[i],col_block_sizes[j]);
@@ -515,25 +519,26 @@ BlockCompressedSimpleSparsityPattern::reinit (
 
 
 
-BlockCompressedSetSparsityPattern::BlockCompressedSetSparsityPattern ()
+BlockCompressedSimpleSparsityPattern::BlockCompressedSimpleSparsityPattern ()
 {}
 
 
 
-BlockCompressedSetSparsityPattern::
-BlockCompressedSetSparsityPattern (const unsigned int n_rows,
+BlockCompressedSimpleSparsityPattern::
+BlockCompressedSimpleSparsityPattern (const unsigned int n_rows,
                                const unsigned int n_columns)
                :
-               BlockSparsityPatternBase<CompressedSetSparsityPattern>(n_rows,
+               BlockSparsityPatternBase<CompressedSimpleSparsityPattern>(n_rows,
                                                                    n_columns)
 {}
 
 
-BlockCompressedSetSparsityPattern::
-BlockCompressedSetSparsityPattern (const std::vector<unsigned int>& row_indices,
+
+BlockCompressedSimpleSparsityPattern::
+BlockCompressedSimpleSparsityPattern (const std::vector<unsigned int>& row_indices,
                                const std::vector<unsigned int>& col_indices)
                :
-               BlockSparsityPatternBase<CompressedSetSparsityPattern>(row_indices.size(),
+               BlockSparsityPatternBase<CompressedSimpleSparsityPattern>(row_indices.size(),
                                                                    col_indices.size())
 {
   for (unsigned int i=0;i<row_indices.size();++i)
@@ -543,12 +548,13 @@ BlockCompressedSetSparsityPattern (const std::vector<unsigned int>& row_indices,
 }
 
 
+
 void
-BlockCompressedSetSparsityPattern::reinit (
+BlockCompressedSimpleSparsityPattern::reinit (
   const std::vector< unsigned int > &row_block_sizes,
   const std::vector< unsigned int > &col_block_sizes)
 {
-  BlockSparsityPatternBase<CompressedSetSparsityPattern>::reinit(row_block_sizes.size(), col_block_sizes.size());
+  BlockSparsityPatternBase<CompressedSimpleSparsityPattern>::reinit(row_block_sizes.size(), col_block_sizes.size());
   for (unsigned int i=0;i<row_block_sizes.size();++i)
     for (unsigned int j=0;j<col_block_sizes.size();++j)
       this->block(i,j).reinit(row_block_sizes[i],col_block_sizes[j]);
@@ -557,6 +563,81 @@ BlockCompressedSetSparsityPattern::reinit (
 
 
 
+#ifdef DEAL_II_USE_TRILINOS
+namespace TrilinosWrappers
+{
+
+  BlockSparsityPattern::BlockSparsityPattern ()
+  {}
+
+
+
+  BlockSparsityPattern::
+  BlockSparsityPattern (const unsigned int n_rows,
+                       const unsigned int n_columns)
+                :
+                dealii::BlockSparsityPatternBase<SparsityPattern>(n_rows,
+                                                                 n_columns)
+  {}
+
+
+
+  BlockSparsityPattern::
+  BlockSparsityPattern (const std::vector<unsigned int>& row_indices,
+                       const std::vector<unsigned int>& col_indices)
+               :
+               BlockSparsityPatternBase<SparsityPattern>(row_indices.size(),
+                                                         col_indices.size())
+  {
+    for (unsigned int i=0;i<row_indices.size();++i)
+      for (unsigned int j=0;j<col_indices.size();++j)
+       this->block(i,j).reinit(row_indices[i],col_indices[j]);
+    this->collect_sizes();
+  }
+
+
+
+  BlockSparsityPattern::
+  BlockSparsityPattern (const std::vector<Epetra_Map>& input_maps)
+               :
+               BlockSparsityPatternBase<SparsityPattern>(input_maps.size(),
+                                                         input_maps.size())
+  {
+    for (unsigned int i=0;i<input_maps.size();++i)
+      for (unsigned int j=0;j<input_maps.size();++j)
+       this->block(i,j).reinit(input_maps[i],input_maps[j]);
+    this->collect_sizes();
+  }
+
+
+
+  void
+  BlockSparsityPattern::reinit (const std::vector<unsigned int> &row_block_sizes,
+                               const std::vector<unsigned int> &col_block_sizes)
+  {
+    dealii::BlockSparsityPatternBase<SparsityPattern>::reinit(row_block_sizes.size(), col_block_sizes.size());
+    for (unsigned int i=0;i<row_block_sizes.size();++i)
+      for (unsigned int j=0;j<col_block_sizes.size();++j)
+       this->block(i,j).reinit(row_block_sizes[i],col_block_sizes[j]);
+    this->collect_sizes();  
+  }
+
+
+
+  void
+  BlockSparsityPattern::reinit (const std::vector<Epetra_Map> &input_maps)
+  {
+    dealii::BlockSparsityPatternBase<SparsityPattern>::reinit(input_maps.size(), 
+                                                             input_maps.size());
+    for (unsigned int i=0;i<input_maps.size();++i)
+      for (unsigned int j=0;j<input_maps.size();++j)
+       this->block(i,j).reinit(input_maps[i],input_maps[j]);
+    this->collect_sizes();  
+  }
+
+}
+
+#endif
 
 // Remark: The explicit instantiations for "BlockSparsityPatternBase" were moved
 // to the top of this source file. The reason is a slightly buggy version
index 8fc610188087486b351473788b92b61557495b23..eea3d942dcfd30b7c8f84aec8c8833200ca8478b 100644 (file)
@@ -164,6 +164,31 @@ namespace TrilinosWrappers
 
 
 
+  template <>  
+  void
+  BlockSparseMatrix::
+  reinit (const BlockSparsityPattern    &block_sparsity_pattern)
+  {
+  
+                                    // Call the other basic reinit function, ...
+    reinit (block_sparsity_pattern.n_block_rows(),
+           block_sparsity_pattern.n_block_cols());
+
+                                    // ... set the correct sizes, ...
+    this->row_block_indices    = block_sparsity_pattern.get_row_indices();
+    this->column_block_indices = block_sparsity_pattern.get_column_indices();
+       
+                                    // ... and then assign the correct
+                                    // data to the blocks.
+    for (unsigned int r=0; r<this->n_block_rows(); ++r)
+      for (unsigned int c=0; c<this->n_block_cols(); ++c)
+        {
+         this->sub_objects[r][c]->reinit (block_sparsity_pattern.block(r,c));
+        }
+  }
+
+
+
   void
   BlockSparseMatrix::
   reinit (const std::vector<Epetra_Map>             &input_maps,
@@ -369,7 +394,7 @@ namespace TrilinosWrappers
   // -------------------- explicit instantiations -----------------------
   //
   template void
-  BlockSparseMatrix::reinit (const BlockSparsityPattern &);
+  BlockSparseMatrix::reinit (const dealii::BlockSparsityPattern &);
   template void
   BlockSparseMatrix::reinit (const BlockCompressedSparsityPattern &);
   template void
@@ -380,7 +405,7 @@ namespace TrilinosWrappers
 
   template void
   BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
-                            const BlockSparsityPattern    &);
+                            const dealii::BlockSparsityPattern    &);
   template void
   BlockSparseMatrix::reinit (const std::vector<Epetra_Map> &,
                             const BlockCompressedSparsityPattern &);
index 5021328f0d93527b5f07a517051477ecb80aca6a..4b5fa34598d3c21ebb8b594b46be08c3a932afae 100644 (file)
@@ -276,7 +276,7 @@ namespace TrilinosWrappers
          break;
         case gmres:
          aztecBlockParams->sublist("Forward Solve")
-           .sublist("AztecOO Settings").set("Aztec Solver", "GMRES");
+           .sublist("AztecOO Settings").set("Aztec Solver", "MinRes");
          aztecBlockParams->sublist("Forward Solve")
            .sublist("AztecOO Settings").set("Size of Krylov Subspace", 
                                             (int)additional_data.gmres_restart_parameter);
index 3d4f2dedaf7912b3615ad6db62f9aeb1d4fea1be..31ced248c3d28329cf14b9406c5c428a325c5673 100755 (executable)
@@ -13,6 +13,7 @@
 
 #include <lac/trilinos_sparse_matrix.h>
 
+#include <lac/trilinos_sparsity_pattern.h>
 #include <lac/sparsity_pattern.h>
 #include <lac/compressed_sparsity_pattern.h>
 #include <lac/compressed_set_sparsity_pattern.h>
@@ -97,6 +98,8 @@ namespace TrilinosWrappers
     matrix->FillComplete();
   }
 
+
+
   SparseMatrix::SparseMatrix (const Epetra_Map  &InputMap,
                              const unsigned int n_max_entries_per_row)
                  :
@@ -109,6 +112,8 @@ namespace TrilinosWrappers
                                        int(n_max_entries_per_row), false)))
   {}
 
+
+
   SparseMatrix::SparseMatrix (const Epetra_Map                &InputMap,
                              const std::vector<unsigned int> &n_entries_per_row)
                  :
@@ -122,6 +127,8 @@ namespace TrilinosWrappers
                                            false)))
   {}
 
+
+
   SparseMatrix::SparseMatrix (const Epetra_Map  &InputRowMap,
                              const Epetra_Map  &InputColMap,
                              const unsigned int n_max_entries_per_row)
@@ -135,6 +142,8 @@ namespace TrilinosWrappers
                                        int(n_max_entries_per_row), false)))
   {}
 
+
+
   SparseMatrix::SparseMatrix (const Epetra_Map                &InputRowMap,
                              const Epetra_Map                &InputColMap,
                              const std::vector<unsigned int> &n_entries_per_row)
@@ -149,6 +158,8 @@ namespace TrilinosWrappers
                                            false)))
   {}
 
+
+
   SparseMatrix::SparseMatrix (const unsigned int m,
                              const unsigned int n,
                              const unsigned int n_max_entries_per_row)
@@ -167,6 +178,8 @@ namespace TrilinosWrappers
                                        int(n_max_entries_per_row), false)))
   {}
 
+
+
   SparseMatrix::SparseMatrix (const unsigned int               m,
                              const unsigned int               n,
                              const std::vector<unsigned int> &n_entries_per_row)
@@ -186,6 +199,25 @@ namespace TrilinosWrappers
                                             false)))
   {}
 
+
+
+  SparseMatrix::SparseMatrix (const SparsityPattern &InputSP)
+                 :
+                  Subscriptor(),
+                  row_map (InputSP.row_map),
+                 col_map (InputSP.col_map),
+                 last_action (Zero),
+                 compressed (true),
+                 matrix (std::auto_ptr<Epetra_FECrsMatrix>
+                         (new Epetra_FECrsMatrix(Copy, *InputSP.graph, false)))
+  {
+    Assert(InputSP.graph->Filled() == true,
+          ExcMessage("The Trilinos sparsity pattern has not been compressed."));
+    compress();
+  }
+
+
+
   SparseMatrix::SparseMatrix (const SparseMatrix &InputMatrix)
                  :
                   Subscriptor(),
@@ -474,6 +506,25 @@ namespace TrilinosWrappers
 
 
 
+  void
+  SparseMatrix::reinit (const SparsityPattern &sparsity_pattern)
+  {
+    matrix.reset();
+
+    row_map = sparsity_pattern.row_map;
+    col_map = sparsity_pattern.col_map;
+
+    Assert (sparsity_pattern.graph->Filled() == true,
+           ExcMessage("The Trilinos sparsity pattern has not been compressed"));
+
+    matrix = std::auto_ptr<Epetra_FECrsMatrix>
+      (new Epetra_FECrsMatrix(Copy, *sparsity_pattern.graph, false));
+
+    compress();
+  }
+
+
+
   void
   SparseMatrix::reinit (const SparseMatrix &sparse_matrix)
   {
@@ -748,20 +799,20 @@ namespace TrilinosWrappers
 
                                      // If the data is not on the
                                      // present processor, we throw
-                                     // an exception. This is on of
+                                     // an exception. This is one of
                                      // the two tiny differences to
                                      // the el(i,j) call, which does
                                      // not throw any assertions.
-    if ((trilinos_i == -1 ) || (trilinos_j == -1))
+    if (trilinos_i == -1)
       {
        Assert (false, ExcAccessToNonLocalElement(i, j, local_range().first,
                                                  local_range().second));
       }
     else
       {
-                                     // Check whether the matrix 
-                                     // already is transformed to
-                                     // local indices.
+                                     // Check whether the matrix has
+                                     // already been transformed to local
+                                     // indices.
        if (matrix->Filled() == false)
          matrix->GlobalAssemble(col_map, row_map, true);
 
@@ -914,7 +965,7 @@ namespace TrilinosWrappers
   unsigned int
   SparseMatrix::n () const
   {
-    int n_cols = matrix -> NumGlobalCols();
+    unsigned int n_cols = matrix -> NumGlobalCols();
     return n_cols;
   }
 
@@ -923,7 +974,7 @@ namespace TrilinosWrappers
   unsigned int
   SparseMatrix::local_size () const
   {
-    int n_rows = matrix -> NumMyRows();
+    unsigned int n_rows = matrix -> NumMyRows();
 
     return n_rows;
   }
@@ -933,7 +984,7 @@ namespace TrilinosWrappers
   std::pair<unsigned int, unsigned int>
   SparseMatrix::local_range () const
   {
-    int begin, end;
+    unsigned int begin, end;
     begin = matrix -> RowMap().MinMyGID();
     end = matrix -> RowMap().MaxMyGID()+1;
     
@@ -945,9 +996,9 @@ namespace TrilinosWrappers
   unsigned int
   SparseMatrix::n_nonzero_elements () const
   {
-    int nnz = matrix->NumGlobalNonzeros();
+    unsigned int nnz = matrix->NumGlobalNonzeros();
 
-    return static_cast<unsigned int>(nnz);
+    return nnz;
   }
 
 
@@ -960,7 +1011,7 @@ namespace TrilinosWrappers
                                  // get a representation of the
                                  // present row
     int ncols = -1;
-    int local_row = matrix->RowMap().LID(row);
+    int local_row = matrix->LRID(row);
 
                                  // on the processor who owns this
                                  // row, we'll have a non-negative
@@ -975,6 +1026,7 @@ namespace TrilinosWrappers
   }
 
 
+
   TrilinosScalar
   SparseMatrix::l1_norm () const
   {
@@ -1374,7 +1426,7 @@ namespace TrilinosWrappers
   // explicit instantiations
   //
   template void
-  SparseMatrix::reinit (const SparsityPattern &);
+  SparseMatrix::reinit (const dealii::SparsityPattern &);
   template void
   SparseMatrix::reinit (const CompressedSparsityPattern &);
   template void
@@ -1385,7 +1437,7 @@ namespace TrilinosWrappers
 
   template void
   SparseMatrix::reinit (const Epetra_Map &,
-                       const SparsityPattern &);
+                       const dealii::SparsityPattern &);
   template void
   SparseMatrix::reinit (const Epetra_Map &,
                        const CompressedSparsityPattern &);
@@ -1400,7 +1452,7 @@ namespace TrilinosWrappers
   template void
   SparseMatrix::reinit (const Epetra_Map &,
                        const Epetra_Map &,
-                       const SparsityPattern &);
+                       const dealii::SparsityPattern &);
   template void
   SparseMatrix::reinit (const Epetra_Map &,
                        const Epetra_Map &,
diff --git a/deal.II/lac/source/trilinos_sparsity_pattern.cc b/deal.II/lac/source/trilinos_sparsity_pattern.cc
new file mode 100755 (executable)
index 0000000..5de4418
--- /dev/null
@@ -0,0 +1,735 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2008 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//---------------------------------------------------------------------------
+
+#include <lac/trilinos_sparsity_pattern.h>
+
+#include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/compressed_set_sparsity_pattern.h>
+#include <lac/compressed_simple_sparsity_pattern.h>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+  namespace SparsityPatternIterators
+  {
+    void
+    SparsityPattern::const_iterator::Accessor::
+    visit_present_row ()
+    {
+                                 // if we are asked to visit the
+                                 // past-the-end line, then simply
+                                 // release all our caches and go on
+                                 // with life
+      if (this->a_row == sparsity_pattern->n_rows())
+       {
+         colnum_cache.reset ();
+
+         return;
+       }
+      
+                                 // otherwise first flush Trilinos caches
+      sparsity_pattern->compress ();
+
+                                 // get a representation of the present
+                                 // row
+      int ncols;
+      int colnums = sparsity_pattern->n_cols();
+      
+      int ierr;
+      ierr = sparsity_pattern->graph->ExtractGlobalRowCopy((int)this->a_row, 
+                                                          colnums,
+                                                          ncols, 
+                                                          (int*)&(*colnum_cache)[0]);
+      AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+                                 // copy it into our caches if the
+                                 // line isn't empty. if it is, then
+                                 // we've done something wrong, since
+                                 // we shouldn't have initialized an
+                                 // iterator for an empty line (what
+                                 // would it point to?)
+      Assert (ncols != 0, ExcInternalError());
+      colnum_cache.reset (new std::vector<unsigned int> (colnums,
+                                                        colnums+ncols));
+    }
+  }
+
+
+                                 // The constructor is actually the
+                                 // only point where we have to check
+                                 // whether we build a serial or a
+                                 // parallel Trilinos matrix.
+                                 // Actually, it does not even matter
+                                 // how many threads there are, but
+                                 // only if we use an MPI compiler or
+                                 // a standard compiler. So, even one
+                                 // thread on a configuration with
+                                 // MPI will still get a parallel
+                                 // interface.
+  SparsityPattern::SparsityPattern ()
+                 :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+                  row_map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+#else
+                  row_map (0, 0, Epetra_SerialComm()),
+#endif
+                 col_map (row_map),
+                 compressed (true),
+                 graph (std::auto_ptr<Epetra_FECrsGraph>
+                                 (new Epetra_FECrsGraph(View, row_map, 0)))
+  {
+    graph->FillComplete();
+  }
+
+  SparsityPattern::SparsityPattern (const Epetra_Map  &InputMap,
+                                   const unsigned int n_entries_per_row)
+                 :
+                  row_map (InputMap),
+                 col_map (row_map),
+                 compressed (false),
+                 graph (std::auto_ptr<Epetra_FECrsGraph>
+                                 (new Epetra_FECrsGraph(Copy, row_map, 
+                                       int(n_entries_per_row), false)))
+  {}
+
+  SparsityPattern::SparsityPattern (const Epetra_Map                &InputMap,
+                                   const std::vector<unsigned int> &n_entries_per_row)
+                 :
+                  row_map (InputMap),
+                 col_map (row_map),
+                 compressed (false),
+                 graph (std::auto_ptr<Epetra_FECrsGraph>
+                   (new Epetra_FECrsGraph(Copy, row_map, 
+                     (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
+                                           false)))
+  {}
+
+  SparsityPattern::SparsityPattern (const Epetra_Map  &InputRowMap,
+                                   const Epetra_Map  &InputColMap,
+                                   const unsigned int n_entries_per_row)
+                 :
+                  row_map (InputRowMap),
+                  col_map (InputColMap),
+                 compressed (false),
+                 graph (std::auto_ptr<Epetra_FECrsGraph>
+                                 (new Epetra_FECrsGraph(Copy, row_map, 
+                                       int(n_entries_per_row), false)))
+  {}
+
+  SparsityPattern::SparsityPattern (const Epetra_Map                &InputRowMap,
+                                   const Epetra_Map                &InputColMap,
+                                   const std::vector<unsigned int> &n_entries_per_row)
+                 :
+                  row_map (InputRowMap),
+                  col_map (InputColMap),
+                 compressed (false),
+                 graph (std::auto_ptr<Epetra_FECrsGraph>
+                   (new Epetra_FECrsGraph(Copy, row_map, 
+                     (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])),
+                                           false)))
+  {}
+
+  SparsityPattern::SparsityPattern (const unsigned int m,
+                                   const unsigned int n,
+                                   const unsigned int n_entries_per_row)
+                 :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+                  row_map (m, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+                  col_map (n, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+#else
+                  row_map (m, 0, Epetra_SerialComm()),
+                  col_map (n, 0, Epetra_SerialComm()),
+#endif
+                 compressed (false),
+                 graph (std::auto_ptr<Epetra_FECrsGraph>
+                               (new Epetra_FECrsGraph(Copy, row_map, 
+                                       int(n_entries_per_row), false)))
+  {}
+
+  SparsityPattern::SparsityPattern (const unsigned int               m,
+                                   const unsigned int               n,
+                                   const std::vector<unsigned int> &n_entries_per_row)
+                 :
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+                  row_map (m, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+                  col_map (n, 0, Epetra_MpiComm(MPI_COMM_WORLD)),
+#else
+                  row_map (m, 0, Epetra_SerialComm()),
+                  col_map (n, 0, Epetra_SerialComm()),
+#endif
+                 compressed (false),
+                 graph (std::auto_ptr<Epetra_FECrsGraph>
+                    (new Epetra_FECrsGraph(Copy, row_map, 
+                       (int*)const_cast<unsigned int*>(&(n_entries_per_row[0])), 
+                                            false)))
+  {}
+
+                                  // Copy function is currently not working
+                                  // because the Trilinos Epetra_FECrsGraph
+                                  // does not implement a reinit function
+                                  // from another graph.
+  /*
+  SparsityPattern::SparsityPattern (const SparsityPattern &InputSP)
+                 :
+                  Subscriptor(),
+                 row_map (InputSP.row_map),
+                 col_map (InputSP.col_map),
+                 compressed (false),
+                 graph (std::auto_ptr<Epetra_FECrsGraph>
+                         (new Epetra_FECrsGraph(*InputSP.graph)))
+  {}
+  */
+
+
+
+  SparsityPattern::~SparsityPattern ()
+  {}
+
+
+
+  void 
+  SparsityPattern::reinit (const Epetra_Map   &input_map,
+                          const unsigned int  n_entries_per_row)
+  {
+    reinit (input_map, input_map, n_entries_per_row);
+  }
+
+
+  void 
+  SparsityPattern::reinit (const unsigned int  m,
+                          const unsigned int  n,
+                          const unsigned int  n_entries_per_row)
+  {
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+    Epetra_MpiComm    trilinos_communicator (MPI_COMM_WORLD);
+#else
+    Epetra_SerialComm trilinos_communicator;
+#endif
+
+    const Epetra_Map rows (m, 0, trilinos_communicator);
+    const Epetra_Map columns (n, 0, trilinos_communicator);
+
+    reinit (rows, columns, n_entries_per_row);
+  }
+
+
+  void 
+  SparsityPattern::reinit (const Epetra_Map   &input_row_map,
+                          const Epetra_Map   &input_col_map,
+                          const unsigned int  n_entries_per_row)
+  {
+    graph.reset();
+
+    row_map = input_row_map;
+    col_map = input_col_map;
+
+    graph = std::auto_ptr<Epetra_FECrsGraph>
+      (new Epetra_FECrsGraph(Copy, row_map, n_entries_per_row, false));
+  }
+
+
+
+  void 
+  SparsityPattern::reinit (const Epetra_Map   &input_map,
+                          const std::vector<unsigned int> &n_entries_per_row)
+  {
+    reinit (input_map, input_map, n_entries_per_row);
+  }
+
+
+
+  void 
+  SparsityPattern::reinit (const unsigned int  m,
+                          const unsigned int  n,
+                          const std::vector<unsigned int> &n_entries_per_row)
+  {
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+    Epetra_MpiComm    trilinos_communicator (MPI_COMM_WORLD);
+#else
+    Epetra_SerialComm trilinos_communicator;
+#endif
+
+    const Epetra_Map rows (m, 0, trilinos_communicator);
+    const Epetra_Map columns (n, 0, trilinos_communicator);
+
+    reinit (rows, columns, n_entries_per_row);
+  }
+
+
+
+  void 
+  SparsityPattern::reinit (const Epetra_Map   &input_row_map,
+                          const Epetra_Map   &input_col_map,
+                          const std::vector<unsigned int> &n_entries_per_row)
+  {
+    graph.reset();
+
+    Assert (n_entries_per_row.size() == 
+             static_cast<unsigned int>(input_row_map.NumGlobalElements()),
+           ExcDimensionMismatch (n_entries_per_row.size(),
+                                 input_row_map.NumGlobalElements()));
+    row_map = input_row_map;
+    col_map = input_col_map;
+
+    graph = std::auto_ptr<Epetra_FECrsGraph>
+      (new Epetra_FECrsGraph(Copy, row_map, 
+                            n_entries_per_row[input_row_map.MinMyGID()], 
+                            false));
+  }
+
+
+
+  template <typename SparsityType>
+  void 
+  SparsityPattern::reinit (const Epetra_Map   &input_map,
+                          const SparsityType &sp)
+  {
+    reinit (input_map, input_map, sp);
+  }
+
+
+
+  template <typename SparsityType>
+  void 
+  SparsityPattern::reinit (const Epetra_Map   &input_row_map,
+                          const Epetra_Map   &input_col_map,
+                          const SparsityType &sp)
+  {
+    graph.reset();
+
+    Assert (sp.n_rows() == 
+             static_cast<unsigned int>(input_row_map.NumGlobalElements()),
+           ExcDimensionMismatch (sp.n_rows(),
+                                 input_row_map.NumGlobalElements()));
+    Assert (sp.n_cols() == 
+             static_cast<unsigned int>(input_col_map.NumGlobalElements()),
+           ExcDimensionMismatch (sp.n_cols(),
+                                 input_col_map.NumGlobalElements()));
+
+    row_map = input_row_map;
+    col_map = input_col_map;
+
+    const unsigned int n_rows = sp.n_rows();
+
+    std::vector<int> n_entries_per_row(n_rows);
+
+    for (unsigned int row=0; row<n_rows; ++row)
+      n_entries_per_row[row] = sp.row_length(row);
+
+    graph = std::auto_ptr<Epetra_FECrsGraph>
+      (new Epetra_FECrsGraph(Copy, row_map, 
+                            n_entries_per_row[input_row_map.MinMyGID()], 
+                            false));
+
+    Assert (graph->NumGlobalRows() == (int)sp.n_rows(),
+           ExcDimensionMismatch (graph->NumGlobalRows(),
+                                 sp.n_rows()));
+
+
+    std::vector<int>   row_indices;
+    
+    for (unsigned int row=0; row<n_rows; ++row)
+      if (row_map.MyGID(row))
+       {
+         const int row_length = sp.row_length(row);
+         row_indices.resize (row_length, -1);
+
+         for (int col=0; col < row_length; ++col)
+           row_indices[col] = sp.column_number (row, col);
+
+         graph->Epetra_CrsGraph::InsertGlobalIndices (row, row_length, 
+                                                      &row_indices[0]);
+       }
+
+    compress();
+  }
+
+
+
+  template<>
+  void 
+  SparsityPattern::reinit (const Epetra_Map   &input_row_map,
+                          const Epetra_Map   &input_col_map,
+                          const CompressedSetSparsityPattern &sp)
+  {
+    graph.reset();
+
+    Assert (sp.n_rows() == 
+             static_cast<unsigned int>(input_row_map.NumGlobalElements()),
+           ExcDimensionMismatch (sp.n_rows(),
+                                 input_row_map.NumGlobalElements()));
+    Assert (sp.n_cols() ==
+             static_cast<unsigned int>(input_col_map.NumGlobalElements()),
+           ExcDimensionMismatch (sp.n_cols(),
+                                 input_col_map.NumGlobalElements()));
+
+    row_map = input_row_map;
+    col_map = input_col_map;
+
+    const unsigned int n_rows = sp.n_rows();
+
+    std::vector<int> n_entries_per_row(n_rows);
+
+    for (unsigned int row=0; row<n_rows; ++row)
+      n_entries_per_row[row] = sp.row_length(row);
+
+    graph = std::auto_ptr<Epetra_FECrsGraph>
+      (new Epetra_FECrsGraph(Copy, row_map, 
+                            n_entries_per_row[input_row_map.MinMyGID()], 
+                            false));
+
+    Assert (graph->NumGlobalRows() == (int)sp.n_rows(),
+           ExcDimensionMismatch (graph->NumGlobalRows(),
+                                 sp.n_rows()));
+
+
+    std::vector<int>   row_indices;
+    
+    for (unsigned int row=0; row<n_rows; ++row)
+      if (row_map.MyGID(row))
+       {
+         const int row_length = sp.row_length(row);
+         row_indices.resize (row_length, -1);
+
+         CompressedSetSparsityPattern::row_iterator col_num = 
+           sp.row_begin (row);
+
+         for (unsigned int col = 0; 
+              col_num != sp.row_end (row); 
+              ++col_num, ++col)
+           row_indices[col] = *col_num;
+
+         graph->Epetra_CrsGraph::InsertGlobalIndices (row, row_length, 
+                                                      &row_indices[0]);
+       }
+
+    compress();
+  }
+
+
+
+  /*  void
+  SparsityPattern::copy_from (const SparsityPattern &sp)
+  {
+    graph.reset();
+    row_map = sp.row_map;
+    col_map = sp.col_map;
+
+    graph = std::auto_ptr<Epetra_FECrsGraph> (new Epetra_FECrsGraph(*sp.graph));
+  }
+  */
+
+
+  template <typename SparsityType>
+  void 
+  SparsityPattern::copy_from (const SparsityType &sp)
+  {
+#ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+    Epetra_MpiComm    trilinos_communicator (MPI_COMM_WORLD);
+#else
+    Epetra_SerialComm trilinos_communicator;
+#endif
+
+    const Epetra_Map rows (sp.n_rows(), 0, trilinos_communicator);
+    const Epetra_Map columns (sp.n_cols(), 0, trilinos_communicator);
+
+    reinit (rows, columns, sp);
+  }
+
+
+
+  void
+  SparsityPattern::clear ()
+  {
+                                 // When we clear the matrix, reset
+                                 // the pointer and generate an
+                                 // empty matrix.
+    graph.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
+
+    col_map = row_map;
+
+    graph = std::auto_ptr<Epetra_FECrsGraph> 
+             (new Epetra_FECrsGraph(View, row_map, 0));
+
+    graph->FillComplete();
+
+    compressed = true;
+  }
+
+
+
+  void
+  SparsityPattern::compress ()
+  {
+                                 // flush buffers
+    int ierr;
+    ierr = graph->GlobalAssemble (col_map, row_map, true);
+    
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    ierr = graph->OptimizeStorage ();
+    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+
+    compressed = true;
+  }
+
+
+
+  bool
+  SparsityPattern::exists (const unsigned int i,
+                          const unsigned int j) const
+  {
+                                     // Extract local indices in
+                                     // the matrix.
+    int trilinos_i = graph->LRID(i), trilinos_j = graph->LRID(j);
+
+                                     // If the data is not on the
+                                     // present processor, we throw
+                                     // an exception. This is on of
+                                     // the two tiny differences to
+                                     // the el(i,j) call, which does
+                                     // not throw any assertions.
+    if (trilinos_i == -1)
+      {
+       return false;
+      }
+    else
+      {
+                                     // Check whether the matrix 
+                                     // already is transformed to
+                                     // local indices.
+       if (graph->Filled() == false)
+         {
+           int nnz_present = graph->NumGlobalIndices(i);
+           int nnz_extracted;
+           int *col_indices;
+
+                                     // Generate the view and make
+                                     // sure that we have not generated
+                                     // an error.
+           int ierr = graph->ExtractGlobalRowView(trilinos_i, nnz_extracted,
+                                                  col_indices);
+           Assert (ierr==0, ExcTrilinosError(ierr));
+           Assert (nnz_present == nnz_extracted,
+                   ExcDimensionMismatch(nnz_present, nnz_extracted));
+
+                                     // Search the index
+           int* el_find = std::find(col_indices, col_indices + nnz_present,
+                                    trilinos_j);
+
+           int local_col_index = (int)(el_find - col_indices);
+
+           if (local_col_index == nnz_present)
+             return false;
+         }
+       else
+         {
+                                     // Prepare pointers for extraction
+                                     // of a view of the row.
+           int nnz_present = graph->NumGlobalIndices(i);
+           int nnz_extracted;
+           int *col_indices;
+
+                                     // Generate the view and make
+                                     // sure that we have not generated
+                                     // an error.
+           int ierr = graph->ExtractMyRowView(trilinos_i, nnz_extracted,
+                                              col_indices);
+           Assert (ierr==0, ExcTrilinosError(ierr));
+
+           Assert (nnz_present == nnz_extracted,
+                   ExcDimensionMismatch(nnz_present, nnz_extracted));
+
+                                     // Search the index
+           int* el_find = std::find(col_indices, col_indices + nnz_present,
+                                    trilinos_j);
+
+           int local_col_index = (int)(el_find - col_indices);
+
+           if (local_col_index == nnz_present)
+             return false;
+         }
+      }
+
+    return true;
+  }
+
+
+
+  unsigned int
+  SparsityPattern::n_rows () const
+  {
+    int n_rows = graph -> NumGlobalRows();
+
+    return n_rows;
+  }
+
+
+
+  unsigned int
+  SparsityPattern::n_cols () const
+  {
+    int n_cols = graph -> NumGlobalCols();
+    return n_cols;
+  }
+
+
+
+  unsigned int
+  SparsityPattern::local_size () const
+  {
+    int n_rows = graph -> NumMyRows();
+
+    return n_rows;
+  }
+
+
+
+  std::pair<unsigned int, unsigned int>
+  SparsityPattern::local_range () const
+  {
+    unsigned int begin, end;
+    begin = graph -> RowMap().MinMyGID();
+    end = graph -> RowMap().MaxMyGID()+1;
+    
+    return std::make_pair (begin, end);
+  }
+
+
+
+  unsigned int
+  SparsityPattern::n_nonzero_elements () const
+  {
+    int nnz = graph->NumGlobalEntries();
+
+    return static_cast<unsigned int>(nnz);
+  }
+
+
+
+  unsigned int
+  SparsityPattern::max_entries_per_row () const
+  {
+    int nnz = graph->MaxRowDim();
+
+    return static_cast<unsigned int>(nnz);
+  }
+
+
+
+  unsigned int
+  SparsityPattern::row_length (const unsigned int row) const
+  {
+    Assert (row < n_rows(), ExcInternalError());
+
+                                 // get a representation of the
+                                 // present row
+    int ncols = -1;
+    int local_row = graph->LRID(row);
+
+                                 // on the processor who owns this
+                                 // row, we'll have a non-negative
+                                 // value.
+    if (local_row >= 0)
+      ncols = graph->NumMyIndices (local_row);
+
+    return static_cast<unsigned int>(ncols);
+  }
+
+
+
+  void
+  SparsityPattern::write_ascii ()
+  {
+    Assert (false, ExcNotImplemented());
+  }
+
+
+
+                                 // As of now, no particularly neat
+                                 // ouput is generated in case of
+                                 // multiple processors.
+  void
+  SparsityPattern::print (std::ostream &out) const
+  {
+    int * indices;
+    int num_entries;
+  
+    for (int i=0; i<graph->NumMyRows(); ++i)
+      {
+       graph->ExtractMyRowView (i, num_entries, indices);
+       for (int j=0; j<num_entries; ++j)
+         out << "(" << i << "," << indices[graph->GRID(j)] << ") " 
+             << std::endl;
+      }
+  
+    AssertThrow (out, ExcIO());
+  }
+
+
+
+
+  // explicit instantiations
+  //
+  template void
+  SparsityPattern::copy_from (const dealii::SparsityPattern &);
+  template void
+  SparsityPattern::copy_from (const CompressedSparsityPattern &);
+  template void
+  SparsityPattern::copy_from (const CompressedSetSparsityPattern &);
+  template void
+  SparsityPattern::copy_from (const CompressedSimpleSparsityPattern &);
+
+
+  template void
+  SparsityPattern::reinit (const Epetra_Map &,
+                          const dealii::SparsityPattern &);
+  template void
+  SparsityPattern::reinit (const Epetra_Map &,
+                          const CompressedSparsityPattern &);
+  template void
+  SparsityPattern::reinit (const Epetra_Map &,
+                          const CompressedSetSparsityPattern &);
+  template void
+  SparsityPattern::reinit (const Epetra_Map &,
+                          const CompressedSimpleSparsityPattern &);
+
+
+  template void
+  SparsityPattern::reinit (const Epetra_Map &,
+                          const Epetra_Map &,
+                          const dealii::SparsityPattern &);
+  template void
+  SparsityPattern::reinit (const Epetra_Map &,
+                          const Epetra_Map &,
+                          const CompressedSparsityPattern &);
+  template void
+  SparsityPattern::reinit (const Epetra_Map &,
+                          const Epetra_Map &,
+                          const CompressedSimpleSparsityPattern &);
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.