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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 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;

//====================================================================
/// Base class for all linear solvers. This merely defines standard 
/// interfaces for linear solvers, so that different solvers can be
/// used in a clean and transparent manner. Note that LinearSolvers
/// are primarily used to solve the linear systems arising in 
/// oomph-lib's Newton iteration. Their primary solve function
/// therefore takes a pointer to the associated problem, construct its 
/// Jacobian matrix and residual vector, and return the solution
/// of the linear system formed by the Jacobian and the residual vector.
///  We also provide broken virtual interfaces
/// to a linear-algebra-type solve function in which the matrix
/// and the rhs can be specified, but this are not guaranteed to
/// implemented for all linear solvers (e.g. for frontal solvers).
//====================================================================
class LinearSolver : public DistributableLinearAlgebraObject
{
  protected:
 
 ///\short Boolean that indicates whether the matrix (or its factors, in
 ///the case of direct solver) should be stored so that the resolve function
 ///can be used.
 bool Enable_resolve;

 /// \short Boolean flag that indicates whether the time taken
 // for the solves should be documented
 bool Doc_time;

 /// \short flag that indicates whether the gradient required for the
 /// globally convergent Newton method should be computed or not
 bool Compute_gradient;

 /// \short flag that indicates whether the gradient was computed or not
 bool Gradient_has_been_computed;

 /// \short DoubleVector storing the gradient for the globally convergent
 /// Newton method
 DoubleVector Gradient_for_glob_conv_newton_solve;

  public:
 
 /// Empty constructor, initialise the member data
  LinearSolver() : Enable_resolve(false), Doc_time(true), 
  Compute_gradient(false), Gradient_has_been_computed(false) {}
 
 /// Broken copy constructor
 LinearSolver(const LinearSolver& dummy) 
  { 
   BrokenCopy::broken_copy("LinearSolver");
  } 
 
 /// Broken assignment operator
 void operator=(const LinearSolver&) 
  {
   BrokenCopy::broken_assign("LinearSolver");
  }
 
 /// Empty virtual destructor
 virtual ~LinearSolver() {}

 /// \short Enable documentation of solve times
 void enable_doc_time() {Doc_time=true;}

 /// \short Disable documentation of solve times
 void disable_doc_time() {Doc_time=false;}

 /// \short Is documentation of solve times enabled?
 bool is_doc_time_enabled() const {return Doc_time;}

 /// \short Boolean flag indicating if resolves are enabled
 bool is_resolve_enabled() const {return Enable_resolve;}

 /// \short Enable resolve (i.e. store matrix and/or LU decomposition, say)
 /// Virtual so it can be overloaded to perform additional tasks
 virtual void enable_resolve() {Enable_resolve=true;}

 /// \short Disable resolve (i.e. store matrix and/or LU decomposition, say)
 /// This function simply resets an internal flag. It's virtual so 
 /// it can be overloaded to perform additional tasks such as 
 /// cleaning up memory that is only required for the resolve.
 virtual void disable_resolve() {Enable_resolve=false;}

 /// \short Solver: Takes pointer to problem and returns the results vector
 /// which contains the solution of the linear system defined by
 /// the problem's fully assembled Jacobian and residual vector.
 virtual void solve(Problem* const &problem_pt, DoubleVector &result)=0;

 /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs 
 /// vector and returns the solution of the linear system. 
 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);
  }

 /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs 
 /// vector and returns the solution of the linear system. 
 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);
  }
 
 /// \short Resolve the system defined by the last assembled jacobian
 /// and the rhs vector. Solution is returned in the vector result.
 /// (broken virtual)
 virtual void resolve(const DoubleVector &rhs, DoubleVector &result)
  {
   throw OomphLibError(
    "Resolve function not implemented for this linear solver",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }

 /// \short Empty virtual function that can be overloaded in specific
 /// linear solvers to clean up any memory that may have been
 /// allocated (e.g. when preparing for a re-solve).
 virtual void clean_up_memory() {}
 
 ///  \short returns the time taken to assemble the Jacobian matrix and 
 /// residual vector (needs to be overloaded for each solver)
 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;
  }

 /// \short return the time taken to solve the linear system (needs to be 
 /// overloaded for each linear solver)
 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;
 }

 /// \short function to enable the computation of the gradient required 
 /// for the globally convergent Newton method
 virtual void enable_computation_of_gradient()
 {
  throw OomphLibError(
   "enable_computation_of_gradient() not implemented for "
   "this linear solver",
   OOMPH_CURRENT_FUNCTION,
   OOMPH_EXCEPTION_LOCATION);
 }

 /// \short function to disable the computation of the gradient required 
 /// for the globally convergent Newton method
 void disable_computation_of_gradient()
 {
  Compute_gradient=false;
 }

 /// \short function to reset the size of the gradient before each Newton
 /// solve
 void reset_gradient()
 {
  Gradient_for_glob_conv_newton_solve.clear();
 }

 /// \short function to access the gradient, provided it has been computed
 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
 }

};

