oomph::BlockPitchForkLinearSolver Class Reference

#include <assembly_handler.h>

Inheritance diagram for oomph::BlockPitchForkLinearSolver:

Public Member Functions
	BlockPitchForkLinearSolver (LinearSolver *const linear_solver_pt)
	Constructor, inherits the original linear solver. More...
	~BlockPitchForkLinearSolver ()
	Destructor: clean up the allocated memory. More...
void	solve (Problem *const &problem_pt, DoubleVector &result)
	The solve function uses the block factorisation. More...
void	solve (DoubleMatrixBase *const &matrix_pt, const DoubleVector &rhs, DoubleVector &result)
void	solve (DoubleMatrixBase *const &matrix_pt, const Vector< double > &rhs, Vector< double > &result)
void	resolve (const DoubleVector &rhs, DoubleVector &result)
	The resolve function also uses the block factorisation. More...
LinearSolver *	linear_solver_pt () const
	Access function to the original linear solver. 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	disable_resolve ()
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 void	clean_up_memory ()
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)

Private Attributes
LinearSolver *	Linear_solver_pt
	Pointer to the original linear solver. More...
Problem *	Problem_pt
	Pointer to the problem, used in the resolve. More...
DoubleVector *	B_pt
	Pointer to the storage for the vector b. More...
DoubleVector *	C_pt
	Pointer to the storage for the vector c. More...
DoubleVector *	D_pt
	Pointer to the storage for the vector d. More...
DoubleVector *	dJy_dparam_pt

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

A custom linear solver class that is used to solve a block-factorised version of the PitchFork bifurcation detection problem.

Constructor & Destructor Documentation

BlockPitchForkLinearSolver()

oomph::BlockPitchForkLinearSolver::BlockPitchForkLinearSolver ( LinearSolver *const  linear_solver_pt )

inline

Constructor, inherits the original linear solver.

      : Linear_solver_pt(linear_solver_pt),
        Problem_pt(0),
        B_pt(0),
        C_pt(0),
        D_pt(0),
        dJy_dparam_pt(0)
    {
    }

~BlockPitchForkLinearSolver()

oomph::BlockPitchForkLinearSolver::~BlockPitchForkLinearSolver

(

)

Destructor: clean up the allocated memory.

Clean up the memory that may have been allocated by the solver.

  {
    if (B_pt != 0)
    {
      delete B_pt;
    }
    if (C_pt != 0)
    {
      delete C_pt;
    }
    if (D_pt != 0)
    {
      delete D_pt;
    }
    if (dJy_dparam_pt != 0)
    {
      delete dJy_dparam_pt;
    }
  }

References B_pt, C_pt, D_pt, and dJy_dparam_pt.

Member Function Documentation

linear_solver_pt()

LinearSolver* oomph::BlockPitchForkLinearSolver::linear_solver_pt ( ) const

inline

Access function to the original linear solver.

    {
      return Linear_solver_pt;
    }

References Linear_solver_pt.

Referenced by oomph::PitchForkHandler::~PitchForkHandler().

resolve()

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

virtual

The resolve function also uses the block factorisation.

