]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add the Trilinos algebraic multigrid preconditioner.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 21 Aug 2008 20:35:47 +0000 (20:35 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 21 Aug 2008 20:35:47 +0000 (20:35 +0000)
git-svn-id: https://svn.dealii.org/trunk@16641 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_preconditioner_amg.h [new file with mode: 0755]
deal.II/lac/source/trilinos_preconditioner_amg.cc [new file with mode: 0755]

diff --git a/deal.II/lac/include/lac/trilinos_preconditioner_amg.h b/deal.II/lac/include/lac/trilinos_preconditioner_amg.h
new file mode 100755 (executable)
index 0000000..9af75d8
--- /dev/null
@@ -0,0 +1,110 @@
+//---------------------------------------------------------------------------
+//    $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_precondition_amg_h
+#define __deal2__trilinos_precondition_amg_h
+
+
+#include <base/config.h>
+#include <lac/trilinos_vector.h>
+#include <lac/trilinos_sparse_matrix.h>
+
+#include <boost/shared_ptr.hpp>
+
+
+#ifdef DEAL_II_USE_TRILINOS
+
+// some forward declarations
+namespace ML_Epetra
+{
+  class MultiLevelPreconditioner;
+}
+class Epetra_Map;
+class Epetra_CrsMatrix;
+
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+                                  // This implements an AMG
+                                  // preconditioner based on the
+                                  // Trilinos ML implementation.
+                                  // What this class does is twofold.
+                                  // When the constructor of the class
+                                  // is invoked, a ML preconditioner
+                                  // object is created based on the
+                                  // DoFHandler and matrix
+                                  // that we want the preconditioner to
+                                  // be based on. A call of
+                                  // the respective <code>vmult</code>
+                                  // function does call the respective
+                                  // operation in the Trilinos package,
+                                  // where it is called 
+                                  // <code>ApplyInverse</code>.
+                                  // There are a few pecularities in
+                                  // the constructor. Since the
+                                  // Trilinos objects we want to use are
+                                  // heavily dependent on Epetra objects,
+                                  // the fundamental construction
+                                  // routines for vectors and 
+                                  // matrices in Trilinos, we do a 
+                                  // copy of our deal.II preconditioner
+                                  // matrix to a Epetra matrix. This 
+                                  // is of course not optimal, but for
+                                  // the time being there is no direct
+                                  // support for our data interface.
+                                  // When doing this time-consuming 
+                                  // operation, we can still profit 
+                                  // from the fact that some of the
+                                  // entries in the preconditioner matrix
+                                  // are zero and hence can be 
+                                  // neglected.
+  class PreconditionerTrilinosAmg : public Subscriptor
+  {
+    public:
+      PreconditionerTrilinosAmg ();
+      ~PreconditionerTrilinosAmg ();
+    
+      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);
+
+      void vmult (dealii::Vector<double>       &dst,
+                 const dealii::Vector<double> &src) const;
+
+    private:
+    
+      boost::shared_ptr<ML_Epetra::MultiLevelPreconditioner> multigrid_operator;
+    
+      Epetra_SerialComm  communicator;
+      boost::shared_ptr<Epetra_Map>       Map;
+      boost::shared_ptr<Epetra_CrsMatrix> Matrix;
+  };
+}
+
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS
+
+/*----------------------------   trilinos_precondition_amg_base.h     ---------------------------*/
+
+#endif
+/*----------------------------   trilinos_precondition_amg_base.h     ---------------------------*/
diff --git a/deal.II/lac/source/trilinos_preconditioner_amg.cc b/deal.II/lac/source/trilinos_preconditioner_amg.cc
new file mode 100755 (executable)
index 0000000..72cb844
--- /dev/null
@@ -0,0 +1,185 @@
+//---------------------------------------------------------------------------
+//    $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_preconditioner_amg.h>
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+
+#include <Epetra_SerialComm.h>
+#include <Epetra_Map.h>
+#include <Epetra_CrsMatrix.h>
+#include <Epetra_Vector.h>
+#include <Teuchos_ParameterList.hpp>
+#include <ml_include.h>
+#include <ml_MultiLevelPreconditioner.h>
+
+
+
+#ifdef DEAL_II_USE_TRILINOS
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace TrilinosWrappers
+{
+
+  PreconditionerTrilinosAmg::PreconditionerTrilinosAmg ()
+  {}
+
+  void PreconditionerTrilinosAmg::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
+  )
+  {
+    Assert (drop_tolerance >= 0,
+           ExcMessage ("Drop tolerance must be a non-negative number."));
+  
+    const unsigned int n_rows = matrix.m();
+    const SparsityPattern *sparsity_pattern = &(matrix.get_sparsity_pattern());
+  
+                                    // 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();
+    }
+  
+                                    // 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);
+       parameter_list.set("smoother: type", "Chebyshev");
+       parameter_list.set("smoother: sweeps", 4);
+      }
+    else
+      {
+       ML_Epetra::SetDefaults("NSSA",parameter_list);
+       parameter_list.set("aggregation: type", "Uncoupled");
+       parameter_list.set("aggregation: block scaling", true);
+      }
+  
+    if (output_details)
+      parameter_list.set("ML output", 10);
+    else
+      parameter_list.set("ML output", 0);
+  
+    if (higher_order_elements)
+      parameter_list.set("aggregation: type", "MIS");
+  
+    Assert (n_rows * null_space_dimension == null_space.size(),
+           ExcDimensionMismatch(n_rows * null_space_dimension,
+                                null_space.size()));
+  
+    if (null_space_dimension > 1)
+      {
+       parameter_list.set("null space: type", "pre-computed");
+       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));
+
+    if (output_details)
+      multigrid_operator->PrintUnused(0);
+  }
+
+                                  // 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
+                                  // to provide a view, i.e., feed 
+                                  // Trilinos with a pointer to the
+                                  // data, so we avoid copying the
+                                  // content of the vectors during
+                                  // the iteration. In the declaration
+                                  // of the right hand side, we need
+                                  // to cast the source vector (that
+                                  // is <code>const</code> in all deal.II 
+                                  // calls) to non-constant value, as
+                                  // this is the way Trilinos wants to
+                                  // have them.
+  void PreconditionerTrilinosAmg::vmult (dealii::Vector<double>       &dst,
+                                        const dealii::Vector<double> &src) const
+  {
+    Epetra_Vector LHS (View, *Map, dst.begin());
+    Epetra_Vector RHS (View, *Map, const_cast<double*>(src.begin()));
+  
+    const int res = multigrid_operator->ApplyInverse (RHS, LHS);
+  
+    Assert (res == 0,
+           ExcMessage ("Trilinos AMG MultiLevel preconditioner returned "
+                       "with an error!"));
+  }
+
+}
+
+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.