//=============================================================================
/// \short Dense LU decomposition-based solve of full assembled linear system.
/// VERY inefficient but useful to illustrate the principle.
/// Only suitable for use with Serial matrices and vectors.
/// This solver will only work with non-distributed matrices and vectors
/// (note: DenseDoubleMatrix is not distributable)
//============================================================================
class DenseLU : public LinearSolver
{
 ///The DenseDoubleMatrix class is a friend
 friend class DenseDoubleMatrix;
 
  public:

 /// Constructor, initialise storage
 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;
  }

 /// Broken copy constructor
 DenseLU(const DenseLU& dummy)
  { 
   BrokenCopy::broken_copy("DenseLU");
  } 
 
 /// Broken assignment operator
 void operator=(const DenseLU&) 
  {
   BrokenCopy::broken_assign("DenseLU");
  }

 /// Destructor, clean up the stored LU factors
 ~DenseLU() {clean_up_memory();}

 /// \short Solver: Takes pointer to problem and returns the results Vector
 /// which contains the solution of the linear system defined by
 /// the problem's fully assembled Jacobian and residual Vector.
 void solve(Problem* const &problem_pt, DoubleVector &result);

 /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs 
 /// vector and returns the solution of the linear system.
 void solve(DoubleMatrixBase* const &matrix_pt,const DoubleVector &rhs,
            DoubleVector &result);

 /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs 
 /// vector and returns the solution of the linear system.
 void solve(DoubleMatrixBase* const &matrix_pt,const Vector<double> &rhs,
            Vector<double> &result);

 ///  \short returns the time taken to assemble the jacobian matrix and 
 /// residual vector
 double jacobian_setup_time() const
  {
   return Jacobian_setup_time;
  }

 /// \short return the time taken to solve the linear system (needs to be 
 /// overloaded for each linear solver)
 virtual double linear_solver_solution_time() const
  {
   return Solution_time;
  }

  protected:

 /// Perform the LU decomposition of the matrix
 void factorise(DoubleMatrixBase* const &matrix_pt);

 /// Do the backsubstitution step to solve the system LU result = rhs
 void backsub(const DoubleVector &rhs,
              DoubleVector &result);

 /// perform back substitution using Vector<double>
 void backsub(const Vector<double> &rhs, Vector<double> &result);

 /// Clean up the stored LU factors
 void clean_up_memory();

 /// Jacobian setup time
 double Jacobian_setup_time;
 
 /// Solution time
 double Solution_time;

 /// \short Sign of the determinant of the matrix 
 /// (obtained during the LU decomposition)
 int Sign_of_determinant_of_matrix;
 
  private:

 /// Pointer to storage for the index of permutations in the LU solve
 long* Index;
 
 /// Pointer to storage for the LU decomposition
 double* LU_factors;
};


//====================================================================
/// \short Dense LU decomposition-based solve of linear system
/// assembled via finite differencing of the residuals Vector.
/// Even more inefficient than DenseLU but excellent sanity check!
//====================================================================
class FD_LU : public DenseLU
{
  public:

 /// Constructor: empty
 FD_LU() : DenseLU() {}

 /// Broken copy constructor
 FD_LU(const FD_LU& dummy) 
  { 
   BrokenCopy::broken_copy("FD_LU");
  } 
 
 /// Broken assignment operator
 void operator=(const FD_LU&) 
  {
   BrokenCopy::broken_assign("FD_LU");
  }

 /// \short Solver: Takes pointer to problem and returns the results Vector
 /// which contains the solution of the linear system defined by
 /// the problem's residual Vector (Jacobian computed by FD approx.)
 void solve(Problem* const &problem_pt, DoubleVector &result);

 /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs 
 /// vector and returns the solution of the linear system. 
 void solve(DoubleMatrixBase* const &matrix_pt,
            const DoubleVector &rhs,
            DoubleVector &result)
  {DenseLU::solve(matrix_pt,rhs,result);}

 /// \short 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);}

};




///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////




//=============================================================================
/// \short. SuperLU Project Solver class. This is a combined wrapper for both 
/// SuperLU and SuperLU Dist.
/// See http://crd.lbl.gov/~xiaoye/SuperLU/
/// Default Behaviour: If this solver is distributed over more than one 
/// processor then SuperLU Dist is used.
/// Member data naming convention: member data associated with the SuperLU
/// Dist solver begins Dist_... and member data associated with the serial
/// SuperLU solver begins Serial_... . 
//=============================================================================
class SuperLUSolver : public LinearSolver
{
 
  public:

