linear_solver.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//====================================================================
 
// This header defines a class for linear solvers
 
// Include guards
#ifndef OOMPH_LINEAR_SOLVER_HEADER
#define OOMPH_LINEAR_SOLVER_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif
 
// oomph-lib headers
#include "Vector.h"
#include "double_vector.h"
#include "matrices.h"
 
namespace oomph
{
  // Forward declaration of problem class
  class Problem;
 
  //====================================================================
  //====================================================================
  class LinearSolver : public DistributableLinearAlgebraObject
  {
  protected:
    bool Enable_resolve;
 
    // for the solves should be documented
    bool Doc_time;
 
    bool Compute_gradient;
 
    bool Gradient_has_been_computed;
 
    DoubleVector Gradient_for_glob_conv_newton_solve;
 
  public:
    LinearSolver()
      : Enable_resolve(false),
        Doc_time(true),
        Compute_gradient(false),
        Gradient_has_been_computed(false)
    {
    }
 
    LinearSolver(const LinearSolver& dummy) = delete;
 
    void operator=(const LinearSolver&) = delete;
 
    virtual ~LinearSolver() {}
 
    void enable_doc_time()
    {
      Doc_time = true;
    }
 
    void disable_doc_time()
    {
      Doc_time = false;
    }
 
    bool is_doc_time_enabled() const
    {
      return Doc_time;
    }
 
    bool is_resolve_enabled() const
    {
      return Enable_resolve;
    }
 
    virtual void enable_resolve()
    {
      Enable_resolve = true;
    }
 
    virtual void disable_resolve()
    {
      Enable_resolve = false;
    }
 
    virtual void solve(Problem* const& problem_pt, DoubleVector& result) = 0;
 
