oomph::SuperLUSolver Class Reference

#include <linear_solver.h>

Inheritance diagram for oomph::SuperLUSolver:

Public Types
enum	Type { Default , Serial , Distributed }

Public Member Functions
	SuperLUSolver ()
	Constructor. Set the defaults. More...
	SuperLUSolver (const SuperLUSolver &dummy)=delete
	Broken copy constructor. More...
void	operator= (const SuperLUSolver &)=delete
	Broken assignment operator. More...
	~SuperLUSolver ()
	Destructor, clean up the stored matrices. More...
void	enable_computation_of_gradient ()
	function to enable the computation of the gradient More...
void	disable_resolve ()
	Overload disable resolve so that it cleans up memory too. More...
void	solve (Problem *const &problem_pt, DoubleVector &result)
void	solve (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &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)
	Resolve the system for a given RHS. More...
void	resolve_transpose (const DoubleVector &rhs, DoubleVector &result)
	Resolve the (transposed) system for a given RHS. More...
void	enable_doc_stats ()
	Enable documentation of solver statistics. More...
void	disable_doc_stats ()
	Disable documentation of solver statistics. More...
double	jacobian_setup_time () const
virtual double	linear_solver_solution_time () const
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 ()
	Clean up the memory allocated by the solver. More...
void	set_solver_type (const Type &t)
void	use_compressed_row_for_superlu_serial ()
	Use the compressed row format in superlu serial. More...
void	use_compressed_column_for_superlu_serial ()
	Use the compressed column format in superlu serial. More...
double	get_memory_usage_for_lu_factors ()
	How much memory do the LU factors take up? In bytes. More...
double	get_total_needed_memory ()
	How much memory was allocated by SuperLU? In bytes. More...
Public Member Functions inherited from oomph::LinearSolver
	LinearSolver ()
	Empty constructor, initialise the member data. More...
	LinearSolver (const LinearSolver &dummy)=delete
	Broken copy constructor. More...
void	operator= (const LinearSolver &)=delete
	Broken assignment operator. More...
virtual	~LinearSolver ()
	Empty virtual destructor. More...
void	enable_doc_time ()
	Enable documentation of solve times. More...
void	disable_doc_time ()
	Disable documentation of solve times. More...
bool	is_doc_time_enabled () const
	Is documentation of solve times enabled? More...
bool	is_resolve_enabled () const
	Boolean flag indicating if resolves are enabled. More...
