trilinos_precondition.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2015 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#include <deal.II/lac/trilinos_precondition.h>

#ifdef DEAL_II_WITH_TRILINOS

#  include <deal.II/lac/vector.h>
#  include <deal.II/lac/sparse_matrix.h>
#  include <deal.II/lac/trilinos_sparse_matrix.h>

DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#  include <Ifpack.h>
#  include <Ifpack_Chebyshev.h>
#  include <Teuchos_ParameterList.hpp>
#  include <Teuchos_RCP.hpp>
#  include <Epetra_MultiVector.h>
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS

DEAL_II_NAMESPACE_OPEN

namespace TrilinosWrappers
{

  PreconditionBase::PreconditionBase()
#ifdef DEAL_II_WITH_MPI
    :
    communicator (MPI_COMM_SELF)
#endif
  {}



  PreconditionBase::PreconditionBase(const PreconditionBase &base)
    :
    Subscriptor (),
    preconditioner (base.preconditioner),
#ifdef DEAL_II_WITH_MPI
    communicator (base.communicator),
#endif
    vector_distributor (new Epetra_Map(*base.vector_distributor))
  {}



  PreconditionBase::~PreconditionBase()
  {}



  void PreconditionBase::clear ()
  {
    preconditioner.reset();
#ifdef DEAL_II_WITH_MPI
    communicator = MPI_COMM_SELF;
#endif
    vector_distributor.reset();
  }



  /* -------------------------- PreconditionJacobi -------------------------- */

  PreconditionJacobi::AdditionalData::
  AdditionalData (const double omega,
                  const double min_diagonal,
                  const unsigned int n_sweeps)
    :
    omega (omega),
    min_diagonal (min_diagonal),
    n_sweeps     (n_sweeps)
  {}