    virtual void solve(DoubleMatrixBase* const& matrix_pt,
                       const DoubleVector& rhs,
                       DoubleVector& result)
    {
      throw OomphLibError(
        "DoubleVector based solve function not implemented for this solver",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void solve(DoubleMatrixBase* const& matrix_pt,
                       const Vector<double>& rhs,
                       Vector<double>& result)
    {
      throw OomphLibError(
        "Vector<double> based solve function not implemented for this solver",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void solve_transpose(Problem* const& problem_pt,
                                 DoubleVector& result)
    {
      // Create an output stream
      std::ostringstream error_message_stream;
 
      // Create the error message
      error_message_stream << "The function to solve the transposed system has "
                           << "not yet been\nimplemented for this solver."
                           << std::endl;
 
      // Throw the error message
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    } // End of solve_transpose
 
    virtual void solve_transpose(DoubleMatrixBase* const& matrix_pt,
                                 const DoubleVector& rhs,
                                 DoubleVector& result)
    {
      throw OomphLibError(
        "DoubleVector based solve function not implemented for this solver",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void solve_transpose(DoubleMatrixBase* const& matrix_pt,
                                 const Vector<double>& rhs,
                                 Vector<double>& result)
    {
      throw OomphLibError(
        "Vector<double> based solve function not implemented for this solver",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void resolve(const DoubleVector& rhs, DoubleVector& result)
    {
      throw OomphLibError(
        "Resolve function not implemented for this linear solver",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void resolve_transpose(const DoubleVector& rhs,
                                   DoubleVector& result)
    {
      // Create an output stream
      std::ostringstream error_message_stream;
 
      // Create the error message
      error_message_stream
        << "The function to resolve the transposed system has "
        << "not yet been\nimplemented for this solver." << std::endl;
 
      // Throw the error message
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    } // End of resolve_transpose
 
    virtual void clean_up_memory() {}
 
    virtual double jacobian_setup_time() const
    {
      throw OomphLibError(
        "jacobian_setup_time has not been implemented for this linear solver",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      return 0;
    }
 
    virtual double linear_solver_solution_time() const
    {
      throw OomphLibError(
        "linear_solver_solution_time() not implemented for this linear solver",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
      return 0;
    }
 
    virtual void enable_computation_of_gradient()
    {
      throw OomphLibError(
        "enable_computation_of_gradient() not implemented for "
        "this linear solver",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    void disable_computation_of_gradient()
    {
      Compute_gradient = false;
    }
 
    void reset_gradient()
    {
      Gradient_for_glob_conv_newton_solve.clear();
    }
 
    void get_gradient(DoubleVector& gradient)
    {
#ifdef PARANOID
      if (Gradient_has_been_computed)
      {
#endif
        gradient = Gradient_for_glob_conv_newton_solve;
#ifdef PARANOID
      }
      else
      {
        throw OomphLibError(
          "The gradient has not been computed for this linear solver!",
          OOMPH_CURRENT_FUNCTION,
          OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
  };
 
  //=============================================================================
  //============================================================================
  class DenseLU : public LinearSolver
  {
    friend class DenseDoubleMatrix;
 
  public:
    DenseLU()
      : Jacobian_setup_time(0.0),
        Solution_time(0.0),
        Sign_of_determinant_of_matrix(0),
        Index(0),
        LU_factors(0)
    {
      // Shut up!
      Doc_time = false;
    }
 
    DenseLU(const DenseLU& dummy) = delete;
 
    void operator=(const DenseLU&) = delete;
 
    ~DenseLU()
    {
      clean_up_memory();
    }
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& result);
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result);
 
    double jacobian_setup_time() const
    {
      return Jacobian_setup_time;
    }
 
    virtual double linear_solver_solution_time() const
    {
      return Solution_time;
    }
 
  protected:
    void factorise(DoubleMatrixBase* const& matrix_pt);
 
    void backsub(const DoubleVector& rhs, DoubleVector& result);
 
    void backsub(const Vector<double>& rhs, Vector<double>& result);
 
    void clean_up_memory();
 
    double Jacobian_setup_time;
 
    double Solution_time;
 
    int Sign_of_determinant_of_matrix;
 
  private:
    long* Index;
 
    double* LU_factors;
  };
 
 
  //====================================================================
  //====================================================================
  class FD_LU : public DenseLU
  {
  public:
    FD_LU() : DenseLU() {}
 
    FD_LU(const FD_LU& dummy) = delete;
 
    void operator=(const FD_LU&) = delete;
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& result)
    {
      DenseLU::solve(matrix_pt, rhs, result);
    }
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const Vector<double>& rhs,
               Vector<double>& result)
    {
      LinearSolver::solve(matrix_pt, rhs, result);
    }
  };
 
 
 
 
  //=============================================================================
  //=============================================================================
  class SuperLUSolver : public LinearSolver
  {
  public:
    enum Type
    {
      Default,
      Serial,
      Distributed
    };
 
    SuperLUSolver() : Serial_f_factors(0)
    {
      // Set solver wide default values and null pointers
      Doc_stats = false;
      Suppress_solve = false;
      Using_dist = false;
      Solver_type = Default;
 
#ifdef OOMPH_HAS_MPI
      // Set default values and nullify pointers for SuperLU Dist
      Dist_use_global_solver = false;
      Dist_delete_matrix_data = false;
      Dist_allow_row_and_col_permutations = true;
      Dist_solver_data_pt = 0;
      Dist_global_solve_data_allocated = false;
      Dist_distributed_solve_data_allocated = false;
      Dist_value_pt = 0;
      Dist_index_pt = 0;
      Dist_start_pt = 0;
#endif
 
      // Set default values and null pointers for SuperLU (serial)
      Serial_compressed_row_flag = true;
      Serial_sign_of_determinant_of_matrix = 0;
      Serial_n_dof = 0;
    }
 
    SuperLUSolver(const SuperLUSolver& dummy) = delete;
 
    void operator=(const SuperLUSolver&) = delete;
 
    ~SuperLUSolver()
    {
      clean_up_memory();
    }
 
    void enable_computation_of_gradient()
    {
      Compute_gradient = true;
    }
 
    // SuperLUSolver methods
 
    void disable_resolve()
    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    }
 
    void solve(Problem* const& problem_pt, DoubleVector& result);
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& result);
 
 
    /*  /// Linear-algebra-type solver: Takes pointer to a matrix */
    /*  /// and rhs vector and returns the solution of the linear system */
    /*  /// Call the broken base-class version. If you want this, please  */
    /*  /// implement it */
    /*  void solve(DoubleMatrixBase* const &matrix_pt, */
    /*                     const Vector<double> &rhs, */
    /*                     Vector<double> &result) */
    /*   {LinearSolver::solve(matrix_pt,rhs,result);} */
 
    void solve_transpose(Problem* const& problem_pt, DoubleVector& result);
 
    void solve_transpose(DoubleMatrixBase* const& matrix_pt,
                         const DoubleVector& rhs,
                         DoubleVector& result);
 
    void resolve(const DoubleVector& rhs, DoubleVector& result);
 
    void resolve_transpose(const DoubleVector& rhs, DoubleVector& result);
 
    void enable_doc_stats()
    {
      Doc_stats = true;
    }
 
    void disable_doc_stats()
    {
      Doc_stats = false;
    }
 
    double jacobian_setup_time() const
    {
      return Jacobian_setup_time;
    }
 
    virtual double linear_solver_solution_time() const
    {
      return Solution_time;
    }
 
    void factorise(DoubleMatrixBase* const& matrix_pt);
 
    void backsub(const DoubleVector& rhs, DoubleVector& result);
 
    void backsub_transpose(const DoubleVector& rhs, DoubleVector& result);
 
    void clean_up_memory();
 
    void set_solver_type(const Type& t)
    {
      this->clean_up_memory();
      Solver_type = t;
    }
 
    // SuperLU (serial) methods
 
    void use_compressed_row_for_superlu_serial()
    {
      Serial_compressed_row_flag = true;
    }
 
    void use_compressed_column_for_superlu_serial()
    {
      Serial_compressed_row_flag = false;
    }
 
#ifdef OOMPH_HAS_MPI
 
    // SuperLU Dist methods
 
    void enable_delete_matrix_data_in_superlu_dist()
    {
      Dist_delete_matrix_data = true;
    }
 
    void disable_delete_matrix_data_in_superlu_dist()
    {
      Dist_delete_matrix_data = false;
    }
 
    void enable_row_and_col_permutations_in_superlu_dist()
    {
      Dist_allow_row_and_col_permutations = true;
    }
 
    void disable_row_and_col_permutations_in_superlu_dist()
    {
      Dist_allow_row_and_col_permutations = false;
    }
 
    void use_global_solve_in_superlu_dist()
    {
      if (!Dist_use_global_solver)
      {
        clean_up_memory();
        Dist_use_global_solver = true;
      }
    }
 
    void use_distributed_solve_in_superlu_dist()
    {
      if (Dist_use_global_solver)
      {
        clean_up_memory();
        Dist_use_global_solver = false;
      }
    }
 
#endif
 
  private:
    void factorise_serial(DoubleMatrixBase* const& matrix_pt);
 
    void backsub_serial(const DoubleVector& rhs, DoubleVector& result);
 
    void backsub_transpose_serial(const DoubleVector& rhs,
                                  DoubleVector& result);
 
#ifdef OOMPH_HAS_MPI
    void factorise_distributed(DoubleMatrixBase* const& matrix_pt);
 
    void backsub_distributed(const DoubleVector& rhs, DoubleVector& result);
 
    void backsub_transpose_distributed(const DoubleVector& rhs,
                                       DoubleVector& result);
#endif
 
    // SuperLUSolver member data
 
    double Jacobian_setup_time;
 
    double Solution_time;
 
    bool Suppress_solve;
 
    bool Doc_stats;
 
    Type Solver_type;
 
    bool Using_dist;
 
    // SuperLU (serial) member data
 
    void* Serial_f_factors;
 
    int Serial_info;
 
    unsigned long Serial_n_dof;
 
    int Serial_sign_of_determinant_of_matrix;
 
    bool Serial_compressed_row_flag;
 
  public:
    double get_memory_usage_for_lu_factors();
 
    double get_total_needed_memory();
 
  private:
#ifdef OOMPH_HAS_MPI
 
    // SuperLU Dist member data
  public:
    static bool Suppress_incorrect_rhs_distribution_warning_in_resolve;
 
  private:
    bool Dist_use_global_solver;
 
    bool Dist_global_solve_data_allocated;
 
    bool Dist_distributed_solve_data_allocated;
 
    void* Dist_solver_data_pt;
 
    int Dist_nprow;
 
    int Dist_npcol;
 
    int Dist_info;
 
    bool Dist_allow_row_and_col_permutations;
 
    bool Dist_delete_matrix_data;
 
    double* Dist_value_pt;
 
    int* Dist_index_pt;
 
    // required by SuperLU_DIST
    int* Dist_start_pt;
 
#endif
  }; // end of SuperLUSolver
} // namespace oomph
#endif