Implement new initialize method for Trilinos AMG preconditioner: It can be based...

diff --git a/include/deal.II/lac/trilinos_precondition.h b/include/deal.II/lac/trilinos_precondition.h
index 11ecf25339f347acf2510b3fe6cf4d8e8126059e..97c015c08fe3284b1972343a639f3314f3e5a4a8 100644 (file)
--- a/include/deal.II/lac/trilinos_precondition.h
+++ b/include/deal.II/lac/trilinos_precondition.h
@@ -36,7 +36,7 @@
 #  include <Epetra_Map.h>
 
 #  include <Teuchos_ParameterList.hpp>
-#  include <Epetra_Operator.h>
+#  include <Epetra_RowMatrix.h>
 #  include <Epetra_Vector.h>
 
 // forward declarations
@@ -1467,9 +1467,19 @@ namespace TrilinosWrappers
      * linear system with the given matrix. The function uses the matrix
      * format specified in TrilinosWrappers::SparseMatrix.
      */
-    void initialize (const SparseMatrix                    &matrix,
+    void initialize (const SparseMatrix   &matrix,
                      const AdditionalData &additional_data = AdditionalData());
 
+    /**
+     * Let Trilinos compute a multilevel hierarchy for the solution of a
+     * linear system with the given matrix. As opposed to the other initialize
+     * function above, this function uses an abstract interface to an object
+     * of type Epetra_RowMatrix which allows a user to pass by rather
+     * arbitrary objects to the ML preconditioner.
+     */
+    void initialize (const Epetra_RowMatrix &matrix,
+                     const AdditionalData   &additional_data = AdditionalData());
+
     /**
      * Let Trilinos compute a multilevel hierarchy for the solution of a
      * linear system with the given matrix. The function uses the matrix
@@ -1485,6 +1495,16 @@ namespace TrilinosWrappers
     void initialize (const SparseMatrix           &matrix,
                      const Teuchos::ParameterList &ml_parameters);
 
+    /**
+     * Let Trilinos compute a multilevel hierarchy for the solution of a
+     * linear system with the given matrix. As opposed to the other initialize
+     * function above, this function uses an abstract interface to an object
+     * of type Epetra_RowMatrix which allows a user to pass by rather
+     * arbitrary objects to the ML preconditioner.
+     */
+    void initialize (const Epetra_RowMatrix       &matrix,
+                     const Teuchos::ParameterList &ml_parameters);
+
     /**
      * Let Trilinos compute a multilevel hierarchy for the solution of a
      * linear system with the given matrix. This function takes a deal.ii

diff --git a/source/lac/trilinos_precondition.cc b/source/lac/trilinos_precondition.cc
index 13d093e154645883fad5a4ea1045089b2800f769..17451985dd75fdd6184cca5e0199e73348327cd0 100644 (file)
--- a/source/lac/trilinos_precondition.cc
+++ b/source/lac/trilinos_precondition.cc
@@ -1,7 +1,7 @@
 // ---------------------------------------------------------------------
 // $Id$
 //
-// Copyright (C) 2008 - 2013 by the deal.II authors
+// Copyright (C) 2008 - 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -37,6 +37,11 @@ namespace TrilinosWrappers
   namespace
   {
 #ifndef DEAL_II_WITH_64BIT_INDICES
+    long long int n_global_rows (const Epetra_RowMatrix &matrix)
+    {
+      return matrix.NumGlobalRows();
+    }
+
     int global_length (const Epetra_MultiVector &vector)
     {
       return vector.GlobalLength();
@@ -47,6 +52,11 @@ namespace TrilinosWrappers
       return map.GID(i);
     }
 #else
+    long long int n_global_rows (const Epetra_RowMatrix &matrix)
+    {
+      return matrix.NumGlobalRows64();
+    }
+
     long long int global_length (const Epetra_MultiVector &vector)
     {
       return vector.GlobalLength64();
@@ -728,12 +738,20 @@ namespace TrilinosWrappers
   }
 
 
+
   void
   PreconditionAMG:: initialize (const SparseMatrix   &matrix,
                                 const AdditionalData &additional_data)
   {
-    const size_type n_rows = matrix.m();
+    initialize(matrix.trilinos_matrix(), additional_data);
+  }
+
+
 
+  void
+  PreconditionAMG:: initialize (const Epetra_RowMatrix &matrix,
+                                const AdditionalData   &additional_data)
+  {
     // Build the AMG preconditioner.
     Teuchos::ParameterList parameter_list;
 
@@ -749,15 +767,8 @@ namespace TrilinosWrappers
         // standard choice uncoupled. if higher order, right now we also just
         // use Uncoupled, but we should be aware that maybe MIS might be
         // needed
-        //
-        // TODO: Maybe there are any other/better options?
         if (additional_data.higher_order_elements)
-          {
-            //if (matrix.m()/matrix.matrix->Comm().NumProc() < 50000)
-            parameter_list.set("aggregation: type", "Uncoupled");
-            //else
-            //parameter_list.set("aggregation: type", "MIS");
-          }
+          parameter_list.set("aggregation: type", "Uncoupled");
       }
     else
       {
@@ -790,7 +801,7 @@ namespace TrilinosWrappers
     else
       parameter_list.set("ML output", 0);
 
-    const Epetra_Map &domain_map = matrix.domain_partitioner();
+    const Epetra_Map &domain_map = matrix.OperatorDomainMap();
 
     const size_type constant_modes_dimension =
       additional_data.constant_modes.size();
@@ -801,6 +812,7 @@ namespace TrilinosWrappers
 
     if (constant_modes_dimension > 0)
       {
+        const size_type n_rows = n_global_rows(matrix);
         const bool constant_modes_are_global =
           additional_data.constant_modes[0].size() == n_rows;
         const size_type n_relevant_rows =
@@ -855,10 +867,19 @@ namespace TrilinosWrappers
   void
   PreconditionAMG::initialize (const SparseMatrix           &matrix,
                                const Teuchos::ParameterList &ml_parameters)
+  {
+    initialize(matrix.trilinos_matrix(), ml_parameters);
+  }
+
+
+
+  void
+  PreconditionAMG::initialize (const Epetra_RowMatrix       &matrix,
+                               const Teuchos::ParameterList &ml_parameters)
   {
     preconditioner.reset ();
     preconditioner.reset (new ML_Epetra::MultiLevelPreconditioner
-                          (matrix.trilinos_matrix(), ml_parameters));
+                          (matrix, ml_parameters));
   }