 /// \short enum type to specify the solver behaviour.
 /// Default - will employ superlu dist if more than 1 processor.
 /// Serial - will always try use superlu (serial).
 /// Distributed - will always try to use superlu dist.
 enum Type { Default, Serial, Distributed };

 /// Constructor. Set the defaults.
 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;
  }

 /// Broken copy constructor
 SuperLUSolver(const SuperLUSolver& dummy) 
  { 
   BrokenCopy::broken_copy("SuperLUSolver");
  } 
 
 /// Broken assignment operator
 void operator=(const SuperLUSolver&) 
  {
   BrokenCopy::broken_assign("SuperLUSolver");
  }

 ///Destructor, clean up the stored matrices
 ~SuperLUSolver()
  {
   clean_up_memory();
  }

 /// function to enable the computation of the gradient
 void enable_computation_of_gradient()
 {
  Compute_gradient=true;
 }

 // SuperLUSolver methods
 ////////////////////////
 
 /// Overload disable resolve so that it cleans up memory too
 void disable_resolve()
  {
   LinearSolver::disable_resolve();
   clean_up_memory();
  }
 
 /// \short Solver: Takes pointer to problem and returns the results Vector
 /// which contains the solution of the linear system defined by
 /// the problem's fully assembled Jacobian and residual Vector.
 void solve(Problem* const &problem_pt, DoubleVector &result);

 /// \short Linear-algebra-type solver: Takes pointer to a matrix and rhs 
 /// vector and returns the solution of the linear system.
 /// The function returns the global result Vector.
 /// Note: if Delete_matrix_data is true the function 
 /// matrix_pt->clean_up_memory() will be used to wipe the matrix data.
 void solve(DoubleMatrixBase* const &matrix_pt,
            const DoubleVector &rhs,
            DoubleVector &result);


/*  /// \short 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);} */

 
 /// \short Resolve the system defined by the last assembled jacobian
 /// and the specified rhs vector if resolve has been enabled.
 /// Note: returns the global result Vector.
 void resolve(const DoubleVector &rhs, DoubleVector &result);

 /// Enable documentation of solver statistics
 void enable_doc_stats() {Doc_stats = true;}

 /// Disable documentation of solver statistics
 void disable_doc_stats() {Doc_stats = false;}

 /// \short returns the time taken to assemble the jacobian matrix and 
 /// residual vector
 double jacobian_setup_time() const
  {
   return Jacobian_setup_time;
  }

 /// \short return the time taken to solve the linear system (needs to be 
 /// overloaded for each linear solver)
 virtual double linear_solver_solution_time() const
  {
   return Solution_time;
  }

 /// \short Do the factorisation stage
 /// Note: if Delete_matrix_data is true the function 
 /// matrix_pt->clean_up_memory() will be used to wipe the matrix data.
 void factorise(DoubleMatrixBase* const &matrix_pt);
  
 /// \short Do the backsubstitution for SuperLU solver
 /// Note: returns the global result Vector.
 void backsub(const DoubleVector &rhs,
              DoubleVector &result);
 
 /// Clean up the memory allocated by the solver
 void clean_up_memory();

 /// \short Specify the solve type. Either default, serial or distributed.
 /// See enum SuperLU_solver_type for more details.
 void set_solver_type(const Type& t)
  {
   this->clean_up_memory();
   Solver_type = t;
  }

 // SuperLU (serial) methods
 ///////////////////////////

 /// Use the compressed row format in superlu serial
 void use_compressed_row_for_superlu_serial()
 {Serial_compressed_row_flag=true;}

 /// Use the compressed column format in superlu serial
 void use_compressed_column_for_superlu_serial()
 {Serial_compressed_row_flag=false;}

#ifdef OOMPH_HAS_MPI

 // SuperLU Dist methods
 ///////////////////////

 /// \short Set Delete_matrix_data flag. SuperLU_dist needs its own copy 
 /// of the input matrix, therefore a copy must be made if any matrix 
 /// used with this solver is to be preserved. If the input matrix can be 
 /// deleted the flag can be set to true to reduce the amount of memory 
 /// required, and the matrix data will be wiped using its clean_up_memory()
 /// function.  Default value is false.
 void enable_delete_matrix_data_in_superlu_dist()
 {Dist_delete_matrix_data=true;}
  
 /// \short Unset Delete_matrix_data flag. SuperLU_dist needs its own copy 
 /// of the input matrix, therefore a copy must be made if any matrix 
 /// used with this solver is to be preserved. If the input matrix can be 
 void disable_delete_matrix_data_in_superlu_dist()
 {Dist_delete_matrix_data=false;}

 /// \short Set flag so that SuperLU_DIST is allowed to permute matrix rows
 /// and columns during factorisation. This is the default for SuperLU_DIST, 
 /// and can lead to significantly faster solves.
 void enable_row_and_col_permutations_in_superlu_dist()
 {Dist_allow_row_and_col_permutations=true;}

 /// \short Set flag so that SuperLU_DIST is not allowed to permute matrix rows
 /// and columns during factorisation. 
 void disable_row_and_col_permutations_in_superlu_dist()
 {Dist_allow_row_and_col_permutations=false;}

  /// \short Calling this method will ensure that when the problem based
  /// solve interface is used, a global (serial) jacobian will be 
  /// assembled.
  /// Note: calling this function will delete any distributed solve data.
  void use_global_solve_in_superlu_dist() 
  {
   if (!Dist_use_global_solver)
    {
     clean_up_memory();
     Dist_use_global_solver=true;
    }
  }

  /// \short Calling this method will ensure that when the problem based
  /// solve interface is used, a distributed jacobian will be 
  /// assembled.
  /// Note: calling this function will delete any global solve data.
  void use_distributed_solve_in_superlu_dist() 
  {
   if (Dist_use_global_solver)
    {
     clean_up_memory();
     Dist_use_global_solver=false;
    }
  }