virtual void	enable_resolve ()
virtual void	solve (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
virtual void	solve_transpose (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
void	disable_computation_of_gradient ()
void	reset_gradient ()
void	get_gradient (DoubleVector &gradient)
	function to access the gradient, provided it has been computed More...
Public Member Functions inherited from oomph::DistributableLinearAlgebraObject
	DistributableLinearAlgebraObject ()
	Default constructor - create a distribution. More...
	DistributableLinearAlgebraObject (const DistributableLinearAlgebraObject &matrix)=delete
	Broken copy constructor. More...
void	operator= (const DistributableLinearAlgebraObject &)=delete
	Broken assignment operator. More...
virtual	~DistributableLinearAlgebraObject ()
	Destructor. More...
LinearAlgebraDistribution *	distribution_pt () const
	access to the LinearAlgebraDistribution More...
unsigned	nrow () const
	access function to the number of global rows. More...
unsigned	nrow_local () const
	access function for the num of local rows on this processor. More...
unsigned	nrow_local (const unsigned &p) const
	access function for the num of local rows on this processor. More...
unsigned	first_row () const
	access function for the first row on this processor More...
unsigned	first_row (const unsigned &p) const
	access function for the first row on this processor More...
bool	distributed () const
	distribution is serial or distributed More...
bool	distribution_built () const
void	build_distribution (const LinearAlgebraDistribution *const dist_pt)
void	build_distribution (const LinearAlgebraDistribution &dist)

Private Member Functions
void	factorise_serial (DoubleMatrixBase *const &matrix_pt)
	factorise method for SuperLU (serial) More...
void	backsub_serial (const DoubleVector &rhs, DoubleVector &result)
	backsub method for SuperLU (serial) More...
void	backsub_transpose_serial (const DoubleVector &rhs, DoubleVector &result)
	backsub method for SuperLU (serial) More...

Private Attributes
double	Jacobian_setup_time
	Jacobian setup time. More...
double	Solution_time
	Solution time. More...
bool	Suppress_solve
	Suppress solve? More...
bool	Doc_stats
	Set to true to output statistics (false by default). More...
Type	Solver_type
	the solver type. see SuperLU_solver_type for details. More...
bool	Using_dist
	boolean flag indicating whether superlu dist is being used More...
void *	Serial_f_factors
	Storage for the LU factors as required by SuperLU. More...
int	Serial_info
	Info flag for the SuperLU solver. More...
unsigned long	Serial_n_dof
	The number of unknowns in the linear system. More...
int	Serial_sign_of_determinant_of_matrix
	Sign of the determinant of the matrix. More...
bool	Serial_compressed_row_flag
	Use compressed row version? More...

Additional Inherited Members
Protected Member Functions inherited from oomph::DistributableLinearAlgebraObject
void	clear_distribution ()
Protected Attributes inherited from oomph::LinearSolver
bool	Enable_resolve
bool	Doc_time
	Boolean flag that indicates whether the time taken. More...
bool	Compute_gradient
bool	Gradient_has_been_computed
	flag that indicates whether the gradient was computed or not More...
DoubleVector	Gradient_for_glob_conv_newton_solve

Detailed Description

//////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// 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_... .

Member Enumeration Documentation

Type

enum oomph::SuperLUSolver::Type

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.

Enumerator
Default
Serial
Distributed

    {
      Default,
      Serial,
      Distributed
    };

Constructor & Destructor Documentation

SuperLUSolver() [1/2]

oomph::SuperLUSolver::SuperLUSolver ( )

inline

Constructor. Set the defaults.

                    : 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;
    }

References Default, Doc_stats, Serial_compressed_row_flag, Serial_n_dof, Serial_sign_of_determinant_of_matrix, Solver_type, Suppress_solve, and Using_dist.

SuperLUSolver() [2/2]

oomph::SuperLUSolver::SuperLUSolver ( const SuperLUSolver &  dummy )

delete

Broken copy constructor.

~SuperLUSolver()

oomph::SuperLUSolver::~SuperLUSolver ( )

inline

Destructor, clean up the stored matrices.

    {
      clean_up_memory();
    }

References clean_up_memory().

Member Function Documentation

backsub()

void oomph::SuperLUSolver::backsub	(	const DoubleVector &	rhs,
		DoubleVector &	result
	)

Do the backsubstitution for SuperLU solver Note: returns the global result Vector.

Do the backsubstitution for SuperLUSolver. Note - this method performs no paranoid checks - these are all performed in solve(...) and resolve(...)

  {
#ifdef OOMPH_HAS_MPI
    if (Using_dist)
    {
      backsub_distributed(rhs, result);
    }
    else
#endif
    {
      backsub_serial(rhs, result);
    }
  }

References backsub_serial(), and Using_dist.

Referenced by resolve(), and solve().

backsub_serial()

void oomph::SuperLUSolver::backsub_serial ( const DoubleVector &  rhs, DoubleVector &  result  )

private

backsub method for SuperLU (serial)

Do the backsubstitution for SuperLU.

  {
    // Find the number of unknowns
    int n = rhs.nrow();
 
#ifdef PARANOID
    // PARANOID check that this rhs distribution is setup
    if (!rhs.built())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector distribution must be setup.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // PARANOID check that the rhs has the right number of global rows
    if (static_cast<int>(Serial_n_dof) != n)
    {
      throw OomphLibError(
        "RHS does not have the same dimension as the linear system",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
    // PARANOID check that the rhs is not distributed
    if (rhs.distribution_pt()->distributed())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector must not be distributed.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // PARANOID check that if the result is setup it matches the distribution
    // of the rhs
    if (result.built())
    {
      if (!(*rhs.distribution_pt() == *result.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream << "If the result distribution is setup then it "
                                "must be the same as the "
                             << "rhs distribution";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // copy result to rhs
    result = rhs;
 
    // Number of RHSs
    int nrhs = 1;
 
    // Cast the boolean flags to ints for SuperLU
    int transpose = Serial_compressed_row_flag;
    int doc = Doc_stats;
 
    // Do the backsubsitition phase
    int i = 2;
    superlu(&i,
            &n,
            0,
            &nrhs,
            0,
            0,
            0,
            result.values_pt(),
            &n,
            &transpose,
            &doc,
            &Serial_f_factors,
            &Serial_info);
 
    // Throw an error if superLU returned an error status in info.
    if (Serial_info != 0)
    {
      std::ostringstream error_msg;
      error_msg << "SuperLU returned the error status code " << Serial_info
                << " . See the SuperLU documentation for what this means.";
      throw OomphLibError(
        error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }

References oomph::DoubleVector::built(), oomph::LinearAlgebraDistribution::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), Doc_stats, i, n, oomph::DistributableLinearAlgebraObject::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Serial_compressed_row_flag, Serial_f_factors, Serial_info, Serial_n_dof, oomph::superlu(), anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose(), and oomph::DoubleVector::values_pt().

Referenced by backsub().

backsub_transpose()

void oomph::SuperLUSolver::backsub_transpose	(	const DoubleVector &	rhs,
		DoubleVector &	result
	)

Do the back-substitution for transposed system of the SuperLU solver Note: Returns the global result Vector.

Do the backsubstitution of the transposed system for SuperLUSolver. Note - this method performs no paranoid checks - these are all performed in solve(...) and resolve(...)

  {
#ifdef OOMPH_HAS_MPI
    if (Using_dist)
    {
      backsub_transpose_distributed(rhs, result);
    }
    else
#endif
    {
      backsub_transpose_serial(rhs, result);
    }
  }

References backsub_transpose_serial(), and Using_dist.

Referenced by resolve_transpose(), and solve_transpose().

backsub_transpose_serial()

void oomph::SuperLUSolver::backsub_transpose_serial ( const DoubleVector &  rhs, DoubleVector &  result  )

private

backsub method for SuperLU (serial)

Do the backsubstitution for SuperLU.

  {
    // Find the number of unknowns
    int n = rhs.nrow();
 
#ifdef PARANOID
    // PARANOID check that this rhs distribution is setup
    if (!rhs.built())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector distribution must be setup.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // PARANOID check that the rhs has the right number of global rows
    if (static_cast<int>(Serial_n_dof) != n)
    {
      throw OomphLibError(
        "RHS does not have the same dimension as the linear system",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
    // PARANOID check that the rhs is not distributed
    if (rhs.distribution_pt()->distributed())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The rhs vector must not be distributed.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
    // PARANOID check that if the result is setup it matches the distribution
    // of the rhs
    if (result.built())
    {
      if (!(*rhs.distribution_pt() == *result.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream << "If the result distribution is setup then it "
                                "must be the same as the "
                             << "rhs distribution";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // copy result to rhs
    result = rhs;
 
    // Number of RHSs
    int nrhs = 1;
 
    // Cast the boolean flags to ints for SuperLU
    int transpose = (!Serial_compressed_row_flag);
    int doc = Doc_stats;
 
    // Do the backsubsitition phase
    int i = 2;
    superlu(&i,
            &n,
            0,
            &nrhs,
            0,
            0,
            0,
            result.values_pt(),
            &n,
            &transpose,
            &doc,
            &Serial_f_factors,
            &Serial_info);
 
    // Throw an error if superLU returned an error status in info.
    if (Serial_info != 0)
    {
      std::ostringstream error_msg;
      error_msg << "SuperLU returned the error status code " << Serial_info
                << " . See the SuperLU documentation for what this means.";
      throw OomphLibError(
        error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }

References oomph::DoubleVector::built(), oomph::LinearAlgebraDistribution::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), Doc_stats, i, n, oomph::DistributableLinearAlgebraObject::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Serial_compressed_row_flag, Serial_f_factors, Serial_info, Serial_n_dof, oomph::superlu(), anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose(), and oomph::DoubleVector::values_pt().

Referenced by backsub_transpose().

clean_up_memory()

void oomph::SuperLUSolver::clean_up_memory ( )

virtual

Clean up the memory allocated by the solver.

Clean up the memory.

Reimplemented from oomph::LinearSolver.

  {
    // If we have non-zero LU factors stored
    if (Serial_f_factors != 0)
    {
      // Clean up those factors
      int i = 3;
      int transpose = Serial_compressed_row_flag;
      superlu(&i,
              0,
              0,
              0,
              0,
              0,
              0,
              0,
              0,
              &transpose,
              0,
              &Serial_f_factors,
              &Serial_info);
 
      // Set the F_factors to zero
      Serial_f_factors = 0;
      Serial_n_dof = 0;
    }
 
#ifdef OOMPH_HAS_MPI
    // If we have non-zero LU factors stored
    if (Dist_solver_data_pt != 0)
    {
      // Clean up any stored solver data
 
      // Doc (0/1) = (true/false)
      int doc = !Doc_stats;
 
      // Reset Info flag
      Dist_info = 0;
 
      // number of DOFs
      int ndof = this->distribution_pt()->nrow();
 
      if (Dist_distributed_solve_data_allocated)
      {
        superlu_dist_distributed_matrix(
          3,
          -1,
          ndof,
          0,
          0,
          0,
          0,
          0,
          0,
          0,
          Dist_nprow,
          Dist_npcol,
          doc,
          &Dist_solver_data_pt,
          &Dist_info,
          this->distribution_pt()->communicator_pt()->mpi_comm());
        Dist_distributed_solve_data_allocated = false;
      }
      if (Dist_global_solve_data_allocated)
      {
        superlu_dist_global_matrix(
          3,
          -1,
          ndof,
          0,
          0,
          0,
          0,
          0,
          Dist_nprow,
          Dist_npcol,
          doc,
          &Dist_solver_data_pt,
          &Dist_info,
          this->distribution_pt()->communicator_pt()->mpi_comm());
        Dist_global_solve_data_allocated = false;
      }
 
      Dist_solver_data_pt = 0;
 
      // Delete internal copy of the matrix
      delete[] Dist_value_pt;
      delete[] Dist_index_pt;
      delete[] Dist_start_pt;
      Dist_value_pt = 0;
      Dist_index_pt = 0;
      Dist_start_pt = 0;
 
      // and the distribution
      this->clear_distribution();
    }
#endif
  }

References oomph::DistributableLinearAlgebraObject::clear_distribution(), oomph::DistributableLinearAlgebraObject::distribution_pt(), Doc_stats, i, oomph::LinearAlgebraDistribution::nrow(), Serial_compressed_row_flag, Serial_f_factors, Serial_info, Serial_n_dof, oomph::superlu(), and anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose().

Referenced by oomph::SuperLUPreconditioner::clean_up_memory(), disable_resolve(), factorise(), factorise_serial(), set_solver_type(), solve(), solve_transpose(), and ~SuperLUSolver().

disable_doc_stats()

void oomph::SuperLUSolver::disable_doc_stats ( )

inline

Disable documentation of solver statistics.

    {
      Doc_stats = false;
    }

References Doc_stats.

Referenced by oomph::SuperLUPreconditioner::disable_doc_stats(), and oomph::SuperLUPreconditioner::SuperLUPreconditioner().

disable_resolve()

void oomph::SuperLUSolver::disable_resolve ( )

inlinevirtual

Overload disable resolve so that it cleans up memory too.

Reimplemented from oomph::LinearSolver.

    {
      LinearSolver::disable_resolve();
      clean_up_memory();
    }

References clean_up_memory(), and oomph::LinearSolver::disable_resolve().

enable_computation_of_gradient()

void oomph::SuperLUSolver::enable_computation_of_gradient ( )

inlinevirtual

function to enable the computation of the gradient

Reimplemented from oomph::LinearSolver.

    {
      Compute_gradient = true;
    }

References oomph::LinearSolver::Compute_gradient.

enable_doc_stats()

void oomph::SuperLUSolver::enable_doc_stats ( )

inline

Enable documentation of solver statistics.

    {
      Doc_stats = true;
    }

References Doc_stats.

Referenced by oomph::SuperLUPreconditioner::enable_doc_stats().

factorise()

void oomph::SuperLUSolver::factorise

(

DoubleMatrixBase *const &

matrix_pt

)

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.

LU decompose the matrix addressed by matrix_pt by using the SuperLU solver. The resulting matrix factors are stored internally.

  {
    // wipe memory
    this->clean_up_memory();
 
    // if we have mpi and the solver is distributed or default and nproc
    // gt 1
#ifdef OOMPH_HAS_MPI
    DistributableLinearAlgebraObject* dist_matrix_pt =
      dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
    unsigned nproc = 1;
    if (dist_matrix_pt != 0)
    {
      nproc = dist_matrix_pt->distribution_pt()->communicator_pt()->nproc();
    }
    if (Solver_type == Distributed || (Solver_type == Default && nproc > 1 &&
                                       MPI_Helpers::mpi_has_been_initialised()))
    {
      // if the matrix is a distributed linear algebra object then use SuperLU
      // dist
      if (dist_matrix_pt != 0)
      {
        factorise_distributed(matrix_pt);
        Using_dist = true;
      }
      else
      {
        factorise_serial(matrix_pt);
        Using_dist = false;
      }
    }
    else
#endif
    {
      factorise_serial(matrix_pt);
      Using_dist = false;
    }
  }

References clean_up_memory(), oomph::LinearAlgebraDistribution::communicator_pt(), Default, Distributed, oomph::DistributableLinearAlgebraObject::distribution_pt(), factorise_serial(), oomph::MPI_Helpers::mpi_has_been_initialised(), oomph::OomphCommunicator::nproc(), Solver_type, and Using_dist.

Referenced by oomph::SuperLUPreconditioner::setup(), solve(), and solve_transpose().

factorise_serial()

void oomph::SuperLUSolver::factorise_serial ( DoubleMatrixBase *const &  matrix_pt )

private

factorise method for SuperLU (serial)

LU decompose the matrix addressed by matrix_pt by using the SuperLU solver. The resulting matrix factors are stored internally.

  {
#ifdef PARANOID
    // PARANOID check that if the matrix is distributable then it should not be
    // then it should not be distributed
    if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt) != 0)
    {
      if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
            ->distributed())
      {
        std::ostringstream error_message_stream;
        error_message_stream << "The matrix must not be distributed.";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // Find # of degrees of freedom (variables)
    int n = matrix_pt->nrow();
 
    // Check that we have a square matrix
#ifdef PARANOID
    int m = matrix_pt->ncol();
    if (n != m)
    {
      std::ostringstream error_message_stream;
      error_message_stream << "Can only solve for square matrices\n"
                           << "N, M " << n << " " << m << std::endl;
 
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Storage for the values, rows and column indices
    // required by SuplerLU
    double* value = 0;
    int *index = 0, *start = 0;
 
    // Integer used to represent compressed row or column format
    // Default compressed row
    int transpose = 0;
 
    // Number of non-zero entries in the matrix
    int nnz = 0;
 
    // Doc flag (convert to int for SuperLU)
    int doc = Doc_stats;
 
    // Is it a CR matrix
    if (dynamic_cast<CRDoubleMatrix*>(matrix_pt))
    {
      // Set the appropriate row flags
      Serial_compressed_row_flag = true;
      transpose = 1;
      // Get a cast pointer to the matrix
      CRDoubleMatrix* CR_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
 
      // Now set the pointers to the interanally stored values
      // and indices
      nnz = CR_matrix_pt->nnz();
      value = CR_matrix_pt->value();
      index = CR_matrix_pt->column_index();
      start = CR_matrix_pt->row_start();
    }
    // Otherwise is it the compressed column version?
    else if (dynamic_cast<CCDoubleMatrix*>(matrix_pt))
    {
      // Set the compressed row flag to false
      Serial_compressed_row_flag = false;
      // Get a cast pointer to the matrix
      CCDoubleMatrix* CC_matrix_pt = dynamic_cast<CCDoubleMatrix*>(matrix_pt);
 
      // Now set the pointers to the interanally stored values
      // and indices
      nnz = CC_matrix_pt->nnz();
      value = CC_matrix_pt->value();
      index = CC_matrix_pt->row_index();
      start = CC_matrix_pt->column_start();
    }
    // Otherwise throw and error
    else
    {
      throw OomphLibError("SuperLU only works with CR or CC Double matrices",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // Clean up any previous storage so that if this is called twice with
    // the same matrix, we don't get a memory leak
    clean_up_memory();
 
    // Perform the lu decompose phase (i=1)
    int i = 1;
    Serial_sign_of_determinant_of_matrix = superlu(&i,
                                                   &n,
                                                   &nnz,
                                                   0,
                                                   value,
                                                   index,
                                                   start,
                                                   0,
                                                   &n,
                                                   &transpose,
                                                   &doc,
                                                   &Serial_f_factors,
                                                   &Serial_info);
 
    // Throw an error if superLU returned an error status in info.
    if (Serial_info != 0)
    {
      std::ostringstream error_msg;
      error_msg << "SuperLU returned the error status code " << Serial_info
                << " . See the SuperLU documentation for what this means.";
      throw OomphLibError(
        error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
 
    // Set the number of degrees of freedom in the linear system
    Serial_n_dof = n;
  }

References clean_up_memory(), oomph::CRDoubleMatrix::column_index(), oomph::CCMatrix< T >::column_start(), oomph::DistributableLinearAlgebraObject::distributed(), Doc_stats, i, m, n, oomph::DoubleMatrixBase::ncol(), oomph::SparseMatrix< T, MATRIX_TYPE >::nnz(), oomph::CRDoubleMatrix::nnz(), oomph::DoubleMatrixBase::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::CCMatrix< T >::row_index(), oomph::CRDoubleMatrix::row_start(), Serial_compressed_row_flag, Serial_f_factors, Serial_info, Serial_n_dof, Serial_sign_of_determinant_of_matrix, oomph::CumulativeTimings::start(), oomph::superlu(), anonymous_namespace{skew_symmetric_matrix3.cpp}::transpose(), Eigen::value, oomph::SparseMatrix< T, MATRIX_TYPE >::value(), and oomph::CRDoubleMatrix::value().

Referenced by factorise().

get_memory_usage_for_lu_factors()

double oomph::SuperLUSolver::get_memory_usage_for_lu_factors

(

)

How much memory do the LU factors take up? In bytes.

How much memory do the LU factors take up? In bytes NOTE: This has been scraped from dQuerySpace(...) in dmemory.c in external_src/oomph_superlu_4.3

  {
    // If we're using the non-distributed version of SuperLU and the LU
    // factors have also been computed
    if ((Solver_type != Distributed) && (Serial_f_factors != 0))
    {
      return get_lu_factor_memory_usage_in_bytes();
    }
#ifdef OOMPH_HAS_MPI
    // If we're using SuperLU dist and the LU factors have been computed
    if ((Solver_type == Distributed) && (Dist_solver_data_pt != 0))
    {
      return get_lu_factor_memory_usage_in_bytes_dist();
    }
#endif
    // If the factors haven't been computed we can't do anything
    else
    {
      return 0.0;
    }
  } // End of get_memory_usage_for_lu_factors

References Distributed, oomph::get_lu_factor_memory_usage_in_bytes(), Serial_f_factors, and Solver_type.

Referenced by oomph::SuperLUPreconditioner::get_memory_usage_for_lu_factors(), and oomph::SuperLUPreconditioner::get_memory_usage_for_superlu().

get_total_needed_memory()

double oomph::SuperLUSolver::get_total_needed_memory

(

)

How much memory was allocated by SuperLU? In bytes.

How much memory was used in total? In bytes NOTE: This has been scraped from dQuerySpace(...) in dmemory.c in external_src/oomph_superlu_4.3

  {
    // If we're using the non-distributed version of SuperLU and the LU
    // factors have also been computed
    if ((Solver_type != Distributed) && (Serial_f_factors != 0))
    {
      return get_total_memory_usage_in_bytes();
    }
#ifdef OOMPH_HAS_MPI
    // If we're using SuperLU dist and the LU factors have been computed
    if ((Solver_type == Distributed) && (Dist_solver_data_pt != 0))
    {
      return get_total_memory_usage_in_bytes_dist();
    }
#endif
    // If the factors haven't been computed we can't do anything
    else
    {
      return 0.0;
    }
  } // End of get_total_needed_memory

References Distributed, oomph::get_total_memory_usage_in_bytes(), Serial_f_factors, and Solver_type.

Referenced by oomph::SuperLUPreconditioner::get_memory_usage_for_superlu(), oomph::SuperLUPreconditioner::get_total_memory_needed_for_superlu(), solve(), and solve_transpose().

jacobian_setup_time()

double oomph::SuperLUSolver::jacobian_setup_time ( ) const

inlinevirtual

returns the time taken to assemble the jacobian matrix and residual vector

Reimplemented from oomph::LinearSolver.

    {
      return Jacobian_setup_time;
    }

References Jacobian_setup_time.

linear_solver_solution_time()

virtual double oomph::SuperLUSolver::linear_solver_solution_time ( ) const

inlinevirtual

return the time taken to solve the linear system (needs to be overloaded for each linear solver)

Reimplemented from oomph::LinearSolver.

    {
      return Solution_time;
    }

References Solution_time.

operator=()

void oomph::SuperLUSolver::operator= ( const SuperLUSolver &  )

delete

Broken assignment operator.

resolve()

void oomph::SuperLUSolver::resolve ( const DoubleVector &  rhs, DoubleVector &  result  )

virtual

Resolve the system for a given RHS.

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.

Reimplemented from oomph::LinearSolver.

  {
    // Store starting time for solve
    double t_start = TimingHelpers::timer();
 
    // backsub
    backsub(rhs, result);
 
    // Doc time for solve
    double t_end = TimingHelpers::timer();
    Solution_time = t_end - t_start;
    if (Doc_time)
    {
      oomph_info << "Time for SuperLUSolver solve (ndof=" << rhs.nrow() << "): "
                 << TimingHelpers::convert_secs_to_formatted_string(t_end -
                                                                    t_start)
                 << std::endl;
    }
  }

References backsub(), oomph::TimingHelpers::convert_secs_to_formatted_string(), oomph::LinearSolver::Doc_time, oomph::DistributableLinearAlgebraObject::nrow(), oomph::oomph_info, Solution_time, and oomph::TimingHelpers::timer().

Referenced by oomph::SuperLUPreconditioner::preconditioner_solve().

resolve_transpose()

void oomph::SuperLUSolver::resolve_transpose ( const DoubleVector &  rhs, DoubleVector &  result  )

virtual

Resolve the (transposed) system for a given RHS.

Resolve the (transposed) system defined by the last assembled Jacobian and the specified rhs vector if resolve has been enabled.

Reimplemented from oomph::LinearSolver.

  {
    // Store starting time for solve
    double t_start = TimingHelpers::timer();
 
    // Backsub (but solve the transposed system)
    backsub_transpose(rhs, result);
 
    // Doc time for solve
    double t_end = TimingHelpers::timer();
    Solution_time = t_end - t_start;
    if (Doc_time)
    {
      oomph_info << "Time for SuperLUSolver solve (ndof=" << rhs.nrow() << "): "
                 << TimingHelpers::convert_secs_to_formatted_string(t_end -
                                                                    t_start)
                 << std::endl;
    }
  }

References backsub_transpose(), oomph::TimingHelpers::convert_secs_to_formatted_string(), oomph::LinearSolver::Doc_time, oomph::DistributableLinearAlgebraObject::nrow(), oomph::oomph_info, Solution_time, and oomph::TimingHelpers::timer().

Referenced by oomph::SuperLUPreconditioner::preconditioner_solve_transpose().

set_solver_type()

void oomph::SuperLUSolver::set_solver_type ( const Type &  t )

inline

Specify the solve type. Either default, serial or distributed. See enum SuperLU_solver_type for more details.

    {
      this->clean_up_memory();
      Solver_type = t;
    }

References clean_up_memory(), Solver_type, and plotPSD::t.

solve() [1/2]

void oomph::SuperLUSolver::solve ( DoubleMatrixBase *const &  matrix_pt, const DoubleVector &  rhs, DoubleVector &  result  )

virtual

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.

Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the linear system. Problem pointer defaults to NULL and can be omitted. 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.

Reimplemented from oomph::LinearSolver.

  {
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // Pointer used in various places
    CRDoubleMatrix* cr_pt = 0;
 
 
#ifdef PARANOID
    // check that the rhs vector is setup
    if (!rhs.built())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The vectors rhs must be setup";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // check that the matrix is square
    if (matrix_pt->nrow() != matrix_pt->ncol())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The matrix at matrix_pt must be square.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // check that the matrix has some entries, and so has a values_pt that
    // makes sense (only for CR because CC is never used I think dense
    // matrices will be safe since they don't use a values pointer).
    cr_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
    if (cr_pt != 0)
    {
      if (cr_pt->nnz() == 0)
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "Attempted to call SuperLu on a CRDoubleMatrix with no entries, "
          << "SuperLU would segfault (because the values array pt is "
          << "uninitialised or null).";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    // check that the matrix and the rhs vector have the same nrow()
    if (matrix_pt->nrow() != rhs.nrow())
    {
      std::ostringstream error_message_stream;
      error_message_stream
        << "The matrix and the rhs vector must have the same number of rows.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // if the matrix is distributable then should have the same distribution
    // as the rhs vector
    DistributableLinearAlgebraObject* dist_matrix_pt =
      dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
    if (dist_matrix_pt != 0)
    {
      if (!(*dist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The matrix matrix_pt must have the same distribution as the "
          << "rhs vector.";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
    // if the matrix is not distributable then it the rhs vector should not be
    // distributed
    else
    {
      if (rhs.distribution_pt()->distributed())
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The matrix (matrix_pt) is not distributable and therefore the rhs"
          << " vector must not be distributed";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
    // if the result vector is setup then check it has the same distribution
    // as the rhs
    if (result.built())
    {
      if (!(*result.distribution_pt() == *rhs.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The result vector distribution has been setup; it must have the "
          << "same distribution as the rhs vector.";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // set the distribution
    if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt))
    {
      // the solver has the same distribution as the matrix if possible
      this->build_distribution(
        dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
          ->distribution_pt());
    }
    else
    {
      // the solver has the same distribution as the RHS
      this->build_distribution(rhs.distribution_pt());
    }
 
    // Doc time for solve
    double t_factorise_start = TimingHelpers::timer();
 
    // Factorise the matrix
    factorise(matrix_pt);
 
    // Doc the end time
    double t_factorise_end = TimingHelpers::timer();
 
    // How long did the factorisation take?
    double factorise_time = t_factorise_end - t_factorise_start;
 
    // Try and upcast the matrix to a CRDoubleMatrix
    // CRDoubleMatrix*
    cr_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
 
    // If the input matrix is a CRDoubleMatrix
    if (cr_pt != 0)
    {
      // ...and actually has an entry
      if (cr_pt->nnz() != 0)
      {
        // Find out how many rows there are in the global Jacobian
        unsigned n_row = cr_pt->nrow();
 
        // And how many non-zeros there are in the global Jacobian
        unsigned n_nnz = cr_pt->nnz();
 
        // Get the memory usage (in bytes) for the global Jacobian storage
        double memory_usage_for_jacobian =
          ((2 * ((n_row + 1) * sizeof(int))) +
           (n_nnz * (sizeof(int) + sizeof(double))));
 
        // Get the memory usage (in bytes) for storage of the LU factors in
        // SuperLU
        double memory_usage_for_lu_storage = get_total_needed_memory();
 
        // Get the memory usage (in bytes) for storage of the LU factors in
        // SuperLU
        double total_memory_usage =
          memory_usage_for_jacobian + memory_usage_for_lu_storage;
 
        // How much memory have we used in the subsidiary preconditioners?
        oomph_info << "\nMemory statistics:"
                   << "\n - Memory used to store the Jacobian (MB)   : "
                   << memory_usage_for_jacobian / 1.0e+06
                   << "\n - Memory used to store the LU factors (MB) : "
                   << memory_usage_for_lu_storage / 1.0e+06
                   << "\n - Total memory used for matrix storage (MB): "
                   << total_memory_usage / 1.0e+06 << "\n"
                   << std::endl;
      }
    } // if (cr_pt!=0)
 
    // Doc the start time
    double t_backsub_start = TimingHelpers::timer();
 
    // Now do the back solve
    backsub(rhs, result);
 
    // Doc the end time
    double t_backsub_end = TimingHelpers::timer();
 
    // How long did the back substitution take?
    double backsub_time = t_backsub_end - t_backsub_start;
 
    // Doc time for solve
    double t_end = TimingHelpers::timer();
    Solution_time = t_end - t_start;
    if (Doc_time)
    {
      oomph_info
        << "Time for LU factorisation : "
        << TimingHelpers::convert_secs_to_formatted_string(factorise_time)
        << "\nTime for back-substitution: "
        << TimingHelpers::convert_secs_to_formatted_string(backsub_time)
        << "\nTime for SuperLUSolver solve (ndof=" << matrix_pt->nrow() << "): "
        << TimingHelpers::convert_secs_to_formatted_string(Solution_time)
        << std::endl;
    }
 
    // If we are not storing the solver data for resolves, delete it
    if (!Enable_resolve)
    {
      clean_up_memory();
    }
  }

References backsub(), oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DoubleVector::built(), clean_up_memory(), oomph::TimingHelpers::convert_secs_to_formatted_string(), oomph::LinearAlgebraDistribution::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::Doc_time, oomph::LinearSolver::Enable_resolve, factorise(), get_total_needed_memory(), int(), oomph::DoubleMatrixBase::ncol(), oomph::CRDoubleMatrix::nnz(), oomph::DistributableLinearAlgebraObject::nrow(), oomph::CRDoubleMatrix::nrow(), oomph::DoubleMatrixBase::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, Solution_time, and oomph::TimingHelpers::timer().

solve() [2/2]

void oomph::SuperLUSolver::solve ( Problem *const &  problem_pt, DoubleVector &  result  )

virtual

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.

Implements oomph::LinearSolver.

  {
    // wipe memory
    this->clean_up_memory();
 
#ifdef OOMPH_HAS_MPI
    // USING SUPERLU DIST
    if (Solver_type == Distributed ||
        (Solver_type == Default && problem_pt->communicator_pt()->nproc() > 1))
    {
      // init the timers
      double t_start = TimingHelpers::timer();
 
      // number of dofs
      unsigned n_dof = problem_pt->ndof();
 
      // set the distribution
      LinearAlgebraDistribution dist(
        problem_pt->communicator_pt(), n_dof, !Dist_use_global_solver);
      this->build_distribution(dist);
 
      // Take a copy of Delete_matrix_data
      bool copy_of_Delete_matrix_data = Dist_delete_matrix_data;
 
      // Set Delete_matrix to true
      Dist_delete_matrix_data = true;
 
      // Use the distributed version of SuperLU_DIST?
      if (!Dist_use_global_solver)
      {
        // Initialise timer
        double t_start = TimingHelpers::timer();
 
        // Storage for the residuals vector
        DoubleVector residuals(this->distribution_pt(), 0.0);
 
        // Get the sparse jacobian and residuals of the problem
        CRDoubleMatrix jacobian(this->distribution_pt());
        problem_pt->get_jacobian(residuals, jacobian);
 
        // Doc time for setup
        double t_end = TimingHelpers::timer();
        Jacobian_setup_time = t_end - t_start;
        if (Doc_time)
        {
          oomph_info << "Time to set up CRDoubleMatrix Jacobian         : "
                     << TimingHelpers::convert_secs_to_formatted_string(
                          Jacobian_setup_time)
                     << std::endl;
        }
 
        // Now call the linear algebra solve, if desired
        if (!Suppress_solve)
        {
          // If the distribution of the result has been build and
          // does not match that of
          // the solver then redistribute before the solve and return
          // to the incoming distribution afterwards.
          if ((result.built()) &&
              (!(*result.distribution_pt() == *this->distribution_pt())))
          {
            LinearAlgebraDistribution temp_global_dist(
              result.distribution_pt());
            result.build(this->distribution_pt(), 0.0);
            solve(&jacobian, residuals, result);
            result.redistribute(&temp_global_dist);
          }
          else
          {
            solve(&jacobian, residuals, result);
          }
        }
      }
      // Otherwise its the global solve version
      else
      {
        // Storage for the residuals vector
        // A non-distriubted residuals vector
        LinearAlgebraDistribution dist(
          problem_pt->communicator_pt(), problem_pt->ndof(), false);
        DoubleVector residuals(&dist, 0.0);
        CRDoubleMatrix jacobian(&dist);
 
        // Get the sparse jacobian and residuals of the problem
        problem_pt->get_jacobian(residuals, jacobian);
 
        // Doc time for setup
        double t_end = TimingHelpers::timer();
        Jacobian_setup_time = t_end - t_start;
        if (Doc_time)
        {
          oomph_info << "Time to set up CR Jacobian    : "
                     << TimingHelpers::convert_secs_to_formatted_string(
                          Jacobian_setup_time)
                     << std::endl;
        }
 
        // Now call the linear algebra solve, if desired
        if (!Suppress_solve)
        {
          // If the result distribution has been built and
          // does not match the global distribution
          // the redistribute before the solve and then return to the
          // distributed version afterwards
          if ((result.built()) && (!(*result.distribution_pt() == dist)))
          {
            LinearAlgebraDistribution temp_global_dist(
              result.distribution_pt());
            result.build(&dist, 0.0);
            solve(&jacobian, residuals, result);
            result.redistribute(&temp_global_dist);
          }
          else
          {
            solve(&jacobian, residuals, result);
          }
        }
      }
      // Set Delete_matrix back to original value
      Dist_delete_matrix_data = copy_of_Delete_matrix_data;
    }
 
    // OTHERWISE WE ARE USING SUPERLU (SERIAL)
    else
#endif
    {
      // set the solver distribution
      LinearAlgebraDistribution dist(
        problem_pt->communicator_pt(), problem_pt->ndof(), false);
      this->build_distribution(dist);
 
      // Allocate storage for the residuals vector
      DoubleVector residuals(dist, 0.0);
 
      // Use the compressed row version?
      if (Serial_compressed_row_flag)
      {
        // Initialise timer
        double t_start = TimingHelpers::timer();
 
        // Get the sparse jacobian and residuals of the problem
        CRDoubleMatrix CR_jacobian(this->distribution_pt());
        problem_pt->get_jacobian(residuals, CR_jacobian);
 
        // If we want to compute the gradient for the globally convergent
        // Newton method, then do it here
        if (Compute_gradient)
        {
          // Compute it
          CR_jacobian.multiply_transpose(residuals,
                                         Gradient_for_glob_conv_newton_solve);
          // Set the flag
          Gradient_has_been_computed = true;
        }
 
        // Doc time for setup
        double t_end = TimingHelpers::timer();
        Jacobian_setup_time = t_end - t_start;
        if (Doc_time)
        {
          oomph_info << std::endl
                     << "Time to set up CRDoubleMatrix Jacobian : "
                     << TimingHelpers::convert_secs_to_formatted_string(
                          Jacobian_setup_time)
                     << std::endl;
        }
 
        // Now call the linear algebra solve, if desired
        if (!Suppress_solve)
        {
          // If the result vector is built and distributed
          // then need to redistribute into the same form as the
          // RHS (non-distributed)
          if ((result.built()) &&
              (!(*result.distribution_pt() == *this->distribution_pt())))
          {
            LinearAlgebraDistribution temp_global_dist(
              result.distribution_pt());
            result.build(this->distribution_pt(), 0.0);
            solve(&CR_jacobian, residuals, result);
            result.redistribute(&temp_global_dist);
          }
          // Otherwise just solve
          else
          {
            solve(&CR_jacobian, residuals, result);
          }
        }
      }
      // Otherwise its the compressed column version
      else
      {
        // Initialise timer
        double t_start = TimingHelpers::timer();
 
        // Get the sparse jacobian and residuals of the problem
        CCDoubleMatrix CC_jacobian;
        problem_pt->get_jacobian(residuals, CC_jacobian);
 
        // If we want to compute the gradient for the globally convergent
        // Newton method, then do it here
        if (Compute_gradient)
        {
          // Compute it
          CC_jacobian.multiply_transpose(residuals,
                                         Gradient_for_glob_conv_newton_solve);
          // Set the flag
          Gradient_has_been_computed = true;
        }
 
        // Doc time for setup
        double t_end = TimingHelpers::timer();
        Jacobian_setup_time = t_end - t_start;
        if (Doc_time)
        {
          oomph_info << "\nTime to set up CCDoubleMatrix Jacobian: "
                     << TimingHelpers::convert_secs_to_formatted_string(
                          Jacobian_setup_time)
                     << std::endl;
        }
 
        // Now call the linear algebra solve, if desired
        if (!Suppress_solve)
        {
          // If the result vector is built and distributed
          // then need to redistribute into the same form as the
          // RHS
          if ((result.built()) &&
              (!(*result.distribution_pt() == *this->distribution_pt())))
          {
            LinearAlgebraDistribution temp_global_dist(
              result.distribution_pt());
            result.build(this->distribution_pt(), 0.0);
            solve(&CC_jacobian, residuals, result);
            result.redistribute(&temp_global_dist);
          }
          // Otherwise just solve
          else
          {
            solve(&CC_jacobian, residuals, result);
          }
        }
      }
 
      // Set the sign of the jacobian
      //(this is computed in the LU decomposition phase)
      problem_pt->sign_of_jacobian() = Serial_sign_of_determinant_of_matrix;
    }
  }

References oomph::DoubleVector::build(), oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DoubleVector::built(), clean_up_memory(), oomph::Problem::communicator_pt(), oomph::LinearSolver::Compute_gradient, oomph::TimingHelpers::convert_secs_to_formatted_string(), Default, Distributed, oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::Doc_time, oomph::Problem::get_jacobian(), oomph::LinearSolver::Gradient_for_glob_conv_newton_solve, oomph::LinearSolver::Gradient_has_been_computed, Jacobian_setup_time, oomph::CRDoubleMatrix::multiply_transpose(), oomph::CCDoubleMatrix::multiply_transpose(), oomph::Problem::ndof(), oomph::OomphCommunicator::nproc(), oomph::oomph_info, oomph::DoubleVector::redistribute(), Serial_compressed_row_flag, Serial_sign_of_determinant_of_matrix, oomph::Problem::sign_of_jacobian(), Solver_type, Suppress_solve, and oomph::TimingHelpers::timer().

solve_transpose() [1/2]

void oomph::SuperLUSolver::solve_transpose ( DoubleMatrixBase *const &  matrix_pt, const DoubleVector &  rhs, DoubleVector &  result  )

virtual

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.

Linear-algebra-type solver: Takes pointer to a matrix and rhs vector and returns the solution of the linear system. Problem pointer defaults to NULL and can be omitted. 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.

Reimplemented from oomph::LinearSolver.

  {
    // Initialise timer
    double t_start = TimingHelpers::timer();
 
    // Pointer used in various places
    CRDoubleMatrix* cr_pt = 0;
 
#ifdef PARANOID
    // check that the rhs vector is setup
    if (!rhs.built())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The vectors rhs must be setup";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // check that the matrix is square
    if (matrix_pt->nrow() != matrix_pt->ncol())
    {
      std::ostringstream error_message_stream;
      error_message_stream << "The matrix at matrix_pt must be square.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // check that the matrix has some entries, and so has a values_pt that
    // makes sense (only for CR because CC is never used I think dense
    // matrices will be safe since they don't use a values pointer).
    cr_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
    if (cr_pt != 0)
    {
      if (cr_pt->nnz() == 0)
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "Attempted to call SuperLu on a CRDoubleMatrix with no entries, "
          << "SuperLU would segfault (because the values array pt is "
          << "uninitialised or null).";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    // check that the matrix and the rhs vector have the same nrow()
    if (matrix_pt->nrow() != rhs.nrow())
    {
      std::ostringstream error_message_stream;
      error_message_stream
        << "The matrix and the rhs vector must have the same number of rows.";
      throw OomphLibError(error_message_stream.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // if the matrix is distributable then should have the same distribution
    // as the rhs vector
    DistributableLinearAlgebraObject* dist_matrix_pt =
      dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
    if (dist_matrix_pt != 0)
    {
      if (!(*dist_matrix_pt->distribution_pt() == *rhs.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The matrix matrix_pt must have the same distribution as the "
          << "rhs vector.";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
    // if the matrix is not distributable then it the rhs vector should not be
    // distributed
    else
    {
      if (rhs.distribution_pt()->distributed())
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The matrix (matrix_pt) is not distributable and therefore the rhs"
          << " vector must not be distributed";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
    // if the result vector is setup then check it has the same distribution
    // as the rhs
    if (result.built())
    {
      if (!(*result.distribution_pt() == *rhs.distribution_pt()))
      {
        std::ostringstream error_message_stream;
        error_message_stream
          << "The result vector distribution has been setup; it must have the "
          << "same distribution as the rhs vector.";
        throw OomphLibError(error_message_stream.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // set the distribution
    if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt))
    {
      // the solver has the same distribution as the matrix if possible
      this->build_distribution(
        dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
          ->distribution_pt());
    }
    else
    {
      // the solver has the same distribution as the RHS
      this->build_distribution(rhs.distribution_pt());
    }
 
    // Doc time for solve
    double t_factorise_start = TimingHelpers::timer();
 
    // Factorise the matrix
    factorise(matrix_pt);
 
    // Doc the end time
    double t_factorise_end = TimingHelpers::timer();
 
    // How long did the factorisation take?
    double factorise_time = t_factorise_end - t_factorise_start;
 
    // Try and upcast the matrix to a CRDoubleMatrix
    // CRDoubleMatrix*
    cr_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
 
    // If the input matrix is a CRDoubleMatrix
    if (cr_pt != 0)
    {
      // ...and actually has an entry
      if (cr_pt->nnz() != 0)
      {
        // Find out how many rows there are in the global Jacobian
        unsigned n_row = cr_pt->nrow();
 
        // And how many non-zeros there are in the global Jacobian
        unsigned n_nnz = cr_pt->nnz();
 
        // Get the memory usage (in bytes) for the global Jacobian storage
        double memory_usage_for_jacobian =
          ((2 * ((n_row + 1) * sizeof(int))) +
           (n_nnz * (sizeof(int) + sizeof(double))));
 
        // Get the memory usage (in bytes) for storage of the LU factors in
        // SuperLU
        double memory_usage_for_lu_storage = get_total_needed_memory();
 
        // Get the memory usage (in bytes) for storage of the LU factors in
        // SuperLU
        double total_memory_usage =
          memory_usage_for_jacobian + memory_usage_for_lu_storage;
 
        // How much memory have we used in the subsidiary preconditioners?
        oomph_info << "\nMemory statistics:"
                   << "\n - Memory used to store the Jacobian (MB): "
                   << memory_usage_for_jacobian / 1.0e+06
                   << "\n - Memory used to store the LU factors (MB): "
                   << memory_usage_for_lu_storage / 1.0e+06
                   << "\n - Total memory used for matrix storage (MB): "
                   << total_memory_usage / 1.0e+06 << "\n"
                   << std::endl;
      }
    } // if (cr_pt!=0)
 
    // Doc the start time
    double t_backsub_start = TimingHelpers::timer();
 
    // Now do the back solve
    backsub_transpose(rhs, result);
 
    // Doc the end time
    double t_backsub_end = TimingHelpers::timer();
 
    // How long did the back substitution take?
    double backsub_time = t_backsub_end - t_backsub_start;
 
    // Doc time for solve
    double t_end = TimingHelpers::timer();
    Solution_time = t_end - t_start;
    if (Doc_time)
    {
      oomph_info
        << "Time for LU factorisation : "
        << TimingHelpers::convert_secs_to_formatted_string(factorise_time)
        << "\nTime for back-substitution: "
        << TimingHelpers::convert_secs_to_formatted_string(backsub_time)
        << "\nTime for SuperLUSolver solve (ndof=" << matrix_pt->nrow() << "): "
        << TimingHelpers::convert_secs_to_formatted_string(Solution_time)
        << std::endl;
    }
 
    // If we are not storing the solver data for resolves, delete it
    if (!Enable_resolve)
    {
      clean_up_memory();
    }
  } // End of solve_transpose

References backsub_transpose(), oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DoubleVector::built(), clean_up_memory(), oomph::TimingHelpers::convert_secs_to_formatted_string(), oomph::LinearAlgebraDistribution::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::Doc_time, oomph::LinearSolver::Enable_resolve, factorise(), get_total_needed_memory(), int(), oomph::DoubleMatrixBase::ncol(), oomph::CRDoubleMatrix::nnz(), oomph::DistributableLinearAlgebraObject::nrow(), oomph::CRDoubleMatrix::nrow(), oomph::DoubleMatrixBase::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, Solution_time, and oomph::TimingHelpers::timer().

solve_transpose() [2/2]

void oomph::SuperLUSolver::solve_transpose ( Problem *const &  problem_pt, DoubleVector &  result  )

virtual

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.

Reimplemented from oomph::LinearSolver.

  {
    // wipe memory
    this->clean_up_memory();
 
#ifdef OOMPH_HAS_MPI
    // USING SUPERLU DIST
    if (Solver_type == Distributed ||
        (Solver_type == Default && problem_pt->communicator_pt()->nproc() > 1))
    {
      // init the timers
      double t_start = TimingHelpers::timer();
 
      // number of dofs
      unsigned n_dof = problem_pt->ndof();
 
      // set the distribution
      LinearAlgebraDistribution dist(
        problem_pt->communicator_pt(), n_dof, !Dist_use_global_solver);
      this->build_distribution(dist);
 
      // Take a copy of Delete_matrix_data
      bool copy_of_Delete_matrix_data = Dist_delete_matrix_data;
 
      // Set Delete_matrix to true
      Dist_delete_matrix_data = true;
 
      // Use the distributed version of SuperLU_DIST?
      if (!Dist_use_global_solver)
      {
        // Initialise timer
        double t_start = TimingHelpers::timer();
 
        // Storage for the residuals vector
        DoubleVector residuals(this->distribution_pt(), 0.0);
 
        // Get the sparse jacobian and residuals of the problem
        CRDoubleMatrix jacobian(this->distribution_pt());
        problem_pt->get_jacobian(residuals, jacobian);
 
        // Doc time for setup
        double t_end = TimingHelpers::timer();
        Jacobian_setup_time = t_end - t_start;
        if (Doc_time)
        {
          oomph_info << "Time to set up CRDoubleMatrix Jacobian         : "
                     << TimingHelpers::convert_secs_to_formatted_string(
                          Jacobian_setup_time)
                     << std::endl;
        }
 
        // Now call the linear algebra solve, if desired
        if (!Suppress_solve)
        {
          // If the distribution of the result has been build and
          // does not match that of
          // the solver then redistribute before the solve and return
          // to the incoming distribution afterwards.
          if ((result.built()) &&
              (!(*result.distribution_pt() == *this->distribution_pt())))
          {
            LinearAlgebraDistribution temp_global_dist(
              result.distribution_pt());
            result.build(this->distribution_pt(), 0.0);
            solve_transpose(&jacobian, residuals, result);
            result.redistribute(&temp_global_dist);
          }
          else
          {
            solve_transpose(&jacobian, residuals, result);
          }
        }
      }
      // Otherwise its the global solve version
      else
      {
        // Storage for the residuals vector
        // A non-distriubted residuals vector
        LinearAlgebraDistribution dist(
          problem_pt->communicator_pt(), problem_pt->ndof(), false);
        DoubleVector residuals(&dist, 0.0);
        CRDoubleMatrix jacobian(&dist);
 
        // Get the sparse jacobian and residuals of the problem
        problem_pt->get_jacobian(residuals, jacobian);
 
        // Doc time for setup
        double t_end = TimingHelpers::timer();
        Jacobian_setup_time = t_end - t_start;
        if (Doc_time)
        {
          oomph_info << "Time to set up CR Jacobian    : "
                     << TimingHelpers::convert_secs_to_formatted_string(
                          Jacobian_setup_time)
                     << std::endl;
        }
 
        // Now call the linear algebra solve, if desired
        if (!Suppress_solve)
        {
          // If the result distribution has been built and
          // does not match the global distribution
          // the redistribute before the solve and then return to the
          // distributed version afterwards
          if ((result.built()) && (!(*result.distribution_pt() == dist)))
          {
            LinearAlgebraDistribution temp_global_dist(
              result.distribution_pt());
            result.build(&dist, 0.0);
            solve_transpose(&jacobian, residuals, result);
            result.redistribute(&temp_global_dist);
          }
          else
          {
            solve_transpose(&jacobian, residuals, result);
          }
        }
      }
      // Set Delete_matrix back to original value
      Dist_delete_matrix_data = copy_of_Delete_matrix_data;
    }
 
    // OTHERWISE WE ARE USING SUPERLU (SERIAL)
    else
#endif
    {
      // set the solver distribution
      LinearAlgebraDistribution dist(
        problem_pt->communicator_pt(), problem_pt->ndof(), false);
      this->build_distribution(dist);
 
      // Allocate storage for the residuals vector
      DoubleVector residuals(dist, 0.0);
 
      // Use the compressed row version?
      if (Serial_compressed_row_flag)
      {
        // Initialise timer
        double t_start = TimingHelpers::timer();
 
        // Get the sparse jacobian and residuals of the problem
        CRDoubleMatrix CR_jacobian(this->distribution_pt());
        problem_pt->get_jacobian(residuals, CR_jacobian);
 
        // If we want to compute the gradient for the globally convergent
        // Newton method, then do it here
        if (Compute_gradient)
        {
          // Compute it
          CR_jacobian.multiply_transpose(residuals,
                                         Gradient_for_glob_conv_newton_solve);
          // Set the flag
          Gradient_has_been_computed = true;
        }
 
        // Doc time for setup
        double t_end = TimingHelpers::timer();
        Jacobian_setup_time = t_end - t_start;
        if (Doc_time)
        {
          oomph_info << std::endl
                     << "Time to set up CRDoubleMatrix Jacobian: "
                     << TimingHelpers::convert_secs_to_formatted_string(
                          Jacobian_setup_time)
                     << std::endl;
        }
 
        // Now call the linear algebra solve, if desired
        if (!Suppress_solve)
        {
          // If the result vector is built and distributed
          // then need to redistribute into the same form as the
          // RHS (non-distributed)
          if ((result.built()) &&
              (!(*result.distribution_pt() == *this->distribution_pt())))
          {
            LinearAlgebraDistribution temp_global_dist(
              result.distribution_pt());
            result.build(this->distribution_pt(), 0.0);
            solve_transpose(&CR_jacobian, residuals, result);
            result.redistribute(&temp_global_dist);
          }
          // Otherwise just solve
          else
          {
            solve_transpose(&CR_jacobian, residuals, result);
          }
        }
      }
      // Otherwise its the compressed column version
      else
      {
        // Initialise timer
        double t_start = TimingHelpers::timer();
 
        // Get the sparse jacobian and residuals of the problem
        CCDoubleMatrix CC_jacobian;
        problem_pt->get_jacobian(residuals, CC_jacobian);
 
        // If we want to compute the gradient for the globally convergent
        // Newton method, then do it here
        if (Compute_gradient)
        {
          // Compute it
          CC_jacobian.multiply_transpose(residuals,
                                         Gradient_for_glob_conv_newton_solve);
          // Set the flag
          Gradient_has_been_computed = true;
        }
 
        // Doc time for setup
        double t_end = TimingHelpers::timer();
        Jacobian_setup_time = t_end - t_start;
        if (Doc_time)
        {
          oomph_info << "\nTime to set up CCDoubleMatrix Jacobian: "
                     << TimingHelpers::convert_secs_to_formatted_string(
                          Jacobian_setup_time)
                     << std::endl;
        }
 
        // Now call the linear algebra solve, if desired
        if (!Suppress_solve)
        {
          // If the result vector is built and distributed
          // then need to redistribute into the same form as the
          // RHS
          if ((result.built()) &&
              (!(*result.distribution_pt() == *this->distribution_pt())))
          {
            LinearAlgebraDistribution temp_global_dist(
              result.distribution_pt());
            result.build(this->distribution_pt(), 0.0);
            solve_transpose(&CC_jacobian, residuals, result);
            result.redistribute(&temp_global_dist);
          }
          // Otherwise just solve
          else
          {
            solve_transpose(&CC_jacobian, residuals, result);
          }
        }
      }
 
      // Set the sign of the jacobian
      //(this is computed in the LU decomposition phase)
      problem_pt->sign_of_jacobian() = Serial_sign_of_determinant_of_matrix;
    }
  }

References oomph::DoubleVector::build(), oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DoubleVector::built(), clean_up_memory(), oomph::Problem::communicator_pt(), oomph::LinearSolver::Compute_gradient, oomph::TimingHelpers::convert_secs_to_formatted_string(), Default, Distributed, oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::Doc_time, oomph::Problem::get_jacobian(), oomph::LinearSolver::Gradient_for_glob_conv_newton_solve, oomph::LinearSolver::Gradient_has_been_computed, Jacobian_setup_time, oomph::CRDoubleMatrix::multiply_transpose(), oomph::CCDoubleMatrix::multiply_transpose(), oomph::Problem::ndof(), oomph::OomphCommunicator::nproc(), oomph::oomph_info, oomph::DoubleVector::redistribute(), Serial_compressed_row_flag, Serial_sign_of_determinant_of_matrix, oomph::Problem::sign_of_jacobian(), Solver_type, Suppress_solve, and oomph::TimingHelpers::timer().

use_compressed_column_for_superlu_serial()

void oomph::SuperLUSolver::use_compressed_column_for_superlu_serial ( )

inline

Use the compressed column format in superlu serial.

    {
      Serial_compressed_row_flag = false;
    }

References Serial_compressed_row_flag.

use_compressed_row_for_superlu_serial()

void oomph::SuperLUSolver::use_compressed_row_for_superlu_serial ( )

inline

Use the compressed row format in superlu serial.

    {
      Serial_compressed_row_flag = true;
    }

References Serial_compressed_row_flag.

Member Data Documentation

Doc_stats

bool oomph::SuperLUSolver::Doc_stats

private

Set to true to output statistics (false by default).

Referenced by backsub_serial(), backsub_transpose_serial(), clean_up_memory(), disable_doc_stats(), enable_doc_stats(), factorise_serial(), and SuperLUSolver().

Jacobian_setup_time

double oomph::SuperLUSolver::Jacobian_setup_time

private

Jacobian setup time.

Referenced by jacobian_setup_time(), solve(), and solve_transpose().

Serial_compressed_row_flag

bool oomph::SuperLUSolver::Serial_compressed_row_flag

private

Use compressed row version?

Referenced by backsub_serial(), backsub_transpose_serial(), clean_up_memory(), factorise_serial(), solve(), solve_transpose(), SuperLUSolver(), use_compressed_column_for_superlu_serial(), and use_compressed_row_for_superlu_serial().

Serial_f_factors

void* oomph::SuperLUSolver::Serial_f_factors

private

Storage for the LU factors as required by SuperLU.

Referenced by backsub_serial(), backsub_transpose_serial(), clean_up_memory(), factorise_serial(), get_memory_usage_for_lu_factors(), and get_total_needed_memory().

Serial_info

int oomph::SuperLUSolver::Serial_info

private

Info flag for the SuperLU solver.

Referenced by backsub_serial(), backsub_transpose_serial(), clean_up_memory(), and factorise_serial().

Serial_n_dof

unsigned long oomph::SuperLUSolver::Serial_n_dof

private

The number of unknowns in the linear system.

Referenced by backsub_serial(), backsub_transpose_serial(), clean_up_memory(), factorise_serial(), and SuperLUSolver().

Serial_sign_of_determinant_of_matrix

int oomph::SuperLUSolver::Serial_sign_of_determinant_of_matrix

private

Sign of the determinant of the matrix.

Referenced by factorise_serial(), solve(), solve_transpose(), and SuperLUSolver().

Solution_time

double oomph::SuperLUSolver::Solution_time

private

Solution time.

Referenced by linear_solver_solution_time(), resolve(), resolve_transpose(), solve(), and solve_transpose().

Solver_type

Type oomph::SuperLUSolver::Solver_type

private

the solver type. see SuperLU_solver_type for details.

Referenced by factorise(), get_memory_usage_for_lu_factors(), get_total_needed_memory(), set_solver_type(), solve(), solve_transpose(), and SuperLUSolver().

Suppress_solve

bool oomph::SuperLUSolver::Suppress_solve

private

Suppress solve?

Referenced by solve(), solve_transpose(), and SuperLUSolver().

Using_dist

bool oomph::SuperLUSolver::Using_dist

private

boolean flag indicating whether superlu dist is being used

Referenced by backsub(), backsub_transpose(), factorise(), and SuperLUSolver().

---

The documentation for this class was generated from the following files:

• linear_solver.h
• linear_solver.cc