general_purpose_preconditioners.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Include guards
#ifndef OOMPH_GENERAL_PRECONDITION_HEADER
#define OOMPH_GENERAL_PRECONDITION_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#include "preconditioner.h"
#include "iterative_linear_solver.h"
#include "matrices.h"
#include "problem.h"
#include <algorithm>
 
 
namespace oomph
{
  //=====================================================================
  //=====================================================================
  class MatrixBasedDiagPreconditioner : public Preconditioner
  {
  public:
    MatrixBasedDiagPreconditioner(){};
 
    ~MatrixBasedDiagPreconditioner(){};
 
    MatrixBasedDiagPreconditioner(const MatrixBasedDiagPreconditioner&) =
      delete;
 
    void operator=(const MatrixBasedDiagPreconditioner&) = delete;
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    void setup();
 
  private:
    Vector<double> Inv_diag;
  };
 
  //=============================================================================
  //=============================================================================
  template<typename MATRIX>
  class MatrixBasedLumpedPreconditioner : public Preconditioner
  {
  public:
    MatrixBasedLumpedPreconditioner()
    {
      // default the positive matrix boolean to false
      Positive_matrix = false;
 
      // set the pointers to the lumped vector to 0
      Inv_lumped_diag_pt = 0;
    };
 
    ~MatrixBasedLumpedPreconditioner()
    {
      this->clean_up_memory();
    }
 
    MatrixBasedLumpedPreconditioner(const MatrixBasedDiagPreconditioner&) =
      delete;
 
    void operator=(const MatrixBasedLumpedPreconditioner&) = delete;
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    void setup();
 
    using Preconditioner::setup;
 