Reimplemented from oomph::LinearSolver.

  {
    std::cout << "Block pitchfork resolve" << std::endl;
    // Check that the factors have been stored
    if (B_pt == 0)
    {
      throw OomphLibError("The required vectors have not been stored",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
 
    // Cache pointer to the problem
    Problem* const problem_pt = Problem_pt;
 
    // Locally cache pointer to the handler
    PitchForkHandler* handler_pt =
      static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());
 
    // Get the augmented distribution from the handler
    LinearAlgebraDistribution aug_dist =
      handler_pt->Augmented_dof_distribution_pt;
    // this->build_distribution(aug_dist);
 
    // Find the number of dofs of the augmented system
    // const unsigned n_aug_dof = problem_pt->ndof();
 
    // Storage for the result distribution
    LinearAlgebraDistribution result_dist;
 
    // if the result vector is not setup then rebuild with distribution
    // = natural distribution of augmented solver
    if (!result.built())
    {
      result.build(&aug_dist, 0.0);
    }
    // Otherwise store the incoming distribution and redistribute
    else
    {
      result_dist.build(result.distribution_pt());
      result.redistribute(&aug_dist);
    }
 
 
    // Locally cache a pointer to the parameter
    // double* const parameter_pt = handler_pt->Parameter_pt;
 
    // Switch things to our block solver
    handler_pt->solve_block_system();
    // We need to find out the number of dofs
    // unsigned n_dof = problem_pt->ndof();
 
    // create the linear algebra distribution for this solver
    // currently only global (non-distributed) distributions are allowed
    // LinearAlgebraDistribution
    // dist(problem_pt->communicator_pt(),n_dof,false);
    // this->build_distribution(dist);
 
    // if the result vector is not setup then rebuild with distribution = global
    // if (!result.built())
    // {
    //  result.build(this->distribution_pt(),0.0);
    // }
 
 
    // Copy the rhs into local storage
    DoubleVector rhs_local = rhs;
    // and redistribute into the augmented distribution
    rhs_local.redistribute(&aug_dist);
 
    // Setup storage for output
    DoubleVector x1(Linear_solver_pt->distribution_pt(), 0.0);
    DoubleVector x3(Linear_solver_pt->distribution_pt(), 0.0);
 
    // Create input for RHS with the natural distribution
    DoubleVector input_x1(handler_pt->Dof_distribution_pt);
    const unsigned n_dof_local = input_x1.nrow_local();
 
    // Set the values of the a vector
    for (unsigned n = 0; n < n_dof_local; n++)
    {
      input_x1[n] = rhs_local[n];
    }
    // Need to redistribute into the linear algebra distribution
    input_x1.redistribute(Linear_solver_pt->distribution_pt());
 
    // We want to resolve, of course
    Linear_solver_pt->enable_resolve();
    // Now solve for the first vector
    Linear_solver_pt->resolve(input_x1, x1);
 
    // Get the symmetry vector from the handler
    DoubleVector psi = handler_pt->Psi;
    // redistribute local copy
    psi.redistribute(Linear_solver_pt->distribution_pt());
 
    // We can now construct various dot products
    double psi_d = psi.dot(*D_pt);
    double psi_c = psi.dot(*C_pt);
    double psi_x1 = psi.dot(x1);
 
    // Calculate another intermediate constant
    const double Psi = psi_d / psi_c;
 
    // Construct the second intermediate value,
    // assumption is that rhs has been set to the current value of the residuals
    double parameter_residual = rhs_local[n_dof_local];
#ifdef OOMPH_HAS_MPI
    // Broadcast to all others, if we have a distributed problem
    if (problem_pt->distributed())
    {
      MPI_Bcast(&parameter_residual,
                1,
                MPI_DOUBLE,
                0,
                problem_pt->communicator_pt()->mpi_comm());
    }
#endif
 
    double x2 = (psi_x1 - parameter_residual) / psi_c;
 
    // Now construct the vectors that multiply the jacobian terms
    // Vector<double> X1(n_dof/*+1*/);
    Vector<DoubleVectorWithHaloEntries> X1(1);
    X1[0].build(Linear_solver_pt->distribution_pt(), 0.0);
 
    const unsigned n_dof_local_linear_solver =
      Linear_solver_pt->distribution_pt()->nrow_local();
    for (unsigned n = 0; n < n_dof_local_linear_solver; n++)
    {
      X1[0][n] = x1[n] - x2 * (*C_pt)[n];
    }
 
    // Local storage for the product term
    Vector<DoubleVectorWithHaloEntries> Jprod_X1(1);
 
    // Redistribute the inputs to have the Dof distribution pt
    X1[0].redistribute(handler_pt->Dof_distribution_pt);
 
    // Set up the halo scheme
#ifdef OOMPH_HAS_MPI
    X1[0].build_halo_scheme(problem_pt->Halo_scheme_pt);
#endif
 
    // Get the product from the problem
    problem_pt->get_hessian_vector_products(handler_pt->Y, X1, Jprod_X1);
 
    // In standard cases the offset here will be one
    unsigned offset = 1;
#ifdef OOMPH_HAS_MPI
    // If we are distributed and not on the first processor
    // then there is no offset and the eigenfunction is
    // directly after the standard dofs
    int my_rank = problem_pt->communicator_pt()->my_rank();
    if ((problem_pt->distributed()) && (my_rank != 0))
    {
      offset = 0;
    }
#endif
 
    // OK, now we can formulate the next vectors
    //(again assuming result contains residuals)
    // Local storage for the product terms
    DoubleVector Mod_Jprod_X1(handler_pt->Dof_distribution_pt, 0.0);
 
    for (unsigned n = 0; n < n_dof_local; n++)
    {
      Mod_Jprod_X1[n] = rhs_local[n_dof_local + offset + n] - Jprod_X1[0][n] -
                        x2 * (*dJy_dparam_pt)[n];
    }
 
    // Redistribute back to the linear solver distribution
    Mod_Jprod_X1.redistribute(Linear_solver_pt->distribution_pt());
 
    // Liner solve to get x3
    Linear_solver_pt->resolve(Mod_Jprod_X1, x3);
 
    // Construst a couple of additional products
    double l_x3 = psi.dot(x3);
    double l_b = psi.dot(*B_pt);
 
    // get the last intermediate variable
    // The required parameter is that corresponding to the dof, which
    // is only stored on the root processor
    double sigma_residual = rhs_local[2 * (n_dof_local + offset) - 1];
#ifdef OOMPH_HAS_MPI
    // Broadcast to all others, if we have a distributed problem
    if (problem_pt->distributed())
    {
      MPI_Bcast(&sigma_residual,
                1,
                MPI_DOUBLE,
                0,
                problem_pt->communicator_pt()->mpi_comm());
    }
#endif
 
    // get the last intermediate variable
    const double delta_sigma = (l_x3 - sigma_residual) / l_b;
    const double delta_lambda = x2 - delta_sigma * Psi;
 
    // Create temporary DoubleVectors to hold the results
    DoubleVector res1(Linear_solver_pt->distribution_pt());
    DoubleVector res2(Linear_solver_pt->distribution_pt());
 
    for (unsigned n = 0; n < n_dof_local_linear_solver; n++)
    {
      res1[n] = x1[n] - delta_lambda * (*C_pt)[n] - delta_sigma * (*D_pt)[n];
      res2[n] = x3[n] - delta_sigma * (*B_pt)[n];
    }
 
    // Now redistribute these into the Dof distribution pointer
    res1.redistribute(handler_pt->Dof_distribution_pt);
    res2.redistribute(handler_pt->Dof_distribution_pt);
 
    // Now can copy over into results into the result vector
    for (unsigned n = 0; n < n_dof_local; n++)
    {
      result[n] = res1[n];
      result[n_dof_local + offset + n] = res2[n];
    }
 
    // Add the final contributions to the residuals
    // only on the root processor if we have a distributed problem
#ifdef OOMPH_HAS_MPI
    if ((!problem_pt->distributed()) || (my_rank == 0))
#endif
    {
      result[n_dof_local] = delta_lambda;
      result[2 * n_dof_local + 1] = delta_sigma;
    }
 
    // Redistribute the result into its incoming distribution (if it had one)
    if (result_dist.built())
    {
      result.redistribute(&result_dist);
    }
 
    Linear_solver_pt->disable_resolve();
 
    // Switch things to our block solver
    handler_pt->solve_full_system();
  }

References oomph::Problem::assembly_handler_pt(), oomph::PitchForkHandler::Augmented_dof_distribution_pt, B_pt, oomph::DoubleVector::build(), oomph::LinearAlgebraDistribution::build(), oomph::DoubleVector::built(), oomph::LinearAlgebraDistribution::built(), C_pt, oomph::Problem::communicator_pt(), D_pt, oomph::LinearSolver::disable_resolve(), oomph::Problem::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::PitchForkHandler::Dof_distribution_pt, oomph::DoubleVector::dot(), oomph::LinearSolver::enable_resolve(), oomph::Problem::get_hessian_vector_products(), Linear_solver_pt, oomph::OomphCommunicator::my_rank(), n, oomph::LinearAlgebraDistribution::nrow_local(), oomph::DistributableLinearAlgebraObject::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, Problem_pt, oomph::PitchForkHandler::Psi, oomph::DoubleVector::redistribute(), oomph::LinearSolver::resolve(), oomph::PitchForkHandler::solve_block_system(), oomph::PitchForkHandler::solve_full_system(), Global_parameters::x1(), Global_parameters::x2(), and oomph::PitchForkHandler::Y.

solve() [1/3]

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

inlinevirtual

The linear-algebra-type solver does not make sense. The interface is deliberately broken

Reimplemented from oomph::LinearSolver.

    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

solve() [2/3]

void oomph::BlockPitchForkLinearSolver::solve ( DoubleMatrixBase *const &  matrix_pt, const Vector< double > &  rhs, Vector< double > &  result  )

inlinevirtual

The linear-algebra-type solver does not make sense. The interface is deliberately broken

Reimplemented from oomph::LinearSolver.

    {
      throw OomphLibError(
        "Linear-algebra interface does not make sense for this linear solver\n",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }

References OOMPH_CURRENT_FUNCTION, and OOMPH_EXCEPTION_LOCATION.

solve() [3/3]

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

virtual

The solve function uses the block factorisation.

Use a block factorisation to solve the augmented system associated with a PitchFork bifurcation.

Implements oomph::LinearSolver.

  {
    std::cout << "Block pitchfork solve" << std::endl;
    // Locally cache the pointer to the handler.
    PitchForkHandler* handler_pt =
      static_cast<PitchForkHandler*>(problem_pt->assembly_handler_pt());
 
    // Get the augmented distribution from the handler and use it as
    // the distribution of the linear solver
    LinearAlgebraDistribution aug_dist =
      handler_pt->Augmented_dof_distribution_pt;
    this->build_distribution(aug_dist);
 
    // If the result vector is not setup then die
    if (!result.built())
    {
      throw OomphLibError("Result vector must be built\n",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    // Store the distribution of the result, which is probably not in
    // the natural distribution of the augmented solver
    LinearAlgebraDistribution result_dist(result.distribution_pt());
 
    // Locally cache a pointer to the parameter
    double* const parameter_pt = handler_pt->Parameter_pt;
 
    // Firstly get the derivatives of the full augmented system with
    // respect to the parameter
 
    // Allocate storage for the derivatives of the residuals with respect
    // to the global parameter using the augmented distribution
    DoubleVector dRdparam;
    // Then get the appropriate derivatives which will come back in
    // some distribution
    problem_pt->get_derivative_wrt_global_parameter(parameter_pt, dRdparam);
    // Redistribute into the augmented distribution
    dRdparam.redistribute(&aug_dist);
 
    // Switch the handler to "block solver" mode (sort out distribution)
    handler_pt->solve_block_system();
 
    // Temporary vector for the result (I shouldn't have to set this up)
    DoubleVector x1;
 
    // We are going to do resolves using the underlying linear solver
    Linear_solver_pt->enable_resolve();
    // Solve the first (standard) system Jx1 = R
    Linear_solver_pt->solve(problem_pt, x1);
 
    // Allocate storage for B, C and D which can be used in the resolve
    // and the derivatives of the jacobian/eigenvector product with
    // respect to the parameter
    if (B_pt != 0)
    {
      delete B_pt;
    }
    B_pt = new DoubleVector(Linear_solver_pt->distribution_pt(), 0.0);
    if (C_pt != 0)
    {
      delete C_pt;
    }
    C_pt = new DoubleVector(Linear_solver_pt->distribution_pt(), 0.0);
    if (D_pt != 0)
    {
      delete D_pt;
    }
    D_pt = new DoubleVector(Linear_solver_pt->distribution_pt(), 0.0);
    // Need this to be in the distribution of the dofs
    if (dJy_dparam_pt != 0)
    {
      delete dJy_dparam_pt;
    }
    dJy_dparam_pt = new DoubleVector(handler_pt->Dof_distribution_pt, 0.0);
 
    // Get the symmetry vector from the handler should have the
    // Dof distribution
    DoubleVector psi = handler_pt->Psi;
 
    // Temporary vector for the rhs that is dR/dparam (augmented distribution)
    // DoubleVector F(Linear_solver_pt->distribution_pt(),0.0);
    DoubleVector F(handler_pt->Dof_distribution_pt);
    const unsigned n_dof_local = F.nrow_local();
 
    // f.nrow local copied from dRdparam
    for (unsigned n = 0; n < n_dof_local; n++)
    {
      F[n] = dRdparam[n];
    }
    // Fill in the rhs that is dJy/dparam  //dJy_dparam nrow local
 
    // In standard cases the offset here will be one
    unsigned offset = 1;
#ifdef OOMPH_HAS_MPI
    // If we are distributed and not on the first processor
    // then there is no offset and the eigenfunction is
    // directly after the standard dofs
    int my_rank = problem_pt->communicator_pt()->my_rank();
    if ((problem_pt->distributed()) && (my_rank != 0))
    {
      offset = 0;
    }
#endif
    for (unsigned n = 0; n < n_dof_local; n++)
    {
      (*dJy_dparam_pt)[n] = dRdparam[n_dof_local + offset + n];
    }
 
    // Now resolve to find c and d
    // First, redistribute F and psi into the linear algebra distribution
    F.redistribute(Linear_solver_pt->distribution_pt());
    psi.redistribute(Linear_solver_pt->distribution_pt());
 
    Linear_solver_pt->resolve(F, *C_pt);
    Linear_solver_pt->resolve(psi, *D_pt);
 
    // We can now construct various dot products
    double psi_d = psi.dot(*D_pt);
    double psi_c = psi.dot(*C_pt);
    double psi_x1 = psi.dot(x1);
 
    // Calculate another intermediate constant
    const double Psi = psi_d / psi_c;
 
    // Construct the second intermediate value,
    // assumption is that result has been
    // set to the current value of the residuals
    // First, redistribute into the Natural distribution of the augmented system
    result.redistribute(&aug_dist);
    // The required parameter is that corresponding to the dof, which
    // is only stored on the root processor
    double parameter_residual = result[n_dof_local];
#ifdef OOMPH_HAS_MPI
    // Broadcast to all others, if we have a distributed problem
    if (problem_pt->distributed())
    {
      MPI_Bcast(&parameter_residual,
                1,
                MPI_DOUBLE,
                0,
                problem_pt->communicator_pt()->mpi_comm());
    }
#endif
    // Now construct the value
    double x2 = (psi_x1 - parameter_residual) / psi_c;
 
    // Now construct the vectors that multiply the jacobian terms
    Vector<DoubleVectorWithHaloEntries> D_and_X1(2);
    D_and_X1[0].build(Linear_solver_pt->distribution_pt(), 0.0);
    D_and_X1[1].build(Linear_solver_pt->distribution_pt(), 0.0);
    // Fill in the appropriate terms
    // Get the number of local dofs from the Linear_solver_pt distribution
    const unsigned n_dof_local_linear_solver =
      Linear_solver_pt->distribution_pt()->nrow_local();
    for (unsigned n = 0; n < n_dof_local_linear_solver; n++)
    {
      const double C_ = (*C_pt)[n];
      D_and_X1[0][n] = (*D_pt)[n] - Psi * C_;
      D_and_X1[1][n] = x1[n] - x2 * C_;
    }
 
    // Local storage for the result of the product terms
    Vector<DoubleVectorWithHaloEntries> Jprod_D_and_X1(2);
 
    // Redistribute the inputs to have the Dof distribution pt
    D_and_X1[0].redistribute(handler_pt->Dof_distribution_pt);
    D_and_X1[1].redistribute(handler_pt->Dof_distribution_pt);
 
    // Set up the halo scheme
#ifdef OOMPH_HAS_MPI
    D_and_X1[0].build_halo_scheme(problem_pt->Halo_scheme_pt);
    D_and_X1[1].build_halo_scheme(problem_pt->Halo_scheme_pt);
#endif
 
    // Get the products from the new problem function
    problem_pt->get_hessian_vector_products(
      handler_pt->Y, D_and_X1, Jprod_D_and_X1);
 
    // OK, now we can formulate the next vectors
    //(again assuming result contains residuals)
    // Need to redistribute F to the dof distribution
    // F.redistribute(handler_pt->Dof_distribution_pt);
    DoubleVector G(handler_pt->Dof_distribution_pt);
 
    for (unsigned n = 0; n < n_dof_local; n++)
    {
      G[n] = result[n_dof_local + offset + n] - Jprod_D_and_X1[1][n] -
             x2 * dRdparam[n_dof_local + offset + n];
      Jprod_D_and_X1[0][n] *= -1.0;
      Jprod_D_and_X1[0][n] -= Psi * dRdparam[n_dof_local + offset + n];
    }
 
 
    // Then redistribute back to the linear solver distribution
    G.redistribute(Linear_solver_pt->distribution_pt());
    Jprod_D_and_X1[0].redistribute(Linear_solver_pt->distribution_pt());
    Jprod_D_and_X1[1].redistribute(Linear_solver_pt->distribution_pt());
 
    // Linear solve to get B
    Linear_solver_pt->resolve(Jprod_D_and_X1[0], *B_pt);
    // Liner solve to get x3
    DoubleVector x3(Linear_solver_pt->distribution_pt(), 0.0);
    Linear_solver_pt->resolve(G, x3);
 
    // Construst a couple of additional products
    double l_x3 = psi.dot(x3);
    double l_b = psi.dot(*B_pt);
 
    // get the last intermediate variable
    // The required parameter is that corresponding to the dof, which
    // is only stored on the root processor
    double sigma_residual = result[2 * (n_dof_local + offset) - 1];
#ifdef OOMPH_HAS_MPI
    // Broadcast to all others, if we have a distributed problem
    if (problem_pt->distributed())
    {
      MPI_Bcast(&sigma_residual,
                1,
                MPI_DOUBLE,
                0,
                problem_pt->communicator_pt()->mpi_comm());
    }
#endif
 
 
    const double delta_sigma = (l_x3 - sigma_residual) / l_b;
    const double delta_lambda = x2 - delta_sigma * Psi;
 
    // Need to do some more rearrangements here because result is global
    // but the other vectors are not!
 
    // Create temporary DoubleVectors to hold the results
    DoubleVector res1(Linear_solver_pt->distribution_pt());
    DoubleVector res2(Linear_solver_pt->distribution_pt());
 
    for (unsigned n = 0; n < n_dof_local_linear_solver; n++)
    {
      res1[n] = x1[n] - delta_lambda * (*C_pt)[n] - delta_sigma * (*D_pt)[n];
      res2[n] = x3[n] - delta_sigma * (*B_pt)[n];
    }
 
    // Now redistribute these into the Dof distribution pointer
    res1.redistribute(handler_pt->Dof_distribution_pt);
    res2.redistribute(handler_pt->Dof_distribution_pt);
 
    // Now can copy over into results into the result vector
    for (unsigned n = 0; n < n_dof_local; n++)
    {
      result[n] = res1[n];
      result[n_dof_local + offset + n] = res2[n];
    }
 
    // Add the final contributions to the residuals
    // only on the root processor if we have a distributed problem
#ifdef OOMPH_HAS_MPI
    if ((!problem_pt->distributed()) || (my_rank == 0))
#endif
    {
      result[n_dof_local] = delta_lambda;
      result[2 * n_dof_local + 1] = delta_sigma;
    }
 
 
    // The sign of the determinant is given by the sign of
    // the product psi_c and l_b
    // NOT CHECKED YET!
    problem_pt->sign_of_jacobian() =
      static_cast<int>(std::fabs(psi_c * l_b) / (psi_c * l_b));
 
    // Redistribute the result into its incoming distribution
    result.redistribute(&result_dist);
 
    // Switch things to our block solver
    handler_pt->solve_full_system();
 
    // If we are not storing things, clear the results
    if (!Enable_resolve)
    {
      // We no longer need to store the matrix
      Linear_solver_pt->disable_resolve();
      delete B_pt;
      B_pt = 0;
      delete C_pt;
      C_pt = 0;
      delete D_pt;
      D_pt = 0;
      delete dJy_dparam_pt;
      dJy_dparam_pt = 0;
    }
    // Otherwise also store the pointer to the problem
    else
    {
      Problem_pt = problem_pt;
    }
  }

References oomph::Problem::assembly_handler_pt(), oomph::PitchForkHandler::Augmented_dof_distribution_pt, B_pt, oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DoubleVector::built(), C_pt, oomph::Problem::communicator_pt(), D_pt, oomph::LinearSolver::disable_resolve(), oomph::Problem::distributed(), oomph::DistributableLinearAlgebraObject::distribution_pt(), dJy_dparam_pt, oomph::PitchForkHandler::Dof_distribution_pt, oomph::DoubleVector::dot(), oomph::LinearSolver::Enable_resolve, oomph::LinearSolver::enable_resolve(), oomph::OcTreeNames::F, boost::multiprecision::fabs(), G, oomph::Problem::get_derivative_wrt_global_parameter(), oomph::Problem::get_hessian_vector_products(), Linear_solver_pt, oomph::OomphCommunicator::my_rank(), n, oomph::LinearAlgebraDistribution::nrow_local(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::PitchForkHandler::Parameter_pt, Problem_pt, oomph::PitchForkHandler::Psi, oomph::DoubleVector::redistribute(), oomph::LinearSolver::resolve(), oomph::Problem::sign_of_jacobian(), oomph::LinearSolver::solve(), oomph::PitchForkHandler::solve_block_system(), oomph::PitchForkHandler::solve_full_system(), Global_parameters::x1(), Global_parameters::x2(), and oomph::PitchForkHandler::Y.

Member Data Documentation

B_pt

DoubleVector* oomph::BlockPitchForkLinearSolver::B_pt

private

Pointer to the storage for the vector b.

Referenced by resolve(), solve(), and ~BlockPitchForkLinearSolver().

C_pt

DoubleVector* oomph::BlockPitchForkLinearSolver::C_pt

private

Pointer to the storage for the vector c.

Referenced by resolve(), solve(), and ~BlockPitchForkLinearSolver().

D_pt

DoubleVector* oomph::BlockPitchForkLinearSolver::D_pt

private

Pointer to the storage for the vector d.

Referenced by resolve(), solve(), and ~BlockPitchForkLinearSolver().

dJy_dparam_pt

DoubleVector* oomph::BlockPitchForkLinearSolver::dJy_dparam_pt

private

Pointer to the storage for the vector of derivatives with respect to the bifurcation parameter

Referenced by solve(), and ~BlockPitchForkLinearSolver().

Linear_solver_pt

LinearSolver* oomph::BlockPitchForkLinearSolver::Linear_solver_pt

private

Pointer to the original linear solver.

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

Problem_pt

Problem* oomph::BlockPitchForkLinearSolver::Problem_pt

private

Pointer to the problem, used in the resolve.

Referenced by resolve(), and solve().

---

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

• assembly_handler.h
• assembly_handler.cc