Wrapper to Mumps solver. More...

#include <mumps_solver.h>

Inheritance diagram for oomph::MumpsSolver:

Public Member Functions
	MumpsSolver ()
	Constructor: Call setup. More...

	MumpsSolver (const MumpsSolver &dummy)=delete
	Broken copy constructor. More...

void	operator= (const MumpsSolver &)=delete
	Broken assignment operator. More...

	~MumpsSolver ()
	Destructor: Cleanup. More...

void	disable_resolve ()
	Overload disable resolve so that it cleans up memory too. More...

void	enable_suppress_warning_about_MPI_COMM_WORLD ()

void	disable_suppress_warning_about_MPI_COMM_WORLD ()
	Don't suppress warning about communicator not equal to MPI_COMM_WORLD. More...

void	enable_suppress_mumps_info_during_solve ()

void	disable_suppress_mumps_info_during_solve ()
	Don't suppress info being printed to screen during MUMPS solve. More...

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

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

void	resolve (const DoubleVector &rhs, DoubleVector &result)

void	enable_doc_stats ()
	Enable documentation of statistics. More...

void	disable_doc_stats ()
	Disable documentation of statistics. More...

double	jacobian_setup_time ()

virtual double	linear_solver_solution_time ()

void	enable_suppress_solve ()

void	disable_suppress_solve ()

void	enable_delete_matrix_data ()

void	disable_delete_matrix_data ()

void	factorise (DoubleMatrixBase *const &matrix_pt)

void	backsub (const DoubleVector &rhs, DoubleVector &result)

void	clean_up_memory ()
	Clean up the memory allocated by the mumps solver. More...

void	declare_jacobian_is_unsymmetric ()
	Tell MUMPS that the Jacobian matrix is unsymmetric. More...

void	declare_jacobian_is_symmetric ()
	Tell MUMPS that the Jacobian matrix is general symmetric. More...

void	declare_jacobian_is_symmetric_positive_definite ()
	Tell MUMPS that the Jacobian matrix is symmetric positive-definite. More...

void	use_pord_ordering ()

void	use_metis_ordering ()

void	use_scotch_ordering ()

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 (Problem *const &problem_pt, DoubleVector &result)

virtual void	solve_transpose (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)

virtual void	solve_transpose (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)

virtual void	resolve_transpose (const DoubleVector &rhs, DoubleVector &result)

virtual double	jacobian_setup_time () const

virtual double	linear_solver_solution_time () const

virtual void	enable_computation_of_gradient ()

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)

Static Public Attributes
static bool	Suppress_incorrect_rhs_distribution_warning_in_resolve

static int	Default_workspace_scaling_factor = 5

Private Types
enum	MumpsJacobianSymmetryFlags { Unsymmetric = 0 , Symmetric_positive_definite = 1 , Symmetric = 2 }

enum	MumpsJacobianOrderingFlags { Scotch_ordering = 3 , Pord_ordering = 4 , Metis_ordering = 5 }

Private Member Functions
void	initialise_mumps ()
	Initialise instance of mumps data structure. More...

void	shutdown_mumps ()
	Shutdown mumps. 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 for MumpsSolver to output statistics (false by default). More...

bool	Suppress_warning_about_MPI_COMM_WORLD

bool	Suppress_mumps_info_during_solve
	Boolean to suppress info being printed to screen during MUMPS solve. More...

bool	Mumps_is_initialised
	Has mumps been initialised? More...

unsigned	Workspace_scaling_factor

bool	Delete_matrix_data

Vector< int >	Irn_loc
	Vector for row numbers. More...

Vector< int >	Jcn_loc

Vector< double >	A_loc

DMUMPS_STRUC_C *	Mumps_struc_pt
	Pointer to MUMPS struct that contains the solver data. More...

unsigned	Jacobian_symmetry_flag

unsigned	Jacobian_ordering_flag

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

Wrapper to Mumps solver.

Member Enumeration Documentation

◆ MumpsJacobianOrderingFlags

enum oomph::MumpsSolver::MumpsJacobianOrderingFlags

private

ordering library to use for serial analysis; magic numbers as defined by MUMPS documentation

Enumerator
Scotch_ordering
Pord_ordering
Metis_ordering

     {
       Scotch_ordering = 3,
       Pord_ordering = 4,
       Metis_ordering = 5
     };

◆ MumpsJacobianSymmetryFlags

enum oomph::MumpsSolver::MumpsJacobianSymmetryFlags

private

values of the SYM variable used by the MUMPS solver which dictates the symmetry properties of the Jacobian matrix; magic numbers as defined by MUMPS documentation

Enumerator
Unsymmetric
Symmetric_positive_definite
Symmetric

     {
       Unsymmetric = 0,
       Symmetric_positive_definite = 1,
       Symmetric = 2
     };

Constructor & Destructor Documentation

◆ MumpsSolver() [1/2]

oomph::MumpsSolver::MumpsSolver ( )

Constructor: Call setup.

     : Suppress_solve(false),
       Doc_stats(false),
       Suppress_warning_about_MPI_COMM_WORLD(false),
       Suppress_mumps_info_during_solve(false),
       Mumps_is_initialised(false),
       Workspace_scaling_factor(Default_workspace_scaling_factor),
       Delete_matrix_data(false),
       Mumps_struc_pt(nullptr),
       Jacobian_symmetry_flag(MumpsJacobianSymmetryFlags::Unsymmetric),
       Jacobian_ordering_flag(MumpsJacobianOrderingFlags::Metis_ordering)
   {
   }

◆ MumpsSolver() [2/2]

oomph::MumpsSolver::MumpsSolver ( const MumpsSolver & dummy )

delete

Broken copy constructor.

◆ ~MumpsSolver()

oomph::MumpsSolver::~MumpsSolver ( )

Destructor: Cleanup.

Destructor: Shutdown mumps.

   {
     shutdown_mumps();
   }

References shutdown_mumps().

Member Function Documentation

