]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Updated Trilinos AMG preconditioner to interact with Trilinos wrappers for matrices...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 22 Aug 2008 21:56:05 +0000 (21:56 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 22 Aug 2008 21:56:05 +0000 (21:56 +0000)
git-svn-id: https://svn.dealii.org/trunk@16658 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 3ab5cfbc1b8ead1912c544c7b1a444d739cd8053..1ecf56904799cc56c464e440e9fedb0ac5d42977 100755 (executable)
 
 
 #include <base/config.h>
+#include <lac/exceptions.h>
 #include <lac/trilinos_vector.h>
 #include <lac/trilinos_sparse_matrix.h>
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
 
 #include <boost/shared_ptr.hpp>
 
@@ -24,7 +27,6 @@
 #ifdef DEAL_II_USE_TRILINOS
 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
 #  include <Epetra_MpiComm.h>
-#  include <mpi.h>
 #else
 #  include <Epetra_SerialComm.h>
 #endif
@@ -65,6 +67,12 @@ namespace TrilinosWrappers
  * fact that some of the entries in the preconditioner matrix are zero
  * and hence can be neglected.
  *
+ * The implementation is able to distinguish between matrices from 
+ * elliptic problems and convection dominated problems. We use the standard
+ * options provided by Trilinos ML for elliptic problems, except that we use a 
+ * Chebyshev smoother instead of a symmetric Gauss-Seidel smoother. 
+ * For most elliptic problems, Chebyshev is better than Gauss-Seidel (SSOR).
+ *
  * @author Martin Kronbichler, 2008
  */
   class PreconditionAMG : public Subscriptor
@@ -84,15 +92,59 @@ namespace TrilinosWrappers
                                        * Let Trilinos compute a
                                        * multilevel hierarchy for the
                                        * solution of a linear system
-                                       * with the given matrix.
+                                       * with the given matrix. The
+                                       * function uses the matrix 
+                                       * format specified in
+                                       * TrilinosWrappers::SparseMatrix.
+                                       */
+      void initialize (const SparseMatrix        &matrix,
+                      const std::vector<double> &null_space,
+                      const unsigned int         null_space_dimension,
+                      const bool                 higher_order_elements,
+                      const bool                 elliptic,
+                      const bool                 output_details);
+
+                                      /**
+                                       * Let Trilinos compute a
+                                       * multilevel hierarchy for the
+                                       * solution of a linear system
+                                       * with the given matrix. This
+                                       * function takes a deal.ii matrix
+                                       * and copies the content into a
+                                       * Trilinos matrix, so the function
+                                       * can be considered rather 
+                                       * inefficient.
                                        */
       void initialize (const dealii::SparseMatrix<double> &matrix,
                       const std::vector<double>  &null_space,
                       const unsigned int          null_space_dimension,
                       const bool                  higher_order_elements,
                       const bool                  elliptic,
-                      const bool                  output_details,
-                      const double                drop_tolerance = 1e-13);
+                      const bool                  output_details);
+
+                                      /**
+                                       * This function can be used 
+                                       * for a faster recalculation of
+                                       * the preconditioner construction 
+                                       * when the matrix entries 
+                                       * underlying the preconditioner 
+                                       * have changed, 
+                                       * but the matrix sparsity pattern
+                                       * has remained the same. What this
+                                       * function does is to take the 
+                                       * already generated coarsening 
+                                       * structure, compute the AMG 
+                                       * prolongation and restriction 
+                                       * according to a smoothed aggregation
+                                       * strategy and then builds the whole
+                                       * multilevel hiearchy. This function
+                                       * can be considerably faster in that
+                                       * case, since the coarsening pattern
+                                       * is usually the most difficult thing
+                                       * to do when setting up the
+                                       * AMG ML preconditioner.
+                                       */
+      void reinit ();
 
                                       /**
                                        * Multiply the source vector
@@ -101,6 +153,15 @@ namespace TrilinosWrappers
                                        * and return it in the
                                        * destination vector.
                                        */