#endif

  private:

  /// factorise method for SuperLU (serial)
 void factorise_serial(DoubleMatrixBase* const &matrix_pt);
  
 /// backsub method for SuperLU (serial)
 void backsub_serial(const DoubleVector &rhs,
                     DoubleVector &result);

#ifdef OOMPH_HAS_MPI
  /// factorise method for SuperLU Dist
 void factorise_distributed(DoubleMatrixBase* const &matrix_pt);
  
 /// backsub method for SuperLU Dist
 void backsub_distributed(const DoubleVector &rhs,
              DoubleVector &result);
#endif

 // SuperLUSolver member data
 ////////////////////////////

 /// Jacobian setup time
 double Jacobian_setup_time;

 /// Solution time
 double Solution_time;

 /// Suppress solve?
 bool Suppress_solve;

 /// Set to true to output statistics (false by default).
 bool Doc_stats;

 /// the solver type. see SuperLU_solver_type for details.
 Type Solver_type;

 /// boolean flag indicating whether superlu dist is being used
 bool Using_dist;

 // SuperLU (serial) member data
 ///////////////////////////////

 /// Storage for the LU factors as required by SuperLU
 void *Serial_f_factors;

 /// Info flag for the SuperLU solver
 int Serial_info;

 /// The number of unknowns in the linear system
 unsigned long Serial_n_dof;
 
 /// Sign of the determinant of the matrix
 int Serial_sign_of_determinant_of_matrix;
 
 /// Use compressed row version?
 bool Serial_compressed_row_flag;

#ifdef OOMPH_HAS_MPI

 // SuperLU Dist member data
 ///////////////////////////
public:
 /// \short Static flag that determines whether the warning about
 /// incorrect distribution of RHSs will be printed or not
 static bool Suppress_incorrect_rhs_distribution_warning_in_resolve;

private:
 /// \short Flag that determines whether the MPIProblem based solve function 
 /// uses the global or distributed version of SuperLU_DIST 
 /// (default value is false).
 bool Dist_use_global_solver;
 
 ///  Flag is true if solve data has been generated for a global matrix
 bool Dist_global_solve_data_allocated;

 ///  Flag is true if solve data has been generated for distributed matrix
 bool Dist_distributed_solve_data_allocated;

 /// Storage for the LU factors and other data required by SuperLU
 void *Dist_solver_data_pt;

 /// Number of rows for the process grid 
 int Dist_nprow;

 /// Number of columns for the process grid
 int Dist_npcol;

 /// Info flag for the SuperLU solver
 int Dist_info;
 
 /// \short If true then SuperLU_DIST is allowed to permute matrix rows
 /// and columns during factorisation. This is the default for SuperLU_DIST, 
 /// and can lead to significantly faster solves, but has been known to 
 /// fail, hence the default value is 0.
 bool Dist_allow_row_and_col_permutations;
 
 /// \short Delete_matrix_data flag. SuperLU_dist needs its own copy 
 /// of the input matrix, therefore a copy must be made if any matrix 
 /// used with this solver is to be preserved. If the input matrix can be 
 /// deleted the flag can be set to true to reduce the amount of memory 
 /// required, and the matrix data will be wiped using its clean_up_memory()
 /// function. Default value is false.
 bool Dist_delete_matrix_data;

 /// Pointer for storage of the matrix values required by SuperLU_DIST
 double *Dist_value_pt;

 /// \short Pointer for storage of matrix rows or column indices required 
 /// by SuperLU_DIST
 int *Dist_index_pt;

 /// \short Pointers for storage of matrix column or row starts 
 // required by SuperLU_DIST
 int *Dist_start_pt;

#endif
}; // end of SuperLUSolver 
} // end of oomph namespace
#endif