]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
More flexibility with initializing a Trilinos sparse matrix from a deal.II sparse...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 5 Feb 2010 11:53:27 +0000 (11:53 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 5 Feb 2010 11:53:27 +0000 (11:53 +0000)
git-svn-id: https://svn.dealii.org/trunk@20503 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_precondition.h
deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/source/trilinos_precondition.cc
deal.II/lac/source/trilinos_sparse_matrix.cc

index 227f2fa38386a2d05a2077423ee79ec5bfc8dbf3..8e6df158451cc09cf5afedbfe8e75706b004801f 100644 (file)
@@ -49,6 +49,7 @@ DEAL_II_NAMESPACE_OPEN
 // forward declarations
 template <typename number> class SparseMatrix;
 template <typename number> class Vector;
+class SparsityPattern;
 
 /*! @addtogroup TrilinosWrappers
  *@{
@@ -1465,7 +1466,8 @@ namespace TrilinosWrappers
       template <typename number>
       void initialize (const ::dealii::SparseMatrix<number> &deal_ii_sparse_matrix,
                       const AdditionalData                 &additional_data = AdditionalData(),
-                      const double                          drop_tolerance = 1e-13);
+                      const double                          drop_tolerance = 1e-13,
+                      const ::dealii::SparsityPattern      *use_this_sparsity = 0);
 
                                       /**
                                        * This function can be used for a
@@ -1492,6 +1494,12 @@ namespace TrilinosWrappers
                                        */
       void reinit ();
 
+                                      /**
+                                       * Prints an estimate of the memory
+                                       * consumption of this class.
+                                       */
+      unsigned int memory_consumption () const;
+
     private:
 
                                       /**
index 6fb4356d7fa8905ac9476bcad09bd22a6fe414b2..2fabc8216eaa14bae2fe7e814ce5d16a3d195688 100644 (file)
@@ -45,6 +45,7 @@ DEAL_II_NAMESPACE_OPEN
 template <typename MatrixType> class BlockMatrixBase;
 
 template <typename number> class SparseMatrix;
+class SparsityPattern;
 
 namespace TrilinosWrappers
 {
@@ -474,7 +475,8 @@ namespace TrilinosWrappers
       template <typename number>
       void reinit (const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
                   const double                          drop_tolerance=1e-13,
-                  const bool                            copy_values=true);
+                  const bool                            copy_values=true,
+                  const ::dealii::SparsityPattern      *use_this_sparsity=0);
 
                                        /**
                                        * This reinit function takes as
@@ -485,93 +487,6 @@ namespace TrilinosWrappers
                                        */
       void reinit (const Epetra_CrsMatrix &input_matrix,
                   const bool              copy_values = true);
-
-                                       /**
-                                        * 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.
-                                       *
-                                       * This is a collective operation
-                                       * that needs to be called on all
-                                       * processors in order to avoid a
-                                       * dead lock.
-                                        */
-      void clear ();
-
-                                       /**
-                                       * This command does two things:
-                                       * <ul>
-                                       * <li> If the matrix was initialized
-                                       * without a sparsity pattern,
-                                       * elements have been added manually
-                                       * using the set() command. When this
-                                       * process is completed, a call to
-                                       * compress() reorganizes the
-                                       * internal data structures (aparsity
-                                       * pattern) so that a fast access to
-                                       * data is possible in matrix-vector
-                                       * products.
-                                       * <li> If the matrix structure has
-                                       * already been fixed (either by
-                                       * initialization with a sparsity
-                                       * pattern or by calling compress()
-                                       * during the setup phase), this
-                                       * command does the %parallel
-                                       * exchange of data. This is
-                                       * necessary when we perform assembly
-                                       * on more than one (MPI) process,
-                                       * because then some non-local row
-                                       * data will accumulate on nodes that
-                                       * belong to the current's processor
-                                       * element, but are actually held by
-                                       * another. This command is usually
-                                       * called after all elements have
-                                       * been traversed.
-                                       * </ul>
-                                       *
-                                        * In both cases, this function
-                                        * compresses the data structures and
-                                        * allows the resulting matrix to be
-                                        * used in all other operations like
-                                        * matrix-vector products. This is a
-                                        * collective operation, i.e., it
-                                        * needs to be run on all processors
-                                        * when used in %parallel.
-                                       *
-                                       * See @ref GlossCompress "Compressing distributed objects"
-                                       * for more information.
-                                        */
-      void compress ();
-
-                                      /**
-                                       * Returns the state of the matrix,
-                                       * i.e., whether compress() needs to
-                                       * be called after an operation
-                                       * requiring data exchange. A call to
-                                       * compress() is also needed when the
-                                       * method set() has been called (even
-                                       * when working in serial).
-                                       */
-      bool is_compressed () const;
 //@}
 /**
  * @name Constructors and initialization using an Epetra_Map description
@@ -787,7 +702,8 @@ namespace TrilinosWrappers
       void reinit (const Epetra_Map                     &parallel_partitioning,
                   const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
                   const double                          drop_tolerance=1e-13,
-                  const bool                            copy_values=true);
+                  const bool                            copy_values=true,
+                  const ::dealii::SparsityPattern      *use_this_sparsity=0);
 
                                       /**
                                         * This function is similar to the
@@ -815,7 +731,8 @@ namespace TrilinosWrappers
                   const Epetra_Map                      &col_parallel_partitioning,
                   const ::dealii::SparseMatrix<number>  &dealii_sparse_matrix,
                   const double                           drop_tolerance=1e-13,
-                  const bool                             copy_values=true);
+                  const bool                             copy_values=true,
+                  const ::dealii::SparsityPattern      *use_this_sparsity=0);
 //@}
 /**
  * @name Constructors and initialization using an IndexSet description
@@ -1039,7 +956,8 @@ namespace TrilinosWrappers
                   const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
                   const MPI_Comm                       &communicator = MPI_COMM_WORLD,
                   const double                          drop_tolerance=1e-13,
-                  const bool                            copy_values=true);
+                  const bool                            copy_values=true,
+                  const ::dealii::SparsityPattern      *use_this_sparsity=0);
 
                                       /**
                                         * This function is similar to the
@@ -1068,7 +986,8 @@ namespace TrilinosWrappers
                   const ::dealii::SparseMatrix<number>  &dealii_sparse_matrix,
                   const MPI_Comm                        &communicator = MPI_COMM_WORLD,
                   const double                           drop_tolerance=1e-13,
-                  const bool                             copy_values=true);
+                  const bool                             copy_values=true,
+                  const ::dealii::SparsityPattern      *use_this_sparsity=0);
 //@}
 /**
  * @name Information on the matrix
@@ -1141,6 +1060,17 @@ namespace TrilinosWrappers
                                         */
       unsigned int row_length (const unsigned int row) const;
 
+                                      /**
+                                       * Returns the state of the matrix,
+                                       * i.e., whether compress() needs to
+                                       * be called after an operation
+                                       * requiring data exchange. A call to
+                                       * compress() is also needed when the
+                                       * method set() has been called (even
+                                       * when working in serial).
+                                       */
+      bool is_compressed () const;
+
                                       /**
                                        * Determine an estimate for the memory
                                        * consumption (in bytes) of this
@@ -1156,6 +1086,83 @@ namespace TrilinosWrappers
  * @name Modifying entries
  */
 //@{
+
+                                       /**
+                                        * 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.
+                                       *
+                                       * This is a collective operation
+                                       * that needs to be called on all
+                                       * processors in order to avoid a
+                                       * dead lock.
+                                        */
+      void clear ();
+
+                                       /**
+                                       * This command does two things:
+                                       * <ul>
+                                       * <li> If the matrix was initialized
+                                       * without a sparsity pattern,
+                                       * elements have been added manually
+                                       * using the set() command. When this
+                                       * process is completed, a call to
+                                       * compress() reorganizes the
+                                       * internal data structures (aparsity
+                                       * pattern) so that a fast access to
+                                       * data is possible in matrix-vector
+                                       * products.
+                                       * <li> If the matrix structure has
+                                       * already been fixed (either by
+                                       * initialization with a sparsity
+                                       * pattern or by calling compress()
+                                       * during the setup phase), this
+                                       * command does the %parallel
+                                       * exchange of data. This is
+                                       * necessary when we perform assembly
+                                       * on more than one (MPI) process,
+                                       * because then some non-local row
+                                       * data will accumulate on nodes that
+                                       * belong to the current's processor
+                                       * element, but are actually held by
+                                       * another. This command is usually
+                                       * called after all elements have
+                                       * been traversed.
+                                       * </ul>
+                                       *
+                                        * In both cases, this function
+                                        * compresses the data structures and
+                                        * allows the resulting matrix to be
+                                        * used in all other operations like
+                                        * matrix-vector products. This is a
+                                        * collective operation, i.e., it
+                                        * needs to be run on all processors
+                                        * when used in %parallel.
+                                       *
+                                       * See @ref GlossCompress "Compressing distributed objects"
+                                       * for more information.
+                                        */
+      void compress ();
+
                                        /**
                                         * Set the element (<i>i,j</i>)
                                         * to @p value.
@@ -3091,10 +3098,12 @@ namespace TrilinosWrappers
                             const ::dealii::SparseMatrix<number> &sparse_matrix,
                             const MPI_Comm      &communicator,
                             const double         drop_tolerance,
-                            const bool           copy_values)
+                            const bool           copy_values,
+                            const ::dealii::SparsityPattern *use_this_sparsity)
   {
     Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator, false);
-    reinit (map, map, sparse_matrix, drop_tolerance, copy_values);
+    reinit (map, map, sparse_matrix, drop_tolerance, copy_values, 
+           use_this_sparsity);
   }
 
 
@@ -3106,13 +3115,15 @@ namespace TrilinosWrappers
                             const ::dealii::SparseMatrix<number> &sparse_matrix,
                             const MPI_Comm      &communicator,
                             const double         drop_tolerance,
-                            const bool           copy_values)
+                            const bool           copy_values,
+                            const ::dealii::SparsityPattern *use_this_sparsity)
   {
     Epetra_Map row_map =
       row_parallel_partitioning.make_trilinos_map (communicator, false);
     Epetra_Map col_map =
       col_parallel_partitioning.make_trilinos_map (communicator, false);
-    reinit (row_map, col_map, sparse_matrix, drop_tolerance, copy_values);
+    reinit (row_map, col_map, sparse_matrix, drop_tolerance, copy_values,
+           use_this_sparsity);
   }
 
 
index 25c18e22f4e2b67f44aa6f108f8553fb5f2d7177..39aeb5217de72c92e708e1e491b0515ff938c0a9 100644 (file)
@@ -645,7 +645,8 @@ namespace TrilinosWrappers
   PreconditionAMG::
   initialize (const ::dealii::SparseMatrix<number> &deal_ii_sparse_matrix,
              const AdditionalData                 &additional_data,
-             const double                          drop_tolerance)
+             const double                          drop_tolerance,
+             const ::dealii::SparsityPattern      *use_this_sparsity)
   {
     const unsigned int n_rows = deal_ii_sparse_matrix.m();
 
@@ -658,26 +659,41 @@ namespace TrilinosWrappers
     Matrix.reset();
     Matrix = std_cxx1x::shared_ptr<SparseMatrix> (new SparseMatrix());
 
-    Matrix->reinit (*map, deal_ii_sparse_matrix, drop_tolerance);
-    Matrix->compress();
+    Matrix->reinit (*map, *map, deal_ii_sparse_matrix, drop_tolerance, true,
+                   use_this_sparsity);
 
     initialize (*Matrix, additional_data);
   }
 
 
+
   void PreconditionAMG::reinit ()
   {
     multilevel_operator->ReComputePreconditioner();
   }
 
 
+  unsigned int PreconditionAMG::memory_consumption() const
+  {
+    unsigned int memory = sizeof(this);
+
+                               // todo: find a way to read out ML's data
+                               // sizes
+    if (Matrix != 0)
+      memory += Matrix->memory_consumption();
+    return memory;
+  }
+
+
 
 
                                // explicit instantiations
   template void PreconditionAMG::initialize (const ::dealii::SparseMatrix<double> &,
-                                            const AdditionalData &, const double);
+                                            const AdditionalData &, const double,
+                                            const ::dealii::SparsityPattern*);
   template void PreconditionAMG::initialize (const ::dealii::SparseMatrix<float> &,
-                                            const AdditionalData &, const double);
+                                            const AdditionalData &, const double,
+                                            const ::dealii::SparsityPattern*);
 
 }
 
index 363c8a1788084cff7277e80b6a86df97d2dc5ccd..cc952399c2bcf47d61e23639103b2c89e074f6b2 100644 (file)
@@ -465,12 +465,14 @@ namespace TrilinosWrappers
          const int row_length = sparsity_pattern.row_length(row);
          row_indices.resize (row_length, -1);
 
-         for (int col=0; col < row_length; ++col)
-           row_indices[col] = sparsity_pattern.column_number (row, col);
-
-         graph->InsertGlobalIndices (static_cast<int>(row),
-                                     row_length, &row_indices[0]);
+         typename SparsityType::row_iterator col_num = 
+           sparsity_pattern.row_begin (row),
+           row_end = sparsity_pattern.row_end(row);
+         for (unsigned int col = 0; col_num != row_end; ++col_num, ++col)
+           row_indices[col] = *col_num;
 
+         graph->Epetra_CrsGraph::InsertGlobalIndices (row, row_length,
+                                                      &row_indices[0]);
        }
 
                                  // Eventually, optimize the graph
@@ -514,93 +516,6 @@ namespace TrilinosWrappers
 
 
 
-                                  // The CompressedSetSparsityPattern
-                                  // class stores the columns
-                                  // differently, so we need to
-                                  // explicitly rewrite this function
-                                  // in that case.
-  template<>
-  void
-  SparseMatrix::reinit (const Epetra_Map                    &input_row_map,
-                       const Epetra_Map                    &input_col_map,
-                       const CompressedSetSparsityPattern  &sparsity_pattern,
-                       const bool                           exchange_data)
-  {
-                                  // this function is similar to the other
-                                  // reinit function with sparsity pattern
-                                  // argument. See there for additional
-                                  // comments.
-    matrix.reset();
-
-    const unsigned int n_rows = sparsity_pattern.n_rows();
-
-    Assert (exchange_data == false, ExcInternalError());
-    if (input_row_map.Comm().MyPID() == 0)
-      {
-       Assert (input_row_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
-               ExcDimensionMismatch (input_row_map.NumGlobalElements(),
-                                     sparsity_pattern.n_rows()));
-       Assert (input_col_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
-               ExcDimensionMismatch (input_col_map.NumGlobalElements(),
-                                     sparsity_pattern.n_cols()));
-      }
-
-    column_space_map = std::auto_ptr<Epetra_Map> (new Epetra_Map (input_col_map));
-
-    std::vector<int> n_entries_per_row(n_rows);
-
-    for (unsigned int row=0; row<n_rows; ++row)
-      n_entries_per_row[row] = sparsity_pattern.row_length(row);
-
-    std::auto_ptr<Epetra_CrsGraph> graph;
-    if (input_row_map.Comm().NumProc() > 1)
-      graph = std::auto_ptr<Epetra_CrsGraph>
-       (new Epetra_CrsGraph (Copy, input_row_map,
-                             &n_entries_per_row[input_row_map.MinMyGID()], true));
-    else
-      graph = std::auto_ptr<Epetra_CrsGraph>
-       (new Epetra_CrsGraph (Copy, input_row_map, input_col_map,
-                             &n_entries_per_row[input_row_map.MinMyGID()], true));
-
-    Assert (graph->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
-           ExcDimensionMismatch (graph->NumGlobalRows(),
-                                 sparsity_pattern.n_rows()));
-
-    std::vector<int>   row_indices;
-
-    for (unsigned int row=0; row<n_rows; ++row)
-      if (input_row_map.MyGID(row))
-       {
-         const int row_length = sparsity_pattern.row_length(row);
-         row_indices.resize (row_length, -1);
-
-         CompressedSetSparsityPattern::row_iterator col_num =
-           sparsity_pattern.row_begin (row);
-
-         for (unsigned int col = 0;
-              col_num != sparsity_pattern.row_end (row);
-              ++col_num, ++col)
-           row_indices[col] = *col_num;
-
-         graph->InsertGlobalIndices (row, row_length, &row_indices[0]);
-       }
-
-    graph->FillComplete(input_col_map, input_row_map);
-    graph->OptimizeStorage();
-
-    Assert (graph->NumGlobalCols() == (int)sparsity_pattern.n_cols(),
-           ExcDimensionMismatch (graph->NumGlobalCols(),
-                                 sparsity_pattern.n_cols()));
-
-    matrix =  std::auto_ptr<Epetra_FECrsMatrix>
-                       (new Epetra_FECrsMatrix(Copy, *graph, false));
-
-    last_action = Zero;
-    compress();
-  }
-
-
-
   void
   SparseMatrix::reinit (const SparsityPattern &sparsity_pattern)
   {
@@ -639,7 +554,8 @@ namespace TrilinosWrappers
   void
   SparseMatrix::reinit (const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
                        const double                          drop_tolerance,
-                       const bool                            copy_values)
+                       const bool                            copy_values,
+                       const ::dealii::SparsityPattern      *use_this_sparsity)
   {
     const Epetra_Map rows (dealii_sparse_matrix.m(),
                           0,
@@ -647,7 +563,8 @@ namespace TrilinosWrappers
     const Epetra_Map columns (dealii_sparse_matrix.n(),
                              0,
                              Utilities::Trilinos::comm_self());
-    reinit (rows, columns, dealii_sparse_matrix, drop_tolerance, copy_values);
+    reinit (rows, columns, dealii_sparse_matrix, drop_tolerance, 
+           copy_values, use_this_sparsity);
   }
 
 
@@ -657,10 +574,11 @@ namespace TrilinosWrappers
   SparseMatrix::reinit (const Epetra_Map                     &input_map,
                        const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
                        const double                          drop_tolerance,
-                       const bool                            copy_values)
+                       const bool                            copy_values,
+                       const ::dealii::SparsityPattern      *use_this_sparsity)
   {
     reinit (input_map, input_map, dealii_sparse_matrix, drop_tolerance,
-           copy_values);
+           copy_values, use_this_sparsity);
   }
 
 
@@ -671,7 +589,8 @@ namespace TrilinosWrappers
                        const Epetra_Map                     &input_col_map,
                        const ::dealii::SparseMatrix<number> &dealii_sparse_matrix,
                        const double                          drop_tolerance,
-                       const bool                            copy_values)
+                       const bool                            copy_values,
+                       const ::dealii::SparsityPattern      *use_this_sparsity)
   {
     matrix.reset();
 
@@ -679,8 +598,12 @@ namespace TrilinosWrappers
       {
                                   // in case we do not copy values, just
                                   // call the other function.
-       reinit (input_row_map, input_col_map,
-               dealii_sparse_matrix.get_sparsity_pattern());
+       if (use_this_sparsity == 0)
+         reinit (input_row_map, input_col_map,
+                 dealii_sparse_matrix.get_sparsity_pattern());
+       else
+         reinit (input_row_map, input_col_map,
+                 *use_this_sparsity);
        return;
       }
 
@@ -695,39 +618,114 @@ namespace TrilinosWrappers
 
     column_space_map = std::auto_ptr<Epetra_Map> (new Epetra_Map (input_col_map));
 
+    const ::dealii::SparsityPattern & sparsity_pattern = 
+      (use_this_sparsity!=0)? *use_this_sparsity : 
+      dealii_sparse_matrix.get_sparsity_pattern();
+
     std::vector<int> n_entries_per_row(n_rows);
 
+    unsigned int maximum_row_length = 0;
     for (unsigned int row=0; row<n_rows; ++row)
-      n_entries_per_row[(int)row] =
-         dealii_sparse_matrix.get_sparsity_pattern().row_length(row);
+      {
+       const unsigned int row_length = sparsity_pattern.row_length(row);
+       if (row_length > maximum_row_length)
+         maximum_row_length = row_length;
+       n_entries_per_row[row] = row_length;
+      }
 
-    matrix = std::auto_ptr<Epetra_FECrsMatrix>
-               (new Epetra_FECrsMatrix(Copy, input_row_map,
-                                       &n_entries_per_row[input_row_map.MinMyGID()],
-                                       false));
+    std::auto_ptr<Epetra_CrsGraph> graph;
+    if (input_row_map.Comm().NumProc() > 1)
+      graph = std::auto_ptr<Epetra_CrsGraph>
+       (new Epetra_CrsGraph (Copy, input_row_map,
+                             &n_entries_per_row[input_row_map.MinMyGID()], false));
+    else
+      graph = std::auto_ptr<Epetra_CrsGraph>
+       (new Epetra_CrsGraph (Copy, input_row_map, input_col_map,
+                             &n_entries_per_row[input_row_map.MinMyGID()], false));
 
-    std::vector<TrilinosScalar> values;
-    std::vector<unsigned int>   row_indices;
+    Assert (graph->NumGlobalRows() == static_cast<int>(n_rows),
+           ExcDimensionMismatch (graph->NumGlobalRows(), n_rows));
 
+                                  // now insert the indices
+    std::vector<int>   row_indices (maximum_row_length);
     for (unsigned int row=0; row<n_rows; ++row)
-      {
-       values.resize (n_entries_per_row[row],0.);
-       row_indices.resize (n_entries_per_row[row],
-                           numbers::invalid_unsigned_int);
-
-       unsigned int index = 0;
-       for (typename ::dealii::SparseMatrix<number>::const_iterator
-             p  = dealii_sparse_matrix.begin(row);
-             p != dealii_sparse_matrix.end(row); ++p)
-         if (std::fabs(p->value()) > drop_tolerance)
+      if (input_row_map.MyGID(row))
+       {
+         typename ::dealii::SparsityPattern::row_iterator 
+           col_num = sparsity_pattern.row_begin (row),
+           row_end = sparsity_pattern.row_end(row);
+         typename ::dealii::SparseMatrix<number>::const_iterator
+           matrix_col_num = dealii_sparse_matrix.begin(row),
+           matrix_row_end = dealii_sparse_matrix.end(row);
+
+         unsigned int col = 0;
+         if (sparsity_pattern.optimize_diagonal())
+           {
+             row_indices[col++] = row;
+             ++col_num;
+             ++ matrix_col_num;
+           }
+         while (col_num != row_end && matrix_col_num != matrix_row_end)
            {
-             row_indices[index] = p->column();
-             values[index]      = p->value();
-             ++index;
+             while (*col_num < matrix_col_num->column() && col_num != row_end)
+               ++col_num;
+             while (*col_num > matrix_col_num->column() && matrix_col_num != matrix_row_end)
+               ++matrix_col_num;
+
+             if (std::fabs(matrix_col_num->value()) > drop_tolerance)
+               row_indices[col++] = *col_num;
+             ++col_num;
+             ++matrix_col_num;
            }
 
-       set (row, index, &row_indices[0], &values[0], false);
-      }
+         graph->Epetra_CrsGraph::InsertGlobalIndices (row, col,
+                                                      &row_indices[0]);
+       }
+    graph->FillComplete(input_col_map, input_row_map);
+    graph->OptimizeStorage();
+
+    matrix =  std::auto_ptr<Epetra_FECrsMatrix>
+                       (new Epetra_FECrsMatrix(Copy, *graph, false));
+    compress();
+
+    std::vector<TrilinosScalar> values (maximum_row_length);
+
+    for (unsigned int row=0; row<n_rows; ++row)
+      if (input_row_map.MyGID(row))
+       {
+         typename ::dealii::SparsityPattern::row_iterator 
+           col_num = sparsity_pattern.row_begin (row),
+           row_end = sparsity_pattern.row_end(row);
+         typename ::dealii::SparseMatrix<number>::const_iterator
+           matrix_col_num = dealii_sparse_matrix.begin(row),
+           matrix_row_end = dealii_sparse_matrix.end(row);
+
+         unsigned int col = 0;
+         if (sparsity_pattern.optimize_diagonal())
+           {
+             values[col] = matrix_col_num->value();
+             row_indices[col++] = row;
+             ++col_num;
+             ++ matrix_col_num;
+           }
+         while (col_num != row_end && matrix_col_num != matrix_row_end)
+           {
+             while (*col_num < matrix_col_num->column() && col_num != row_end)
+               ++col_num;
+             while (*col_num > matrix_col_num->column() && matrix_col_num != matrix_row_end)
+               ++matrix_col_num;
+
+             if (std::fabs(matrix_col_num->value()) > drop_tolerance)
+               {
+                 values[col] = matrix_col_num->value();
+                 row_indices[col++] = *col_num;
+               }
+             ++col_num;
+             ++matrix_col_num;
+           }
+         set (row, col, reinterpret_cast<unsigned int*>(&row_indices[0]), 
+              &values[0], false);
+       }
 
     compress();
   }
@@ -1467,54 +1465,73 @@ namespace TrilinosWrappers
                        const Epetra_Map &,
                        const CompressedSparsityPattern &,
                        const bool);
+  template void
+  SparseMatrix::reinit (const Epetra_Map &,
+                       const Epetra_Map &,
+                       const CompressedSimpleSparsityPattern &,
+                       const bool);
+  template void
+  SparseMatrix::reinit (const Epetra_Map &,
+                       const Epetra_Map &,
+                       const CompressedSetSparsityPattern &,
+                       const bool);
 
   template void
   SparseMatrix::reinit (const dealii::SparseMatrix<float> &,
                        const double,
-                       const bool);
+                       const bool,
+                       const dealii::SparsityPattern *);
   template void
   SparseMatrix::reinit (const dealii::SparseMatrix<double> &,
                        const double,
-                       const bool);
+                       const bool,
+                       const dealii::SparsityPattern *);
   template void
   SparseMatrix::reinit (const dealii::SparseMatrix<long double> &,
                        const double,
-                       const bool);
+                       const bool,
+                       const dealii::SparsityPattern *);
 
   template void
   SparseMatrix::reinit (const Epetra_Map &,
                        const dealii::SparseMatrix<float> &,
                        const double,
-                       const bool);
+                       const bool,
+                       const dealii::SparsityPattern *);
   template void
   SparseMatrix::reinit (const Epetra_Map &,
                        const dealii::SparseMatrix<double> &,
                        const double,
-                       const bool);
+                       const bool,
+                       const dealii::SparsityPattern *);
   template void
   SparseMatrix::reinit (const Epetra_Map &,
                        const dealii::SparseMatrix<long double> &,
                        const double,
-                       const bool);
+                       const bool,
+                       const dealii::SparsityPattern *);
 
   template void
   SparseMatrix::reinit (const Epetra_Map &,
                        const Epetra_Map &,
                        const dealii::SparseMatrix<float> &,
                        const double,
-                       const bool);
+                       const bool,
+                       const dealii::SparsityPattern *);
   template void
   SparseMatrix::reinit (const Epetra_Map &,
                        const Epetra_Map &,
                        const dealii::SparseMatrix<double> &,
                        const double,
-                       const bool);
+                       const bool,
+                       const dealii::SparsityPattern *);
   template void
   SparseMatrix::reinit (const Epetra_Map &,
                        const Epetra_Map &,
                        const dealii::SparseMatrix<long double> &,
                        const double,
-                       const bool);
+                       const bool,
+                       const dealii::SparsityPattern *);
 
 
 }

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.