+      void vmult (Vector        &dst,
+                 const Vector  &src) const;
+
+                                      /**
+                                       * Do the same as before, but 
+                                       * now use deal.II vectors instead
+                                       * of the ones provided in the
+                                       * Trilinos wrapper class.
+                                       */
       void vmult (dealii::Vector<double>       &dst,
                  const dealii::Vector<double> &src) const;
 
@@ -122,7 +183,7 @@ namespace TrilinosWrappers
                                        * A copy of the deal.II matrix
                                        * into Trilinos format.
                                        */
-      boost::shared_ptr<Epetra_CrsMatrix> Matrix;
+      boost::shared_ptr<SparseMatrix> Matrix;
   };
 }
 
index 54b078c972315ba21ce1869490deebfd53c7862b..f3e432bc2a9dbc16059628ae309320b24dfab32d 100755 (executable)
@@ -16,7 +16,8 @@
 
 #include <base/config.h>
 #include <base/subscriptor.h>
-#include <lac/compressed_sparsity_pattern.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/sparse_matrix.h>
 #include <lac/exceptions.h>
 #include <lac/trilinos_vector.h>
 
@@ -414,8 +415,8 @@ namespace TrilinosWrappers
                                        * not directly available, use one of the
                                        * other reinit functions.
                                         */
-      void reinit (const CompressedSparsityPattern &sparsity_pattern,
-                  const unsigned int               n_max_entries_per_row);
+      void reinit (const SparsityPattern  &sparsity_pattern,
+                  const unsigned int      n_max_entries_per_row);
 
                                        /**
                                         * This function initializes the
@@ -426,7 +427,22 @@ namespace TrilinosWrappers
                                        * of nonzeros from the sparsity
                                        * pattern internally.
                                         */
-      void reinit (const CompressedSparsityPattern &sparsity_pattern);
+      void reinit (const SparsityPattern &sparsity_pattern);
+
+                                       /**
+                                        * 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 whose 
+                                       * modulus is larger than the 
+                                       * threshold (so zeros in the 
+                                       * deal.II matrix can be filtered
+                                       * away).
+                                        */
+      void reinit (const Epetra_Map                     &input_map,
+                  const ::dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
+                  const double                          drop_tolerance=1e-15);
 
                                       /**
                                         * This function is similar to the
@@ -437,8 +453,8 @@ namespace TrilinosWrappers
                                        * when the matrix structure changes,
                                        * e.g. when the grid is refined.
                                         */
-      void reinit (const Epetra_Map                &input_map,
-                  const CompressedSparsityPattern &sparsity_pattern);
+      void reinit (const Epetra_Map       &input_map,
+                  const SparsityPattern  &sparsity_pattern);
 
                                       /**
                                         * This function is similar to the
@@ -449,9 +465,9 @@ namespace TrilinosWrappers
                                        * To be used e.g. for rectangular 
                                        * matrices after remeshing.
                                         */