◆ backsub()

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

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

Do the backsubstitution for mumps solver. This does not make any assumption about the distribution of the vectors

   {
     double t_start = TimingHelpers::timer();
  
 #ifdef PARANOID
     if (!Suppress_warning_about_MPI_COMM_WORLD)
     {
       if (this->distribution_pt()->communicator_pt()->mpi_comm() !=
           MPI_COMM_WORLD)
       {
         std::ostringstream error_message_stream;
         error_message_stream
           << "Warning: Mumps wrapper assumes that communicator is "
              "MPI_COMM_WORLD\n"
           << "         which is not the case, so mumps may die...\n"
           << "         If it does initialise oomph-lib's mpi without "
              "requesting\n"
           << "         the communicator to be a duplicate of MPI_COMM_WORLD\n"
           << "         (done via an optional boolean to "
              "MPI_Helpers::init(...)\n\n"
           << "         (You can suppress this warning by recompiling without\n"
           << "         paranoia or calling \n\n"
           << "         "
              "MumpsSolver::enable_suppress_warning_about_MPI_COMM_WORLD()\n\n"
           << "         \n";
         OomphLibWarning(error_message_stream.str(),
                         "MumpsSolver::backsub()",
                         OOMPH_EXCEPTION_LOCATION);
       }
     }
  
     // Initialised?
     if (!Mumps_is_initialised)
     {
       std::ostringstream error_message_stream;
       error_message_stream << "Mumps has not been initialised.";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
  
     // check that the rhs vector is setup
     if (!rhs.distribution_pt()->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);
     }
  
 #endif
  
     // Check that the rhs distribution is the same as the distribution as this
     // solver. If not redistribute and issue a warning
     LinearAlgebraDistribution rhs_distribution(rhs.distribution_pt());
     if (!(*rhs.distribution_pt() == *this->distribution_pt()))
     {
       if (!Suppress_incorrect_rhs_distribution_warning_in_resolve)
       {
         std::ostringstream warning_stream;
         warning_stream << "The distribution of rhs vector does not match that "
                           "of the solver.\n";
         warning_stream
           << "The rhs may have to be redistributed but we're not doing this "
              "because\n"
           << "I'm no longer convinced it's necessary. Keep an eye on this...\n";
         warning_stream
           << "To remove this warning you can either:\n"
           << "    i) Ensure that the rhs vector has the correct distribution\n"
           << "       before calling the resolve() function\n"
           << "or ii) Set the flag \n"
           << " MumpsSolver::Suppress_incorrect_rhs_distribution_warning_in_"
              "resolve\n"
           << "       to be true\n\n";
  
         OomphLibWarning(warning_stream.str(),
                         "MumpsSolver::resolve()",
                         OOMPH_EXCEPTION_LOCATION);
       }
  
       // //Have to cast away const-ness (which tells us that we shouldn't really
       // //be doing this!)
       // const_cast<DoubleVector&>(rhs).redistribute(this->distribution_pt());
     }
  
  
 #ifdef PARANOID
     // if the result vector is setup then check it has the same distribution
     // as the rhs
     if (result.distribution_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
  
  
     // Doc stats?
     if (Doc_stats)
     {
       // Output stream for global info on host. Negative value suppresses
       // printing
       Mumps_struc_pt->ICNTL(3) = 6;
     }
  
     // number of DOFs
     int ndof = this->distribution_pt()->nrow();
  
     // Make backup to avoid over-writing
     DoubleVector tmp_rhs;
     tmp_rhs = rhs;
  
     // Now turn this into a global (non-distributed) vector
     // because that's what mumps needs
  
     // Make a global distribution (i.e. one that isn't distributed)
     LinearAlgebraDistribution global_distribution(
       this->distribution_pt()->communicator_pt(), ndof, false);
  
     // Redistribute the tmp_rhs vector with this distribution -- it's
     // now "global", as required for mumps
     tmp_rhs.redistribute(&global_distribution);
  
     // Do the backsubsitution phase -- overwrites the tmp_rhs vector with the
     // solution
     Mumps_struc_pt->rhs = &tmp_rhs[0];
     Mumps_struc_pt->job = 3;
     dmumps_c(Mumps_struc_pt);
  
     // Copy back
     unsigned n = Mumps_struc_pt->n;
     for (unsigned j = 0; j < n; j++)
     {
       tmp_rhs[j] = Mumps_struc_pt->rhs[j];
     }
  
     // Broadcast the result which is only held on root
     MPI_Bcast(&tmp_rhs[0],
               ndof,
               MPI_DOUBLE,
               0,
               this->distribution_pt()->communicator_pt()->mpi_comm());
  
     // If the result vector is distributed, re-distribute the
     // non-distributed tmp_rhs vector to match
     if (result.built())
     {
       tmp_rhs.redistribute(result.distribution_pt());
     }
     else
     {
       tmp_rhs.redistribute(this->distribution_pt());
     }
  
     // Now copy the tmp_rhs vector into the (matching) result
     result = tmp_rhs;
  
     if ((Doc_time) &&
         (this->distribution_pt()->communicator_pt()->my_rank() == 0))
     {
       double t_end = TimingHelpers::timer();
       oomph_info << "Time for MumpsSolver backsub                            : "
                  << TimingHelpers::convert_secs_to_formatted_string(t_end -
                                                                     t_start)
                  << std::endl;
     }
  
     // Switch off docing again by setting output stream for global info on
     // to negative number
     Mumps_struc_pt->ICNTL(3) = -1;
   }

References oomph::DoubleVector::built(), oomph::LinearAlgebraDistribution::built(), oomph::TimingHelpers::convert_secs_to_formatted_string(), oomph::DistributableLinearAlgebraObject::distribution_built(), oomph::DistributableLinearAlgebraObject::distribution_pt(), Doc_stats, oomph::LinearSolver::Doc_time, j, Mumps_is_initialised, Mumps_struc_pt, n, oomph::LinearAlgebraDistribution::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::DoubleVector::redistribute(), Suppress_incorrect_rhs_distribution_warning_in_resolve, Suppress_warning_about_MPI_COMM_WORLD, and oomph::TimingHelpers::timer().

Referenced by resolve(), and solve().

◆ clean_up_memory()

void oomph::MumpsSolver::clean_up_memory ( )

virtual

Clean up the memory allocated by the mumps solver.

Clean up the memory allocated for the MumpsSolver solver.

Reimplemented from oomph::LinearSolver.

   {
     shutdown_mumps();
   }

References shutdown_mumps().

Referenced by oomph::NewMumpsPreconditioner::clean_up_memory(), disable_resolve(), and solve().

◆ declare_jacobian_is_symmetric()

void oomph::MumpsSolver::declare_jacobian_is_symmetric ( )

inline

Tell MUMPS that the Jacobian matrix is general symmetric.

     {
       Jacobian_symmetry_flag = Symmetric;
     }

References Jacobian_symmetry_flag, and Symmetric.

◆ declare_jacobian_is_symmetric_positive_definite()

void oomph::MumpsSolver::declare_jacobian_is_symmetric_positive_definite ( )

inline

Tell MUMPS that the Jacobian matrix is symmetric positive-definite.

     {
       Jacobian_symmetry_flag = Symmetric_positive_definite;
     }

References Jacobian_symmetry_flag, and Symmetric_positive_definite.

◆ declare_jacobian_is_unsymmetric()

void oomph::MumpsSolver::declare_jacobian_is_unsymmetric ( )

inline

Tell MUMPS that the Jacobian matrix is unsymmetric.

     {
       Jacobian_symmetry_flag = Unsymmetric;
     }

References Jacobian_symmetry_flag, and Unsymmetric.

◆ disable_delete_matrix_data()

void oomph::MumpsSolver::disable_delete_matrix_data ( )

inline

Unset Delete_matrix_data flag. MumpsSolver 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.

     {
       Delete_matrix_data = false;
     }

References Delete_matrix_data.

◆ disable_doc_stats()

void oomph::MumpsSolver::disable_doc_stats ( )

inline

Disable documentation of statistics.

     {
       Doc_stats = false;
     }

References Doc_stats.

◆ disable_resolve()

void oomph::MumpsSolver::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().

◆ disable_suppress_mumps_info_during_solve()

void oomph::MumpsSolver::disable_suppress_mumps_info_during_solve ( )

inline

Don't suppress info being printed to screen during MUMPS solve.

     {
       Suppress_mumps_info_during_solve = false;
     }

References Suppress_mumps_info_during_solve.

◆ disable_suppress_solve()

void oomph::MumpsSolver::disable_suppress_solve ( )

inline

Unset the flag so that the system is actually solved again This is the default (obviously)

     {
       Suppress_solve = false;
     }

References Suppress_solve.

◆ disable_suppress_warning_about_MPI_COMM_WORLD()

void oomph::MumpsSolver::disable_suppress_warning_about_MPI_COMM_WORLD ( )

inline

Don't suppress warning about communicator not equal to MPI_COMM_WORLD.

     {
       Suppress_warning_about_MPI_COMM_WORLD = false;
     }

References Suppress_warning_about_MPI_COMM_WORLD.

◆ enable_delete_matrix_data()

void oomph::MumpsSolver::enable_delete_matrix_data ( )

inline

Set Delete_matrix_data flag. MumpsSolver 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.

     {
       Delete_matrix_data = true;
     }

References Delete_matrix_data.

◆ enable_doc_stats()

void oomph::MumpsSolver::enable_doc_stats ( )

inline

Enable documentation of statistics.

     {
       Doc_stats = true;
     }

References Doc_stats.

◆ enable_suppress_mumps_info_during_solve()

void oomph::MumpsSolver::enable_suppress_mumps_info_during_solve ( )

inline

Set boolean to suppress info being printed to screen during MUMPS solve

     {
       Suppress_mumps_info_during_solve = true;
     }

References Suppress_mumps_info_during_solve.

◆ enable_suppress_solve()

void oomph::MumpsSolver::enable_suppress_solve ( )

inline

Set the flag to avoid solution of the system, so just assemble the Jacobian and the rhs. (Used only for timing runs, obviously)

     {
       Suppress_solve = true;
     }

References Suppress_solve.

◆ enable_suppress_warning_about_MPI_COMM_WORLD()

void oomph::MumpsSolver::enable_suppress_warning_about_MPI_COMM_WORLD ( )

inline

Set boolean to suppress warning about communicator not equal to MPI_COMM_WORLD

     {
       Suppress_warning_about_MPI_COMM_WORLD = true;
     }

References Suppress_warning_about_MPI_COMM_WORLD.

◆ factorise()

void oomph::MumpsSolver::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 using mumps. The resulting matrix factors are stored internally. Note: if Delete_matrix_data is true the function matrix_pt->clean_up_memory() will be used to wipe the matrix data.

   {
     // Initialise timer
     double t_start = TimingHelpers::timer();
  
     // set the distribution
     DistributableLinearAlgebraObject* dist_matrix_pt =
       dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
     if (dist_matrix_pt)
     {
       // the solver has the same distribution as the matrix if possible
       this->build_distribution(dist_matrix_pt->distribution_pt());
     }
     else
     {
       throw OomphLibError("Matrix must be distributable",
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
  
  
     // Check that we have a square matrix
 #ifdef PARANOID
     int n = matrix_pt->nrow();
     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);
     }
     if (!Suppress_warning_about_MPI_COMM_WORLD)
     {
       if (this->distribution_pt()->communicator_pt()->mpi_comm() !=
           MPI_COMM_WORLD)
       {
         std::ostringstream error_message_stream;
         error_message_stream
           << "Warning: Mumps wrapper assumes that communicator is "
              "MPI_COMM_WORLD\n"
           << "         which is not the case, so mumps may die...\n"
           << "         If it does initialise oomph-lib's mpi without "
              "requesting\n"
           << "         the communicator to be a duplicate of MPI_COMM_WORLD\n"
           << "         (done via an optional boolean to "
              "MPI_Helpers::init(...)\n\n"
           << "         (You can suppress this warning by recompiling without\n"
           << "         paranoia or calling \n\n"
           << "         "
              "MumpsSolver::enable_suppress_warning_about_MPI_COMM_WORLD()\n\n"
           << "         \n";
         OomphLibWarning(error_message_stream.str(),
                         "MumpsSolver::factorise()",
                         OOMPH_EXCEPTION_LOCATION);
       }
     }
  
 #endif
  
     // Initialise
     if (Mumps_is_initialised)
     {
       shutdown_mumps();
     }
     initialise_mumps();
  
  
     // Doc stats?
     if (Doc_stats)
     {
       // Output stream for global info on host. Negative value suppresses
       // printing
       Mumps_struc_pt->ICNTL(3) = 6;
     }
  
     // number of processors
     unsigned nproc = 1;
     if (dist_matrix_pt != 0)
     {
       nproc = dist_matrix_pt->distribution_pt()->communicator_pt()->nproc();
     }
  
     // Is it a CRDoubleMatrix?
     CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt);
     if (cr_matrix_pt != 0)
     {
 #ifdef PARANOID
       // paranoid check that the matrix has been set up
       if (!cr_matrix_pt->built())
       {
         throw OomphLibError(
           "To apply MumpsSolver to a CRDoubleMatrix - it must be built",
           OOMPH_CURRENT_FUNCTION,
           OOMPH_EXCEPTION_LOCATION);
       }
 #endif
  
       // if the matrix is distributed then set up solver
       if ((nproc == 1) || (cr_matrix_pt->distributed()))
       {
         double t_start_copy = TimingHelpers::timer();
  
         // Find the number of rows and non-zero entries in the matrix
         const int nnz_loc = int(cr_matrix_pt->nnz());
         const int n = matrix_pt->nrow();
  
         // Create mumps storage
  
         // Create vector for row numbers
         Irn_loc.resize(nnz_loc);
  
         // Create vector for column numbers
         Jcn_loc.resize(nnz_loc);
  
         // Vector for entries
         A_loc.resize(nnz_loc);
  
         // First row held on this processor
         int first_row = cr_matrix_pt->first_row();
  
         // Copy into coordinate storage scheme using pointer arithmetic
         double* matrix_value_pt = cr_matrix_pt->value();
         int* matrix_index_pt = cr_matrix_pt->column_index();
         int* matrix_start_pt = cr_matrix_pt->row_start();
         int i_row = 0;
  
         // is the matrix symmetric? If so, we must only provide
         // MUMPS with the upper (or lower) triangular part of the Jacobian
         if (Mumps_struc_pt->sym != MumpsJacobianSymmetryFlags::Unsymmetric)
         {
           oomph_info << "Assuming Jacobian matrix is symmetric "
                      << "(passing MUMPS the upper triangular part)"
                      << std::endl;
         }
  
         for (int count = 0; count < nnz_loc; count++)
         {
           A_loc[count] = matrix_value_pt[count];
           Jcn_loc[count] = matrix_index_pt[count] + 1;
           if (count < matrix_start_pt[i_row + 1])
           {
             Irn_loc[count] = first_row + i_row + 1;
           }
           else
           {
             i_row++;
             Irn_loc[count] = first_row + i_row + 1;
           }
  
           // only pass the upper triangular part (and diagonal)
           // if we have a symmetric matrix (MUMPS sums values,
           // so need to set the lwoer triangle to zero)
           if ((Mumps_struc_pt->sym !=
                MumpsJacobianSymmetryFlags::Unsymmetric) &&
               (Irn_loc[count] > Jcn_loc[count]))
           {
             A_loc[count] = 0.0;
           }
         }
  
         // Now delete the matrix if we are allowed
         if (Delete_matrix_data == true)
         {
           cr_matrix_pt->clear();
         }
  
         if ((Doc_time) &&
             (this->distribution_pt()->communicator_pt()->my_rank() == 0))
         {
           double t_end_copy = TimingHelpers::timer();
           oomph_info << "Time for copying matrix into MumpsSolver data "
                         "structure : "
                      << TimingHelpers::convert_secs_to_formatted_string(
                           t_end_copy - t_start_copy)
                      << std::endl;
         }
  
         // Call mumps factorisation
         //-------------------------
  
         // Specify size of system
         Mumps_struc_pt->n = n;
  
         // Number of nonzeroes
         Mumps_struc_pt->nz_loc = nnz_loc;
  
         // The entries
         Mumps_struc_pt->irn_loc = &Irn_loc[0];
         Mumps_struc_pt->jcn_loc = &Jcn_loc[0];
         Mumps_struc_pt->a_loc = &A_loc[0];
  
         double t_start_analyse = TimingHelpers::timer();
  
         // Do analysis
         Mumps_struc_pt->job = 1;
         dmumps_c(Mumps_struc_pt);
  
  
         if ((Doc_time) &&
             (this->distribution_pt()->communicator_pt()->my_rank() == 0))
         {
           double t_end_analyse = TimingHelpers::timer();
           oomph_info
             << "Time for mumps analysis stage in MumpsSolver            : "
             << TimingHelpers::convert_secs_to_formatted_string(t_end_analyse -
                                                                t_start_analyse)
             << "\n(Ordering generated using: ";
  
           switch (Mumps_struc_pt->INFOG(7))
           {
             case MumpsJacobianOrderingFlags::Scotch_ordering:
               oomph_info << "SCOTCH";
               break;
             case MumpsJacobianOrderingFlags::Pord_ordering:
               oomph_info << "PORD";
               break;
             case MumpsJacobianOrderingFlags::Metis_ordering:
               oomph_info << "METIS";
               break;
             default:
               oomph_info << Mumps_struc_pt->INFOG(7);
           }
  
           oomph_info << ")" << std::endl;
         }
  
  
         int my_rank = this->distribution_pt()->communicator_pt()->my_rank();
  
         // Document estimate for size of workspace
         if (my_rank == 0)
         {
           if (!Suppress_mumps_info_during_solve)
           {
             oomph_info << "Estimated max. workspace in MB: "
                        << Mumps_struc_pt->INFOG(26) << std::endl;
           }
         }
  
         double t_start_factor = TimingHelpers::timer();
  
         // Loop until successfully factorised
         bool factorised = false;
         while (!factorised)
         {
           // Set workspace to multiple of that -- ought to be "significantly
           // larger than infog(26)", according to the manual :(
           Mumps_struc_pt->ICNTL(23) =
             Workspace_scaling_factor * Mumps_struc_pt->INFOG(26);
  
           if (!Suppress_mumps_info_during_solve)
           {
             oomph_info << "Attempting factorisation with workspace estimate: "
                        << Mumps_struc_pt->ICNTL(23) << " MB\n";
             oomph_info << "corresponding to Workspace_scaling_factor = "
                        << Workspace_scaling_factor << "\n";
           }
  
           // Do factorisation
           Mumps_struc_pt->job = 2;
           dmumps_c(Mumps_struc_pt);
  
           // Check for error
           if (Mumps_struc_pt->INFOG(1) != 0)
           {
             if (!Suppress_mumps_info_during_solve)
             {
               oomph_info << "Error during mumps factorisation!\n";
               oomph_info << "Error codes: " << Mumps_struc_pt->INFO(1) << " "
                          << Mumps_struc_pt->INFO(2) << std::endl;
             }
  
             // Increase scaling factor for workspace and run again
             Workspace_scaling_factor *= 2;
  
             if (!Suppress_mumps_info_during_solve)
             {
               oomph_info << "Increasing workspace_scaling_factor to "
                          << Workspace_scaling_factor << std::endl;
             }
           }
           else
           {
             factorised = true;
           }
         }
  
  
         if ((Doc_time) &&
             (this->distribution_pt()->communicator_pt()->my_rank() == 0))
         {
           double t_end_factor = TimingHelpers::timer();
           oomph_info << "Time for actual mumps factorisation in MumpsSolver"
                         "      : "
                      << TimingHelpers::convert_secs_to_formatted_string(
                           t_end_factor - t_start_factor)
                      << std::endl;
         }
       }
       // else the CRDoubleMatrix is not distributed
       else
       {
         std::ostringstream error_message_stream;
         error_message_stream << "MumpsSolver only works for a "
                              << " distributed CRDoubleMatrix\n";
         throw OomphLibError(error_message_stream.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
       }
     }
     // Otherwise throw an error
     else
     {
       std::ostringstream error_message_stream;
       error_message_stream << "MumpsSolver implemented only for "
                            << "distributed CRDoubleMatrix. \n";
       throw OomphLibError(error_message_stream.str(),
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
     }
  
  
     if ((Doc_time) &&
         (this->distribution_pt()->communicator_pt()->my_rank() == 0))
     {
       double t_end = TimingHelpers::timer();
       oomph_info << "Time for MumpsSolver factorisation                      : "
                  << TimingHelpers::convert_secs_to_formatted_string(t_end -
                                                                     t_start)
                  << std::endl;
     }
  
     // Switch off docing again by setting output stream for global info on
     // to negative number
     Mumps_struc_pt->ICNTL(3) = -1;
   }

References A_loc, oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::CRDoubleMatrix::built(), oomph::CRDoubleMatrix::clear(), oomph::CRDoubleMatrix::column_index(), oomph::LinearAlgebraDistribution::communicator_pt(), oomph::TimingHelpers::convert_secs_to_formatted_string(), Delete_matrix_data, oomph::DistributableLinearAlgebraObject::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), Doc_stats, oomph::LinearSolver::Doc_time, oomph::DistributableLinearAlgebraObject::first_row(), initialise_mumps(), int(), Irn_loc, Jcn_loc, m, Mumps_is_initialised, Mumps_struc_pt, oomph::OomphCommunicator::my_rank(), n, oomph::DoubleMatrixBase::ncol(), oomph::CRDoubleMatrix::nnz(), oomph::OomphCommunicator::nproc(), oomph::DoubleMatrixBase::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::CRDoubleMatrix::row_start(), shutdown_mumps(), Suppress_mumps_info_during_solve, Suppress_warning_about_MPI_COMM_WORLD, oomph::TimingHelpers::timer(), oomph::CRDoubleMatrix::value(), and Workspace_scaling_factor.

Referenced by oomph::NewMumpsPreconditioner::setup(), and solve().

◆ initialise_mumps()

void oomph::MumpsSolver::initialise_mumps ( )

private

Initialise instance of mumps data structure.

   {
     // Make new instance of Mumps data structure
     Mumps_struc_pt = new DMUMPS_STRUC_C;
  
     // Mumps' hack to indicate that mpi_comm_world is used. Conversion of any
     // other communicators appears to be non-portable, so we simply
     // issue a warning if we later discover that oomph-lib's communicator
     // is not MPI_COMM_WORLD
     Mumps_struc_pt->comm_fortran = -987654;
  
     // Root processor participates in solution
     Mumps_struc_pt->par = 1;
  
     // Set matrix symmetry properties
     // (unsymmetric / symmetric positive definite / general symmetric)
     Mumps_struc_pt->sym = Jacobian_symmetry_flag;
  
     // First call does the initialise phase
     Mumps_struc_pt->job = -1;
  
     // Call c interface to double precision mumps to initialise
     dmumps_c(Mumps_struc_pt);
     Mumps_is_initialised = true;
  
     // Output stream for global info on host. Negative value suppresses printing
     Mumps_struc_pt->ICNTL(3) = -1;
  
     // Only show error messages and stats
     // (4 for full doc; creates huge amount of output)
     Mumps_struc_pt->ICNTL(4) = 2;
  
     // Write matrix
     // sprintf(Mumps_struc_pt->write_problem,"/work/e173/e173/mheil/matrix");
  
     // Assembled matrix (rather than element-by_element)
     Mumps_struc_pt->ICNTL(5) = 0;
  
     // set the package to use for ordering during (sequential) analysis
     Mumps_struc_pt->ICNTL(7) = Jacobian_ordering_flag;
  
     // Distributed problem with user-specified distribution
     Mumps_struc_pt->ICNTL(18) = 3;
  
     // Dense RHS
     Mumps_struc_pt->ICNTL(20) = 0;
  
     // Non-distributed solution
     Mumps_struc_pt->ICNTL(21) = 0;
   }

References Jacobian_ordering_flag, Jacobian_symmetry_flag, Mumps_is_initialised, and Mumps_struc_pt.

Referenced by factorise().

◆ jacobian_setup_time()

double oomph::MumpsSolver::jacobian_setup_time ( )

inline

Returns the time taken to assemble the Jacobian matrix and residual vector

     {
       return Jacobian_setup_time;
     }

References Jacobian_setup_time.

◆ linear_solver_solution_time()

virtual double oomph::MumpsSolver::linear_solver_solution_time ( )

inlinevirtual

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

     {
       return Solution_time;
     }

References Solution_time.

◆ operator=()

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

delete

Broken assignment operator.

◆ resolve()

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

virtual

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.

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.

   {
 #ifdef PARANOID
     if (!Suppress_warning_about_MPI_COMM_WORLD)
     {
       if (this->distribution_pt()->communicator_pt()->mpi_comm() !=
           MPI_COMM_WORLD)
       {
         std::ostringstream error_message_stream;
         error_message_stream
           << "Warning: Mumps wrapper assumes communicator is MPI_COMM_WORLD\n"
           << "         which is not the case, so mumps may die...\n"
           << "         If it does initialise oomph-lib's mpi without "
              "requesting\n"
           << "         the communicator to be a duplicate of MPI_COMM_WORLD\n"
           << "         (done via an optional boolean to "
              "MPI_Helpers::init(...)\n\n"
           << "         (You can suppress this warning by recompiling without\n"
           << "         paranoia or calling\n\n"
           << "         "
              "MumpsSolver::enable_suppress_warning_about_MPI_COMM_WORLD()\n\n"
           << "         \n";
         OomphLibWarning(error_message_stream.str(),
                         "MumpsSolver::resolve()",
                         OOMPH_EXCEPTION_LOCATION);
       }
     }
 #endif
  
     // Store starting time for solve
     double t_start = TimingHelpers::timer();
  
     // Doc stats?
     if (Doc_stats)
     {
       // Output stream for global info on host. Negative value suppresses
       // printing
       Mumps_struc_pt->ICNTL(3) = 6;
     }
  
     // Now do the back substitution phase
     backsub(rhs, result);
  
     // Doc time for solve
     double t_end = TimingHelpers::timer();
     Solution_time = t_end - t_start;
  
     // Switch off docing again by setting output stream for global info on
     // to negative number
     Mumps_struc_pt->ICNTL(3) = -1;
  
     if ((Doc_time) &&
         (this->distribution_pt()->communicator_pt()->my_rank() == 0))
     {
       oomph_info << "Time for MumpsSolver solve: "
                  << TimingHelpers::convert_secs_to_formatted_string(t_end -
                                                                     t_start)
                  << std::endl;
     }
   }

References backsub(), oomph::TimingHelpers::convert_secs_to_formatted_string(), oomph::DistributableLinearAlgebraObject::distribution_pt(), Doc_stats, oomph::LinearSolver::Doc_time, Mumps_struc_pt, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, Solution_time, Suppress_warning_about_MPI_COMM_WORLD, and oomph::TimingHelpers::timer().

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

◆ shutdown_mumps()

void oomph::MumpsSolver::shutdown_mumps ( )

private

Shutdown mumps.

   {
     if (Mumps_is_initialised)
     {
       // Shut down flag
       Mumps_struc_pt->job = -2;
  
       // Call c interface to double precision mumps to shut down
       dmumps_c(Mumps_struc_pt);
  
       // Cleanup
       delete Mumps_struc_pt;
       Mumps_struc_pt = 0;
  
       Mumps_is_initialised = false;
  
       // Cleanup storage
       Irn_loc.clear();
       Jcn_loc.clear();
       A_loc.clear();
     }
   }

References A_loc, Irn_loc, Jcn_loc, Mumps_is_initialised, and Mumps_struc_pt.

Referenced by clean_up_memory(), factorise(), and ~MumpsSolver().

◆ solve() [1/2]

void oomph::MumpsSolver::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. Matrix must be CRDoubleMatrix. 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.

   {
 #ifdef PARANOID
     if (!Suppress_warning_about_MPI_COMM_WORLD)
     {
       if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt)
             ->distribution_pt()
             ->communicator_pt()
             ->mpi_comm() != MPI_COMM_WORLD)
       {
         std::ostringstream error_message_stream;
         error_message_stream
           << "Warning: Mumps wrapper assumes that communicator is "
              "MPI_COMM_WORLD\n"
           << "         which is not the case, so mumps may die...\n"
           << "         If it does initialise oomph-lib's mpi without "
              "requesting\n"
           << "         the communicator to be a duplicate of MPI_COMM_WORLD\n"
           << "         (done via an optional boolean to "
              "MPI_Helpers::init(...)\n\n"
           << "         (You can suppress this warning by recompiling without\n"
           << "         paranoia or calling \n\n"
           << "         "
              "MumpsSolver::enable_suppress_warning_about_MPI_COMM_WORLD()\n\n"
           << "         \n";
         OomphLibWarning(error_message_stream.str(),
                         "MumpsSolver::solve()",
                         OOMPH_EXCEPTION_LOCATION);
       }
     }
 #endif
  
     // Initialise timer
     double t_start = TimingHelpers::timer();
  
 #ifdef PARANOID
     // check that the rhs vector is setup
     if (!rhs.distribution_pt()->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 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* ddist_matrix_pt =
       dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
     if (ddist_matrix_pt != 0)
     {
       if (!(*ddist_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
     DistributableLinearAlgebraObject* dist_matrix_pt =
       dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt);
     if (dist_matrix_pt)
     {
       // the solver has the same distribution as the matrix if possible
       this->build_distribution(dist_matrix_pt->distribution_pt());
     }
     else
     {
       // the solver has the same distribution as the RHS
       this->build_distribution(rhs.distribution_pt());
     }
  
     // Factorise the matrix
     factorise(matrix_pt);
  
     // Now do the back solve
     backsub(rhs, result);
  
     // Doc time for solve
     double t_end = TimingHelpers::timer();
     Solution_time = t_end - t_start;
  
     if ((Doc_time) &&
         (this->distribution_pt()->communicator_pt()->my_rank() == 0))
     {
       oomph_info << "Time for MumpsSolver solve                              : "
                  << TimingHelpers::convert_secs_to_formatted_string(t_end -
                                                                     t_start)
                  << std::endl;
     }
  
     // Switch off docing again by setting output stream for global info on
     // to negative number
     Mumps_struc_pt->ICNTL(3) = -1;
  
     // 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(), oomph::LinearAlgebraDistribution::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(), Mumps_struc_pt, oomph::DoubleMatrixBase::ncol(), oomph::DistributableLinearAlgebraObject::nrow(), oomph::DoubleMatrixBase::nrow(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, Solution_time, Suppress_warning_about_MPI_COMM_WORLD, and oomph::TimingHelpers::timer().

◆ solve() [2/2]

void oomph::MumpsSolver::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.

   {
 #ifdef PARANOID
     if (!Suppress_warning_about_MPI_COMM_WORLD)
     {
       if (problem_pt->communicator_pt()->mpi_comm() != MPI_COMM_WORLD)
       {
         std::ostringstream error_message_stream;
         error_message_stream
           << "Warning: Mumps wrapper assumes that communicator is "
              "MPI_COMM_WORLD\n"
           << "         which is not the case, so mumps may die...\n"
           << "         If it does initialise oomph-lib's mpi without "
              "requesting\n"
           << "         the communicator to be a duplicate of MPI_COMM_WORLD\n"
           << "         (done via an optional boolean to "
              "MPI_Helpers::init(...)\n\n"
           << "         (You can suppress this warning by recompiling without\n"
           << "         paranoia or calling \n\n"
           << "         "
              "MumpsSolver::enable_suppress_warning_about_MPI_COMM_WORLD()\n\n"
           << "         \n";
         OomphLibWarning(error_message_stream.str(),
                         "MumpsSolver::solve()",
                         OOMPH_EXCEPTION_LOCATION);
       }
     }
 #endif
  
     // Initialise timer
     double t_start = TimingHelpers::timer();
  
     // number of dofs
     unsigned n_dof = problem_pt->ndof();
  
     // Set the distribution for the solver.
     this->distribution_pt()->build(problem_pt->communicator_pt(), n_dof);
  
     // Take a copy of Delete_matrix_data
     bool copy_of_Delete_matrix_data = Delete_matrix_data;
  
     // Set Delete_matrix to true
     Delete_matrix_data = true;
  
     // Initialise timer
     t_start = TimingHelpers::timer();
  
     // Storage for the distributed 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;
     int my_rank = this->distribution_pt()->communicator_pt()->my_rank();
     if ((Doc_time) && (my_rank == 0))
     {
       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);
       }
     }
  
     // Set Delete_matrix back to original value
     Delete_matrix_data = copy_of_Delete_matrix_data;
  
     // Finalise/doc timings
     if ((Doc_time) && (my_rank == 0))
     {
       double t_end = TimingHelpers::timer();
       oomph_info << "Total time for MumpsSolver "
                  << "(np="
                  << this->distribution_pt()->communicator_pt()->nproc()
                  << ",N=" << problem_pt->ndof() << ")                : "
                  << TimingHelpers::convert_secs_to_formatted_string(t_end -
                                                                     t_start)
                  << std::endl;
     }
   }

References oomph::DoubleVector::build(), oomph::LinearAlgebraDistribution::build(), oomph::DoubleVector::built(), oomph::Problem::communicator_pt(), oomph::LinearAlgebraDistribution::communicator_pt(), oomph::TimingHelpers::convert_secs_to_formatted_string(), Delete_matrix_data, oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::LinearSolver::Doc_time, oomph::Problem::get_jacobian(), Jacobian_setup_time, oomph::OomphCommunicator::my_rank(), oomph::Problem::ndof(), oomph::OomphCommunicator::nproc(), OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::DoubleVector::redistribute(), Suppress_solve, Suppress_warning_about_MPI_COMM_WORLD, and oomph::TimingHelpers::timer().

◆ use_metis_ordering()

void oomph::MumpsSolver::use_metis_ordering ( )

inline

     {
       Jacobian_ordering_flag = Metis_ordering;
     }

References Jacobian_ordering_flag, and Metis_ordering.

◆ use_pord_ordering()

void oomph::MumpsSolver::use_pord_ordering ( )

inline

     {
       Jacobian_ordering_flag = Pord_ordering;
     }

References Jacobian_ordering_flag, and Pord_ordering.

◆ use_scotch_ordering()

void oomph::MumpsSolver::use_scotch_ordering ( )

inline

     {
       Jacobian_ordering_flag = Scotch_ordering;
     }

References Jacobian_ordering_flag, and Scotch_ordering.

Member Data Documentation

◆ A_loc

Vector<double> oomph::MumpsSolver::A_loc

private

Referenced by factorise(), and shutdown_mumps().

◆ Default_workspace_scaling_factor

int oomph::MumpsSolver::Default_workspace_scaling_factor = 5

static

Default factor for workspace – static so it can be overwritten globally.

/////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////// Default factor for workspace – static so it can be overwritten globally.

◆ Delete_matrix_data

bool oomph::MumpsSolver::Delete_matrix_data

private

Delete_matrix_data flag. MumpsSolver 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.

Referenced by disable_delete_matrix_data(), enable_delete_matrix_data(), factorise(), and solve().

◆ Doc_stats

bool oomph::MumpsSolver::Doc_stats

private

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

Referenced by backsub(), disable_doc_stats(), enable_doc_stats(), factorise(), and resolve().

◆ Irn_loc

Vector<int> oomph::MumpsSolver::Irn_loc

private

Vector for row numbers.

Referenced by factorise(), and shutdown_mumps().

◆ Jacobian_ordering_flag

unsigned oomph::MumpsSolver::Jacobian_ordering_flag

private

stores an integer from the public enum which specifies which package (PORD, Metis or SCOTCH) is used to reorder the Jacobian matrix during the analysis

Referenced by initialise_mumps(), use_metis_ordering(), use_pord_ordering(), and use_scotch_ordering().

◆ Jacobian_setup_time

double oomph::MumpsSolver::Jacobian_setup_time

private

Jacobian setup time.

Referenced by jacobian_setup_time(), and solve().

◆ Jacobian_symmetry_flag

unsigned oomph::MumpsSolver::Jacobian_symmetry_flag

private

symmetry of the Jacobian matrix we're solving; takes one of the enum values above

Referenced by declare_jacobian_is_symmetric(), declare_jacobian_is_symmetric_positive_definite(), declare_jacobian_is_unsymmetric(), and initialise_mumps().

◆ Jcn_loc

Vector<int> oomph::MumpsSolver::Jcn_loc

private

Referenced by factorise(), and shutdown_mumps().

◆ Mumps_is_initialised

bool oomph::MumpsSolver::Mumps_is_initialised

private

Has mumps been initialised?

Referenced by backsub(), factorise(), initialise_mumps(), and shutdown_mumps().

◆ Mumps_struc_pt

DMUMPS_STRUC_C* oomph::MumpsSolver::Mumps_struc_pt

private

Pointer to MUMPS struct that contains the solver data.

Referenced by backsub(), factorise(), initialise_mumps(), resolve(), shutdown_mumps(), and solve().

◆ Solution_time

double oomph::MumpsSolver::Solution_time

private

Solution time.

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

◆ Suppress_incorrect_rhs_distribution_warning_in_resolve

bool oomph::MumpsSolver::Suppress_incorrect_rhs_distribution_warning_in_resolve

static

Initial value:

=

false

Static flag that determines whether the warning about incorrect distribution of RHSs will be printed or not

Static warning to suppress warnings about incorrect distribution of RHS vector. Default is false

Referenced by backsub().

◆ Suppress_mumps_info_during_solve

bool oomph::MumpsSolver::Suppress_mumps_info_during_solve

private

Boolean to suppress info being printed to screen during MUMPS solve.

Referenced by disable_suppress_mumps_info_during_solve(), enable_suppress_mumps_info_during_solve(), and factorise().

◆ Suppress_solve

bool oomph::MumpsSolver::Suppress_solve

private

Suppress solve?

Referenced by disable_suppress_solve(), enable_suppress_solve(), and solve().

◆ Suppress_warning_about_MPI_COMM_WORLD

bool oomph::MumpsSolver::Suppress_warning_about_MPI_COMM_WORLD

private

Boolean to suppress warning about communicator not equal to MPI_COMM_WORLD

Referenced by backsub(), disable_suppress_warning_about_MPI_COMM_WORLD(), enable_suppress_warning_about_MPI_COMM_WORLD(), factorise(), resolve(), and solve().

◆ Workspace_scaling_factor

unsigned oomph::MumpsSolver::Workspace_scaling_factor

private

Referenced by factorise().

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

Public Member Functions

Static Public Attributes

Private Types

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Member Enumeration Documentation

◆ MumpsJacobianOrderingFlags

◆ MumpsJacobianSymmetryFlags

Constructor & Destructor Documentation

◆ MumpsSolver() [1/2]

◆ MumpsSolver() [2/2]

◆ ~MumpsSolver()

Member Function Documentation

◆ backsub()

◆ clean_up_memory()

◆ declare_jacobian_is_symmetric()

◆ declare_jacobian_is_symmetric_positive_definite()

◆ declare_jacobian_is_unsymmetric()

◆ disable_delete_matrix_data()

◆ disable_doc_stats()

◆ disable_resolve()

◆ disable_suppress_mumps_info_during_solve()

◆ disable_suppress_solve()

◆ disable_suppress_warning_about_MPI_COMM_WORLD()

◆ enable_delete_matrix_data()

◆ enable_doc_stats()

◆ enable_suppress_mumps_info_during_solve()

◆ enable_suppress_solve()

◆ enable_suppress_warning_about_MPI_COMM_WORLD()

◆ factorise()

◆ initialise_mumps()

◆ jacobian_setup_time()

◆ linear_solver_solution_time()

◆ operator=()

◆ resolve()

◆ shutdown_mumps()

◆ solve() [1/2]

◆ solve() [2/2]

◆ use_metis_ordering()

◆ use_pord_ordering()

◆ use_scotch_ordering()

Member Data Documentation

◆ A_loc

◆ Default_workspace_scaling_factor

◆ Delete_matrix_data

◆ Doc_stats

◆ Irn_loc

◆ Jacobian_ordering_flag

◆ Jacobian_setup_time

◆ Jacobian_symmetry_flag

◆ Jcn_loc

◆ Mumps_is_initialised

◆ Mumps_struc_pt

◆ Solution_time

◆ Suppress_incorrect_rhs_distribution_warning_in_resolve

◆ Suppress_mumps_info_during_solve

◆ Suppress_solve

◆ Suppress_warning_about_MPI_COMM_WORLD

◆ Workspace_scaling_factor