    bool positive_matrix() const
    {
#ifdef PARANOID
      if (Inv_lumped_diag_pt == 0)
      {
        throw OomphLibError("The preconditioner has not been setup.",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Positive_matrix;
    }
 
 
    double* inverse_lumped_vector_pt()
    {
#ifdef PARANOID
      if (Inv_lumped_diag_pt == 0)
      {
        throw OomphLibError("The inverse lumped vector has not been created. "
                            "Created in setup(...)",
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Inv_lumped_diag_pt;
    }
 
 
    unsigned& nrow()
    {
      return Nrow;
    }
 
    void clean_up_memory()
    {
      delete[] Inv_lumped_diag_pt;
    }
 
  private:
    double* Inv_lumped_diag_pt;
 
    // boolean to indicate whether the lumped matrix was positive
    bool Positive_matrix;
 
    // number of rows in preconditioner
    unsigned Nrow;
  };
 
 
  //=============================================================================
  //=============================================================================
  class CompressedMatrixCoefficient
  {
  public:
    CompressedMatrixCoefficient() {}
 
    CompressedMatrixCoefficient(const unsigned& index, const double& value)
    {
      // store the index and value
      Index = index;
      Value = value;
    }
 
    ~CompressedMatrixCoefficient() {}
 
    CompressedMatrixCoefficient(const CompressedMatrixCoefficient& a)
    {
      Index = a.index();
      Value = a.value();
    }
 
    void operator=(const CompressedMatrixCoefficient& a)
    {
      Index = a.index();
      Value = a.value();
    }
 
    bool operator<(const CompressedMatrixCoefficient& a) const
    {
      return Index < a.index();
    }
 
    unsigned& index()
    {
      return Index;
    }
 
    double& value()
    {
      return Value;
    }
 
    unsigned index() const
    {
      return Index;
    }
 
    double value() const
    {
      return Value;
    }
 
  private:
    unsigned Index;
 
    double Value;
  };
 
 
  //=============================================================================
  //=============================================================================
  template<typename MATRIX>
  class ILUZeroPreconditioner : public Preconditioner
  {
  };
 
 
  //=============================================================================
  //=============================================================================
  template<>
  class ILUZeroPreconditioner<CCDoubleMatrix> : public Preconditioner
  {
  public:
    ILUZeroPreconditioner(){};
 
    ~ILUZeroPreconditioner(){};
 
 
    ILUZeroPreconditioner(const ILUZeroPreconditioner&) = delete;
 
    void operator=(const ILUZeroPreconditioner&) = delete;
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    void setup();
 
  private:
    Vector<unsigned> U_column_start;
 
    Vector<CompressedMatrixCoefficient> U_row_entry;
 
    Vector<unsigned> L_column_start;
 
    Vector<CompressedMatrixCoefficient> L_row_entry;
  };
 
 
  //=============================================================================
  //=============================================================================
  template<>
  class ILUZeroPreconditioner<CRDoubleMatrix> : public Preconditioner
  {
  public:
    ILUZeroPreconditioner(){};
 
 
    ILUZeroPreconditioner(const ILUZeroPreconditioner&) = delete;
 
    void operator=(const ILUZeroPreconditioner&) = delete;
 
    ~ILUZeroPreconditioner(){};
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    void setup();
 
 
  private:
    Vector<unsigned> U_row_start;
 
    Vector<CompressedMatrixCoefficient> U_row_entry;
 
    Vector<unsigned> L_row_start;
 
    Vector<CompressedMatrixCoefficient> L_row_entry;
  };
 
  //=============================================================================
  //=============================================================================
  template<class SOLVER, class PRECONDITIONER>
  class InnerIterationPreconditioner : public Preconditioner
  {
  public:
    InnerIterationPreconditioner()
    {
      // create the solver
      Solver_pt = new SOLVER;
 
      // create the preconditioner
      Preconditioner_pt = new PRECONDITIONER;
 
#ifdef PARANOID
      // paranoid check that the solver is an iterative solver and
      // the preconditioner is a preconditioner
      if (dynamic_cast<IterativeLinearSolver*>(Solver_pt) == 0)
      {
        throw OomphLibError(
          "The template argument SOLVER must be of type IterativeLinearSolver",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
      if (dynamic_cast<Preconditioner*>(Preconditioner_pt) == 0)
      {
        throw OomphLibError(
          "The template argument PRECONDITIONER must be of type Preconditioner",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // ensure the solver does not re-setup the preconditioner
      Solver_pt->disable_setup_preconditioner_before_solve();
 
      // pass the preconditioner to the solver
      Solver_pt->preconditioner_pt() = Preconditioner_pt;
    }
 
    // destructor
    ~InnerIterationPreconditioner()
    {
      delete Solver_pt;
      delete Preconditioner_pt;
    }
 
    // clean the memory
    void clean_up_memory()
    {
      Preconditioner_pt->clean_up_memory();
      Solver_pt->clean_up_memory();
    }
 
    void setup()
    {
      // set the distribution
      DistributableLinearAlgebraObject* dist_pt =
        dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt());
      if (dist_pt != 0)
      {
        this->build_distribution(dist_pt->distribution_pt());
      }
      else
      {
        LinearAlgebraDistribution dist(comm_pt(), matrix_pt()->nrow(), false);
        this->build_distribution(dist);
      }
 
      // Setup the inner iteration preconditioner (For some reason we need to
      // remind the compiler that there is also a function named setup in the
      // base class.)
      Preconditioner_pt->Preconditioner::setup(matrix_pt());
 
      // setup the solverready for resolve
      unsigned max_iter = Solver_pt->max_iter();
      Solver_pt->max_iter() = 1;
      DoubleVector x(this->distribution_pt(), 0.0);
      DoubleVector y(x);
      Solver_pt->enable_resolve();
      Solver_pt->solve(matrix_pt(), x, y);
      Solver_pt->max_iter() = max_iter;
    }
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z)
    {
      Solver_pt->resolve(r, z);
    }
 
    double& tolerance()
    {
      return Solver_pt->tolerance();
    }
 
    unsigned& max_iter()
    {
      return Solver_pt->max_iter();
    }
 
    SOLVER* solver_pt()
    {
      return Solver_pt;
    }
 
    PRECONDITIONER* preconditioner_pt()
    {
      return Preconditioner_pt;
    }
 
  private:
    SOLVER* Solver_pt;
 
    PRECONDITIONER* Preconditioner_pt;
  };
} // namespace oomph
#endif