-      void reinit (const Epetra_Map                &input_row_map,
-                  const Epetra_Map                &input_col_map,
-                  const CompressedSparsityPattern &sparsity_pattern);
+      void reinit (const Epetra_Map       &input_row_map,
+                  const Epetra_Map       &input_col_map,
+                  const SparsityPattern  &sparsity_pattern);
 
                                        /**
                                         * Release all memory and return
index 332241451fc75d1e907045bfa9587ac0dba2a6f8..ca338ed565c52fd04949cb3c0880d5663b4435e3 100755 (executable)
 
 
 #include <lac/trilinos_precondition_amg.h>
-#include <lac/vector.h>
-#include <lac/sparse_matrix.h>
-
 
 #ifdef DEAL_II_USE_TRILINOS
 
 #include <Epetra_Map.h>
-#include <Epetra_CrsMatrix.h>
 #include <Epetra_Vector.h>
 #include <Teuchos_ParameterList.hpp>
 #include <ml_include.h>
@@ -43,80 +39,30 @@ namespace TrilinosWrappers
   PreconditionAMG::~PreconditionAMG ()
   {}
 
-  
-  
+
+
   void
   PreconditionAMG::
-  initialize (const dealii::SparseMatrix<double> &matrix,
-             const std::vector<double>  &null_space,
-             const unsigned int          null_space_dimension,
-             const bool                  elliptic,
-             const bool                  higher_order_elements,
-             const bool                  output_details,
-             const double                drop_tolerance)
+  initialize (const SparseMatrix        &matrix,
+             const std::vector<double> &null_space,
+             const unsigned int         null_space_dimension,
+             const bool                 elliptic,
+             const bool                 higher_order_elements,
+             const bool                 output_details)
   {
-    Assert (drop_tolerance >= 0,
-           ExcMessage ("Drop tolerance must be a non-negative number."));
-  
     const unsigned int n_rows = matrix.m();
-  
-                                    // Init Epetra Matrix, avoid 
-                                    // storing the nonzero elements.
-    {
-      Map.reset (new Epetra_Map(n_rows, 0, communicator));
     
-      std::vector<int> row_lengths (n_rows);
-      for (dealii::SparseMatrix<double>::const_iterator p = matrix.begin();
-          p != matrix.end(); ++p)
-       if (std::abs(p->value()) > drop_tolerance)
-         ++row_lengths[p->row()];
-  
-      Matrix.reset (new Epetra_CrsMatrix(Copy, *Map, &row_lengths[0], true));
-  
-      const unsigned int max_nonzero_entries
-       = *std::max_element (row_lengths.begin(), row_lengths.end());
-  
-      std::vector<double> values(max_nonzero_entries, 0);
-      std::vector<int> row_indices(max_nonzero_entries);
-  
-      for (unsigned int row=0; row<n_rows; ++row)
-       {
-         unsigned int index = 0;
-         for (dealii::SparseMatrix<double>::const_iterator p = matrix.begin(row);
-              p != matrix.end(row); ++p)
-           if (std::abs(p->value()) > drop_tolerance)
-             {
-               row_indices[index] = p->column();
-               values[index]      = p->value();
-               ++index;
-             }
-
-         Assert (index == static_cast<unsigned int>(row_lengths[row]),
-                 ExcMessage("Filtering out zeros could not "
-                            "be successfully finished!"));
-  
-         Matrix->InsertGlobalValues(row, row_lengths[row],
-                                    &values[0], &row_indices[0]);
-       }
-      
-      Matrix->FillComplete();
-    }
-  
+    if (!Matrix)
+      {
+       Matrix = boost::shared_ptr<SparseMatrix>
+                         (const_cast<SparseMatrix*>(&matrix));
+       Map = boost::shared_ptr<Epetra_Map>
+                         (const_cast<Epetra_Map*>(&matrix.row_map));
+      }
+
                                     // Build the AMG preconditioner.
     Teuchos::ParameterList parameter_list;
   
-                                    // The implementation is able
-                                    // to distinguish between
-                                    // matrices from elliptic problems
-                                    // and convection dominated 
-                                    // problems. We use the standard
-                                    // options for elliptic problems,
-                                    // except that we use a 
-                                    // Chebyshev smoother instead
-                                    // of a symmetric Gauss-Seidel
-                                    // smoother. For most elliptic 
-                                    // problems, Chebyshev is better
-                                    // than Gauss-Seidel (SSOR).
     if (elliptic)
       {
        ML_Epetra::SetDefaults("SA",parameter_list);
@@ -148,20 +94,67 @@ namespace TrilinosWrappers
        parameter_list.set("null space: dimension", int(null_space_dimension));
        parameter_list.set("null space: vectors", (double *)&null_space[0]);
       }
-  
+
     multigrid_operator = boost::shared_ptr<ML_Epetra::MultiLevelPreconditioner>
-                        (new ML_Epetra::MultiLevelPreconditioner(*Matrix, parameter_list, true));
+                        (new ML_Epetra::MultiLevelPreconditioner(
+                                     *matrix.matrix, parameter_list, true));
 
     if (output_details)
       multigrid_operator->PrintUnused(0);
   }
 
+
+
+  void
+  PreconditionAMG::
+  initialize (const dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
+             const std::vector<double>  &null_space,
+             const unsigned int          null_space_dimension,
+             const bool                  elliptic,
+             const bool                  higher_order_elements,
+             const bool                  output_details)
+  {
+    const unsigned int n_rows = deal_ii_sparse_matrix.m();
+  
+                                    // Init Epetra Matrix, avoid 
+                                    // storing the nonzero elements.
+
+    Map.reset (new Epetra_Map(n_rows, 0, communicator));
+
+    Matrix.reset();
+    Matrix = boost::shared_ptr<SparseMatrix> (new SparseMatrix());
+
+    Matrix->reinit (*Map, deal_ii_sparse_matrix);
+    Matrix->compress();
+
+    initialize (*Matrix, null_space, null_space_dimension, elliptic,
+               higher_order_elements, output_details);
+  }
+  
+  
+  void PreconditionAMG::
+  reinit ()
+  {
+    multigrid_operator->ReComputePreconditioner();
+  }
+  
+  
+  
+  void PreconditionAMG::vmult (Vector        &dst,
+                              const Vector  &src) const
+  {
+    const int ierr = multigrid_operator->ApplyInverse (*src.vector, *dst.vector);
+
+    Assert (ierr == 0, ExcTrilinosError(ierr));
+  }
+
+
                                   // For the implementation of the
                                   // <code>vmult</code> function we
                                   // note that invoking a call of 
                                   // the Trilinos preconditioner 
                                   // requires us to use Epetra vectors
-                                  // as well. Luckily, it is sufficient
+                                  // as well. It is faster
                                   // to provide a view, i.e., feed 
                                   // Trilinos with a pointer to the
                                   // data, so we avoid copying the
@@ -176,6 +169,15 @@ namespace TrilinosWrappers
   void PreconditionAMG::vmult (dealii::Vector<double>       &dst,
                               const dealii::Vector<double> &src) const
   {
+    Assert (Map->SameAs(Matrix->row_map),
+           ExcMessage("The sparse matrix given to the preconditioner uses "
+                      "a map that is not compatible. Check ML preconditioner "
+                      "setup."));
+    Assert (Map->SameAs(Matrix->col_map),
+           ExcMessage("The sparse matrix given to the preconditioner uses "
+                      "a map that is not compatible. Check ML preconditioner "
+                      "setup."));
+    
     Epetra_Vector LHS (View, *Map, dst.begin());
     Epetra_Vector RHS (View, *Map, const_cast<double*>(src.begin()));
   
index 0e9789f71c39e8e6e326461151052aa77f1a1d03..e5b4862043c9da1381b0ccd7ddb8055c66aa7746 100755 (executable)
@@ -147,8 +147,8 @@ namespace TrilinosWrappers
 
   
   void
-  SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern,
-                       const unsigned int               n_max_entries_per_row)
+  SparseMatrix::reinit (const SparsityPattern &sparsity_pattern,
+                       const unsigned int     n_max_entries_per_row)
   {
 
     unsigned int n_rows = sparsity_pattern.n_rows();
@@ -187,14 +187,14 @@ namespace TrilinosWrappers
 
 
   void
-  SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern)
+  SparseMatrix::reinit (const SparsityPattern &sparsity_pattern)
   {
     unsigned int n_rows = sparsity_pattern.n_rows();
-    
+
     Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(),
            ExcDimensionMismatch (matrix->NumGlobalRows(),
                                  sparsity_pattern.n_rows()));
-    
+
                                 // Trilinos seems to have a bug for
                                 // rectangular matrices at this point,
                                 // so do not check for consistent 
@@ -204,10 +204,10 @@ namespace TrilinosWrappers
 //                               sparsity_pattern.n_cols()));
 
     std::vector<int> n_entries_per_row(n_rows);
-    
+
     for (unsigned int row=0; row<n_rows; ++row)
       n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
-    
+
     const unsigned int n_max_entries_per_row = *std::max_element (
                    &n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
 
@@ -217,30 +217,31 @@ namespace TrilinosWrappers
 
 
   void
-  SparseMatrix::reinit (const Epetra_Map                &input_map,
-                       const CompressedSparsityPattern &sparsity_pattern)
+  SparseMatrix::reinit (const Epetra_Map       &input_map,
+                       const SparsityPattern  &sparsity_pattern)
   {
+    matrix.reset();
 
     unsigned int n_rows = sparsity_pattern.n_rows();
-    matrix.reset();
-    row_map = input_map;
-    col_map = row_map;
-    
+
     Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
            ExcDimensionMismatch (input_map.NumGlobalElements(),
                                  sparsity_pattern.n_rows()));
     Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_cols(),
            ExcDimensionMismatch (input_map.NumGlobalElements(),
                                  sparsity_pattern.n_cols()));
-    
+
+    row_map = input_map;
+    col_map = row_map;
+
     std::vector<int> n_entries_per_row(n_rows);
-    
+
     for (unsigned int row=0; row<n_rows; ++row)
       n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
-    
+
     matrix = std::auto_ptr<Epetra_FECrsMatrix> (new Epetra_FECrsMatrix (
                                Copy, row_map, &n_entries_per_row[0], true));
-    
+
     const unsigned int n_max_entries_per_row = *std::max_element (
                    &n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
 
@@ -250,32 +251,33 @@ namespace TrilinosWrappers
 
 
   void
-  SparseMatrix::reinit (const Epetra_Map                &input_row_map,
-                       const Epetra_Map                &input_col_map,
-                       const CompressedSparsityPattern &sparsity_pattern)
+  SparseMatrix::reinit (const Epetra_Map       &input_row_map,
+                       const Epetra_Map       &input_col_map,
+                       const SparsityPattern  &sparsity_pattern)
   {
+    matrix.reset();
 
     unsigned int n_rows = sparsity_pattern.n_rows();
-    matrix.reset();
-    row_map = input_row_map;
-    col_map = input_col_map;
-    
+
     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()));
-    
+
+    row_map = input_row_map;
+    col_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[(int)row] = sparsity_pattern.row_length(row);
-    
+
     matrix = std::auto_ptr<Epetra_FECrsMatrix>
              (new Epetra_FECrsMatrix(Copy, row_map, col_map,
                                      &n_entries_per_row[0], true));
-    
+
     const unsigned int n_max_entries_per_row = *std::max_element (
                    &n_entries_per_row[0], &n_entries_per_row[n_rows-1]);
 
@@ -284,6 +286,63 @@ namespace TrilinosWrappers
 
 
 
+  void
+  SparseMatrix::reinit (const Epetra_Map                     &input_map,
+                       const ::dealii::SparseMatrix<double> &deal_ii_sparse_matrix,
+                       const double                          drop_tolerance)
+  {
+    matrix.reset();
+
+    unsigned int n_rows = deal_ii_sparse_matrix.m();
+
+    Assert (input_map.NumGlobalElements() == (int)n_rows,
+           ExcDimensionMismatch (input_map.NumGlobalElements(),
+                                 n_rows));
+    Assert (input_map.NumGlobalElements() == (int)deal_ii_sparse_matrix.n(),
+           ExcDimensionMismatch (input_map.NumGlobalElements(),
+                                 deal_ii_sparse_matrix.n()));
+    
+    row_map = input_map;
+    col_map = row_map;
+
+    std::vector<int> n_entries_per_row(n_rows);
+
+    for (unsigned int row=0; row<n_rows; ++row)
+      n_entries_per_row[(int)row] = 
+         deal_ii_sparse_matrix.get_sparsity_pattern().row_length(row);
+
+    matrix = std::auto_ptr<Epetra_FECrsMatrix>
+             (new Epetra_FECrsMatrix(Copy, row_map, col_map,
+                                     &n_entries_per_row[0], true));
+    
+    std::vector<double> values;
+    std::vector<int> row_indices;
+
+    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],-1);
+       
+       unsigned int index = 0;
+       for (dealii::SparseMatrix<double>::const_iterator  
+             p  = deal_ii_sparse_matrix.begin(row);
+             p != deal_ii_sparse_matrix.end(row); ++p)
+         if (std::abs(p->value()) > drop_tolerance)
+           {
+             row_indices[index] = p->column();
+             values[index]      = p->value();
+             ++index;
+           }
+       const int n_row_entries = index;
+
+       matrix->InsertGlobalValues(row, n_row_entries,
+                                  &values[0], &row_indices[0]);
+      }
+
+  }
+
+
+
   void
   SparseMatrix::clear ()
   {

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.