  void
  PreconditionJacobi::initialize (const SparseMatrix   &matrix,
                                  const AdditionalData &additional_data)
  {
    // release memory before reallocation
    preconditioner.reset ();
    preconditioner.reset (Ifpack().Create
                          ("point relaxation",
                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
                           0));

    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                    (preconditioner.get());
    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
                                     "preconditioner"));

    int ierr;

    Teuchos::ParameterList parameter_list;
    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
    parameter_list.set ("relaxation: type", "Jacobi");
    parameter_list.set ("relaxation: damping factor", additional_data.omega);
    parameter_list.set ("relaxation: min diagonal value",
                        additional_data.min_diagonal);

    ierr = ifpack->SetParameters(parameter_list);
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Initialize();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Compute();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
  }



  /* -------------------------- PreconditionSSOR -------------------------- */

  PreconditionSSOR::AdditionalData::
  AdditionalData (const double       omega,
                  const double       min_diagonal,
                  const unsigned int overlap,
                  const unsigned int n_sweeps)
    :
    omega        (omega),
    min_diagonal (min_diagonal),
    overlap      (overlap),
    n_sweeps     (n_sweeps)
  {}



  void
  PreconditionSSOR::initialize (const SparseMatrix   &matrix,
                                const AdditionalData &additional_data)
  {
    preconditioner.reset ();
    preconditioner.reset (Ifpack().Create
                          ("point relaxation",
                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
                           additional_data.overlap));

    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                    (preconditioner.get());
    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
                                     "preconditioner"));

    int ierr;

    Teuchos::ParameterList parameter_list;
    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
    parameter_list.set ("relaxation: type", "symmetric Gauss-Seidel");
    parameter_list.set ("relaxation: damping factor", additional_data.omega);
    parameter_list.set ("relaxation: min diagonal value",
                        additional_data.min_diagonal);
    parameter_list.set ("schwarz: combine mode", "Add");

    ierr = ifpack->SetParameters(parameter_list);
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Initialize();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Compute();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
  }



  /* -------------------------- PreconditionSOR -------------------------- */

  PreconditionSOR::AdditionalData::
  AdditionalData (const double       omega,
                  const double       min_diagonal,
                  const unsigned int overlap,
                  const unsigned int n_sweeps)
    :
    omega        (omega),
    min_diagonal (min_diagonal),
    overlap      (overlap),
    n_sweeps     (n_sweeps)
  {}



  void
  PreconditionSOR::initialize (const SparseMatrix   &matrix,
                               const AdditionalData &additional_data)
  {
    preconditioner.reset ();
    preconditioner.reset (Ifpack().Create
                          ("point relaxation",
                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
                           additional_data.overlap));

    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                    (preconditioner.get());
    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
                                     "preconditioner"));

    int ierr;

    Teuchos::ParameterList parameter_list;
    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
    parameter_list.set ("relaxation: type", "Gauss-Seidel");
    parameter_list.set ("relaxation: damping factor", additional_data.omega);
    parameter_list.set ("relaxation: min diagonal value",
                        additional_data.min_diagonal);
    parameter_list.set ("schwarz: combine mode", "Add");

    ierr = ifpack->SetParameters(parameter_list);
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Initialize();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Compute();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
  }



  /* ----------------------- PreconditionBlockJacobi ---------------------- */

  PreconditionBlockJacobi::AdditionalData::
  AdditionalData (const unsigned int block_size,
                  const std::string  block_creation_type,
                  const double omega,
                  const double min_diagonal,
                  const unsigned int n_sweeps)
    :
    block_size(block_size),
    block_creation_type(block_creation_type),
    omega (omega),
    min_diagonal (min_diagonal),
    n_sweeps     (n_sweeps)
  {}



  void
  PreconditionBlockJacobi::initialize (const SparseMatrix   &matrix,
                                       const AdditionalData &additional_data)
  {
    // release memory before reallocation
    preconditioner.reset ();
    preconditioner.reset (Ifpack().Create
                          ("block relaxation",
                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
                           0));

    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                    (preconditioner.get());
    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
                                     "preconditioner"));

    int ierr;

    Teuchos::ParameterList parameter_list;
    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
    parameter_list.set ("relaxation: type", "Jacobi");
    parameter_list.set ("relaxation: damping factor", additional_data.omega);
    parameter_list.set ("relaxation: min diagonal value",
                        additional_data.min_diagonal);
    parameter_list.set ("partitioner: type", additional_data.block_creation_type);
    int n_local_parts = (matrix.trilinos_matrix().NumMyRows()+additional_data.
                         block_size-1)/additional_data.block_size;
    parameter_list.set ("partitioner: local parts", n_local_parts);

    ierr = ifpack->SetParameters(parameter_list);
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Initialize();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Compute();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
  }



  /* ----------------------- PreconditionBlockSSOR ------------------------ */

  PreconditionBlockSSOR::AdditionalData::
  AdditionalData (const unsigned int block_size,
                  const std::string  block_creation_type,
                  const double       omega,
                  const double       min_diagonal,
                  const unsigned int overlap,
                  const unsigned int n_sweeps)
    :
    block_size(block_size),
    block_creation_type(block_creation_type),
    omega        (omega),
    min_diagonal (min_diagonal),
    overlap      (overlap),
    n_sweeps     (n_sweeps)
  {}



  void
  PreconditionBlockSSOR::initialize (const SparseMatrix   &matrix,
                                     const AdditionalData &additional_data)
  {
    preconditioner.reset ();
    preconditioner.reset (Ifpack().Create
                          ("block relaxation",
                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
                           additional_data.overlap));

    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                    (preconditioner.get());
    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
                                     "preconditioner"));

    int ierr;

    Teuchos::ParameterList parameter_list;
    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
    parameter_list.set ("relaxation: type", "symmetric Gauss-Seidel");
    parameter_list.set ("relaxation: damping factor", additional_data.omega);
    parameter_list.set ("relaxation: min diagonal value",
                        additional_data.min_diagonal);
    parameter_list.set ("schwarz: combine mode", "Add");
    parameter_list.set ("partitioner: type", additional_data.block_creation_type);
    int n_local_parts = (matrix.trilinos_matrix().NumMyRows()+additional_data.
                         block_size-1)/additional_data.block_size;
    parameter_list.set ("partitioner: local parts", n_local_parts);

    ierr = ifpack->SetParameters(parameter_list);
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Initialize();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Compute();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
  }



  /* ------------------------ PreconditionBlockSOR ------------------------ */

  PreconditionBlockSOR::AdditionalData::
  AdditionalData (const unsigned int block_size,
                  const std::string  block_creation_type,
                  const double       omega,
                  const double       min_diagonal,
                  const unsigned int overlap,
                  const unsigned int n_sweeps)
    :
    block_size(block_size),
    block_creation_type(block_creation_type),
    omega        (omega),
    min_diagonal (min_diagonal),
    overlap      (overlap),
    n_sweeps     (n_sweeps)
  {}



  void
  PreconditionBlockSOR::initialize (const SparseMatrix   &matrix,
                                    const AdditionalData &additional_data)
  {
    preconditioner.reset ();
    preconditioner.reset (Ifpack().Create
                          ("block relaxation",
                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
                           additional_data.overlap));

    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                    (preconditioner.get());
    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
                                     "preconditioner"));

    int ierr;

    Teuchos::ParameterList parameter_list;
    parameter_list.set ("relaxation: sweeps", static_cast<int>(additional_data.n_sweeps));
    parameter_list.set ("relaxation: type", "Gauss-Seidel");
    parameter_list.set ("relaxation: damping factor", additional_data.omega);
    parameter_list.set ("relaxation: min diagonal value",
                        additional_data.min_diagonal);
    parameter_list.set ("schwarz: combine mode", "Add");
    parameter_list.set ("partitioner: type", additional_data.block_creation_type);
    int n_local_parts = (matrix.trilinos_matrix().NumMyRows()+additional_data.
                         block_size-1)/additional_data.block_size;
    parameter_list.set ("partitioner: local parts", n_local_parts);

    ierr = ifpack->SetParameters(parameter_list);
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Initialize();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Compute();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
  }



  /* -------------------------- PreconditionIC -------------------------- */

  PreconditionIC::AdditionalData::
  AdditionalData (const unsigned int ic_fill,
                  const double       ic_atol,
                  const double       ic_rtol,
                  const unsigned int overlap)
    :
    ic_fill (ic_fill),
    ic_atol (ic_atol),
    ic_rtol (ic_rtol),
    overlap (overlap)
  {}



  void
  PreconditionIC::initialize (const SparseMatrix   &matrix,
                              const AdditionalData &additional_data)
  {
    preconditioner.reset ();
    preconditioner.reset (Ifpack().Create
                          ("IC",
                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
                           additional_data.overlap));

    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                    (preconditioner.get());
    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
                                     "preconditioner"));

    int ierr;

    Teuchos::ParameterList parameter_list;
    parameter_list.set ("fact: level-of-fill",(int)additional_data.ic_fill);
    parameter_list.set ("fact: absolute threshold",additional_data.ic_atol);
    parameter_list.set ("fact: relative threshold",additional_data.ic_rtol);
    parameter_list.set ("schwarz: combine mode", "Add");

    ierr = ifpack->SetParameters(parameter_list);
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Initialize();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Compute();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
  }



  /* -------------------------- PreconditionILU -------------------------- */

  PreconditionILU::AdditionalData::
  AdditionalData (const unsigned int ilu_fill,
                  const double       ilu_atol,
                  const double       ilu_rtol,
                  const unsigned int overlap)
    :
    ilu_fill (ilu_fill),
    ilu_atol (ilu_atol),
    ilu_rtol (ilu_rtol),
    overlap  (overlap)
  {}



  void
  PreconditionILU::initialize (const SparseMatrix   &matrix,
                               const AdditionalData &additional_data)
  {
    preconditioner.reset ();
    preconditioner.reset (Ifpack().Create
                          ("ILU",
                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
                           additional_data.overlap));

    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                    (preconditioner.get());
    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
                                     "preconditioner"));

    int ierr;

    Teuchos::ParameterList parameter_list;
    parameter_list.set ("fact: level-of-fill", static_cast<int>(additional_data.ilu_fill));
    parameter_list.set ("fact: absolute threshold", additional_data.ilu_atol);
    parameter_list.set ("fact: relative threshold", additional_data.ilu_rtol);
    parameter_list.set ("schwarz: combine mode", "Add");

    ierr = ifpack->SetParameters(parameter_list);
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Initialize();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Compute();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
  }



  /* -------------------------- PreconditionILUT -------------------------- */

  PreconditionILUT::AdditionalData::
  AdditionalData (const double       ilut_drop,
                  const unsigned int ilut_fill,
                  const double       ilut_atol,
                  const double       ilut_rtol,
                  const unsigned int overlap)
    :
    ilut_drop (ilut_drop),
    ilut_fill (ilut_fill),
    ilut_atol (ilut_atol),
    ilut_rtol (ilut_rtol),
    overlap  (overlap)
  {}



  void
  PreconditionILUT::initialize (const SparseMatrix   &matrix,
                                const AdditionalData &additional_data)
  {
    preconditioner.reset ();
    preconditioner.reset (Ifpack().Create
                          ("ILUT",
                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
                           additional_data.overlap));

    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                    (preconditioner.get());
    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
                                     "preconditioner"));

    int ierr;

    Teuchos::ParameterList parameter_list;
    parameter_list.set ("fact: drop value",additional_data.ilut_drop);
    parameter_list.set ("fact: level-of-fill",(int)additional_data.ilut_fill);
    parameter_list.set ("fact: absolute threshold",additional_data.ilut_atol);
    parameter_list.set ("fact: relative threshold",additional_data.ilut_rtol);
    parameter_list.set ("schwarz: combine mode", "Add");

    ierr = ifpack->SetParameters(parameter_list);
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Initialize();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Compute();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
  }



  /* ---------------------- PreconditionBlockDirect --------------------- */

  PreconditionBlockwiseDirect::AdditionalData::
  AdditionalData (const unsigned int overlap)
    :
    overlap  (overlap)
  {}



  void
  PreconditionBlockwiseDirect::initialize (const SparseMatrix   &matrix,
                                           const AdditionalData &additional_data)
  {
    preconditioner.reset ();
    preconditioner.reset (Ifpack().Create
                          ("Amesos",
                           const_cast<Epetra_CrsMatrix *>(&matrix.trilinos_matrix()),
                           additional_data.overlap));

    Ifpack_Preconditioner *ifpack = static_cast<Ifpack_Preconditioner *>
                                    (preconditioner.get());
    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
                                     "preconditioner"));

    int ierr;

    Teuchos::ParameterList parameter_list;
    parameter_list.set ("schwarz: combine mode", "Add");

    ierr = ifpack->SetParameters(parameter_list);
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Initialize();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Compute();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
  }



  /* ---------------------- PreconditionBlockDirect --------------------- */

  PreconditionChebyshev::AdditionalData::
  AdditionalData (const unsigned int degree,
                  const double       max_eigenvalue,
                  const double       eigenvalue_ratio,
                  const double       min_eigenvalue,
                  const double       min_diagonal,
                  const bool         nonzero_starting)
    :
    degree  (degree),
    max_eigenvalue (max_eigenvalue),
    eigenvalue_ratio (eigenvalue_ratio),
    min_eigenvalue (min_eigenvalue),
    min_diagonal (min_diagonal),
    nonzero_starting (nonzero_starting)
  {}



  void
  PreconditionChebyshev::initialize (const SparseMatrix   &matrix,
                                     const AdditionalData &additional_data)
  {
    preconditioner.reset ();
    preconditioner.reset (new Ifpack_Chebyshev (&matrix.trilinos_matrix()));

    Ifpack_Chebyshev *ifpack = static_cast<Ifpack_Chebyshev *>
                               (preconditioner.get());
    Assert (ifpack != 0, ExcMessage ("Trilinos could not create this "
                                     "preconditioner"));

    int ierr;

    Teuchos::ParameterList parameter_list;
    parameter_list.set ("chebyshev: ratio eigenvalue",
                        additional_data.eigenvalue_ratio);
    parameter_list.set ("chebyshev: min eigenvalue",
                        additional_data.min_eigenvalue);
    parameter_list.set ("chebyshev: max eigenvalue",
                        additional_data.max_eigenvalue);
    parameter_list.set ("chebyshev: degree",
                        (int)additional_data.degree);
    parameter_list.set ("chebyshev: min diagonal value",
                        additional_data.min_diagonal);
    parameter_list.set ("chebyshev: zero starting solution",
                        !additional_data.nonzero_starting);

    ierr = ifpack->SetParameters(parameter_list);
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Initialize();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));

    ierr = ifpack->Compute();
    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
  }








  /* -------------------------- PreconditionIdentity --------------------- */

  void
  PreconditionIdentity::vmult(VectorBase       &dst,
                              const VectorBase &src) const
  {
    dst = src;
  }

  void
  PreconditionIdentity::Tvmult(VectorBase       &dst,
                               const VectorBase &src) const
  {
    dst = src;
  }

  void
  PreconditionIdentity::vmult(::Vector<double>       &dst,
                              const ::Vector<double> &src) const
  {
    dst = src;
  }

  void
  PreconditionIdentity::Tvmult(::Vector<double>       &dst,
                               const ::Vector<double> &src) const
  {
    dst = src;
  }

  void
  PreconditionIdentity::vmult(parallel::distributed::Vector<double>       &dst,
                              const parallel::distributed::Vector<double> &src) const
  {
    dst = src;
  }

  void
  PreconditionIdentity::Tvmult(parallel::distributed::Vector<double>       &dst,
                               const parallel::distributed::Vector<double> &src) const
  {
    dst = src;
  }
}

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_TRILINOS