#include <biharmonic_preconditioner.h>

Inheritance diagram for oomph::InexactSubBiharmonicPreconditioner:

Public Member Functions
	InexactSubBiharmonicPreconditioner (BiharmonicPreconditioner *master_prec_pt, const bool use_amg)

	~InexactSubBiharmonicPreconditioner ()
	destructor - just calls this->clean_up_memory() More...

void	clean_up_memory ()
	clean the memory More...

	InexactSubBiharmonicPreconditioner (const InexactSubBiharmonicPreconditioner &)=delete
	Broken copy constructor. More...

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

void	setup ()
	Setup the preconditioner. More...

void	preconditioner_solve (const DoubleVector &r, DoubleVector &z)
	Apply preconditioner to r. More...

void	compute_inexact_schur_complement ()

Public Member Functions inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
	BlockPreconditioner ()
	Constructor. More...

	BlockPreconditioner (const BlockPreconditioner &)=delete
	Broken copy constructor. More...

virtual	~BlockPreconditioner ()
	Destructor. More...

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

CRDoubleMatrix *	matrix_pt () const

void	turn_on_recursive_debug_flag ()

void	turn_off_recursive_debug_flag ()

void	turn_on_debug_flag ()
	Toggles on the debug flag. More...

void	turn_off_debug_flag ()
	Toggles off the debug flag. More...

void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)

void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse, const Vector< Vector< unsigned >> &doftype_coarsen_map_coarse)

virtual void	block_setup ()

void	block_setup (const Vector< unsigned > &dof_to_block_map)

void	get_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_matrix, const bool &ignore_replacement_block=false) const

CRDoubleMatrix	get_block (const unsigned &i, const unsigned &j, const bool &ignore_replacement_block=false) const

void	set_master_matrix_pt (CRDoubleMatrix *in_matrix_pt)
	Set the matrix_pt in the upper-most master preconditioner. More...

void	get_block_other_matrix (const unsigned &i, const unsigned &j, CRDoubleMatrix *in_matrix_pt, CRDoubleMatrix &output_matrix)

void	get_blocks (DenseMatrix< bool > &required_blocks, DenseMatrix< CRDoubleMatrix * > &block_matrix_pt) const

void	get_dof_level_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_block, const bool &ignore_replacement_block=false) const

void	get_dof_level_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix &output_block, const bool &ignore_replacement_block) const

CRDoubleMatrix	get_concatenated_block (const VectorMatrix< BlockSelector > &selected_block)

void	get_concatenated_block_vector (const Vector< unsigned > &block_vec_number, const DoubleVector &v, DoubleVector &b)

void	return_concatenated_block_vector (const Vector< unsigned > &block_vec_number, const DoubleVector &b, DoubleVector &v) const
	Takes concatenated block ordered vector, b, and copies its. More...

void	get_block_vectors (const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const

void	get_block_vectors (const DoubleVector &v, Vector< DoubleVector > &s) const

void	return_block_vectors (const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const

void	return_block_vectors (const Vector< DoubleVector > &s, DoubleVector &v) const

void	get_block_vector (const unsigned &n, const DoubleVector &v, DoubleVector &b) const

void	return_block_vector (const unsigned &n, const DoubleVector &b, DoubleVector &v) const

void	get_block_ordered_preconditioner_vector (const DoubleVector &v, DoubleVector &w)

void	return_block_ordered_preconditioner_vector (const DoubleVector &w, DoubleVector &v) const

unsigned	nblock_types () const
	Return the number of block types. More...

unsigned	ndof_types () const
	Return the total number of DOF types. More...

const Mesh *	mesh_pt (const unsigned &i) const

unsigned	nmesh () const

int	block_number (const unsigned &i_dof) const
	Return the block number corresponding to a global index i_dof. More...

int	index_in_block (const unsigned &i_dof) const

const LinearAlgebraDistribution *	block_distribution_pt (const unsigned &b) const
	Access function to the block distributions (const version). More...

LinearAlgebraDistribution *	block_distribution_pt (const unsigned b)
	Access function to the block distributions (non-const version). More...

LinearAlgebraDistribution *	dof_block_distribution_pt (const unsigned &b)
	Access function to the dof-level block distributions. More...

const LinearAlgebraDistribution *	master_distribution_pt () const

unsigned	ndof_types_in_mesh (const unsigned &i) const

bool	is_subsidiary_block_preconditioner () const

bool	is_master_block_preconditioner () const

void	set_block_output_to_files (const std::string &basefilename)

void	disable_block_output_to_files ()
	Turn off output of blocks (by clearing the basefilename string). More...

bool	block_output_on () const
	Test if output of blocks is on or not. More...

void	output_blocks_to_files (const std::string &basefilename, const unsigned &precision=8) const

void	post_block_matrix_assembly_partial_clear ()

BlockPreconditioner< CRDoubleMatrix > *	master_block_preconditioner_pt () const
	Access function to the master block preconditioner pt. More...

void	clear_block_preconditioner_base ()

void	document ()

Vector< Vector< unsigned > >	doftype_coarsen_map_fine () const

Vector< unsigned >	get_fine_grain_dof_types_in (const unsigned &i) const

unsigned	nfine_grain_dof_types_in (const unsigned &i) const

MapMatrix< unsigned, CRDoubleMatrix * >	replacement_dof_block_pt () const
	Access function to the replaced dof-level blocks. More...

void	setup_matrix_vector_product (MatrixVectorProduct matvec_prod_pt, CRDoubleMatrix block_pt, const Vector< unsigned > &block_col_indices)

void	setup_matrix_vector_product (MatrixVectorProduct matvec_prod_pt, CRDoubleMatrix block_pt, const unsigned &block_col_index)

void	internal_get_block_ordered_preconditioner_vector (const DoubleVector &v, DoubleVector &w) const

void	internal_return_block_ordered_preconditioner_vector (const DoubleVector &w, DoubleVector &v) const

unsigned	internal_nblock_types () const

unsigned	internal_ndof_types () const

void	internal_return_block_vector (const unsigned &n, const DoubleVector &b, DoubleVector &v) const

void	internal_get_block_vector (const unsigned &n, const DoubleVector &v, DoubleVector &b) const

void	internal_get_block_vectors (const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const

void	internal_get_block_vectors (const DoubleVector &v, Vector< DoubleVector > &s) const

void	internal_return_block_vectors (const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const

void	internal_return_block_vectors (const Vector< DoubleVector > &s, DoubleVector &v) const

void	internal_get_block (const unsigned &i, const unsigned &j, CRDoubleMatrix &output_block) const

void	internal_get_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix &output_block) const

int	internal_block_number (const unsigned &i_dof) const

int	internal_index_in_block (const unsigned &i_dof) const

const LinearAlgebraDistribution *	internal_block_distribution_pt (const unsigned &b) const
	Access function to the internal block distributions. More...

void	insert_auxiliary_block_distribution (const Vector< unsigned > &block_vec_number, LinearAlgebraDistribution *dist_pt)

void	block_matrix_test (const unsigned &i, const unsigned &j, const CRDoubleMatrix *block_matrix_pt) const

int	get_index_of_value (const Vector< myType > &vec, const myType val, const bool sorted=false) const

Public Member Functions inherited from oomph::Preconditioner
	Preconditioner ()
	Constructor. More...

	Preconditioner (const Preconditioner &)=delete
	Broken copy constructor. More...

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

virtual	~Preconditioner ()
	Destructor (empty) More...

virtual void	preconditioner_solve_transpose (const DoubleVector &r, DoubleVector &z)

void	setup (DoubleMatrixBase *matrix_pt)

void	setup (const Problem problem_pt, DoubleMatrixBase matrix_pt)

void	enable_silent_preconditioner_setup ()
	Set up the block preconditioner quietly! More...

void	disable_silent_preconditioner_setup ()
	Be verbose in the block preconditioner setup. More...

virtual void	set_matrix_pt (DoubleMatrixBase *matrix_pt)
	Set the matrix pointer. More...

virtual const OomphCommunicator *	comm_pt () const
	Get function for comm pointer. More...

virtual void	set_comm_pt (const OomphCommunicator *const comm_pt)
	Set the communicator pointer. More...

double	setup_time () const
	Returns the time to setup the preconditioner. More...

virtual void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)

virtual void	turn_into_subsidiary_block_preconditioner (BlockPreconditioner< CRDoubleMatrix > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse, const Vector< Vector< unsigned >> &doftype_coarsen_map_coarse)

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)

Public Attributes
Preconditioner *	S_00_preconditioner_pt
	Pointer to the S00 Schur Complement preconditioner. More...

MatrixBasedLumpedPreconditioner< CRDoubleMatrix > *	Lumped_J_11_preconditioner_pt
	Preconditioner for storing the lumped J_11 matrix. More...

MatrixBasedLumpedPreconditioner< CRDoubleMatrix > *	Lumped_J_22_preconditioner_pt
	Preconditioner for storing the lumped J_22 matrix. More...

DenseMatrix< CRDoubleMatrix * >	Matrix_of_block_pointers

CRDoubleMatrix *	S_00_pt

unsigned	Use_amg

Additional Inherited Members
Protected Member Functions inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
void	set_nmesh (const unsigned &n)

void	set_mesh (const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)

void	set_replacement_dof_block (const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix *replacement_dof_block_pt)

bool	any_mesh_distributed () const

int	internal_dof_number (const unsigned &i_dof) const

unsigned	internal_index_in_dof (const unsigned &i_dof) const

unsigned	internal_block_dimension (const unsigned &b) const

unsigned	internal_dof_block_dimension (const unsigned &i) const

unsigned	master_nrow () const

unsigned	internal_master_dof_number (const unsigned &b) const

const LinearAlgebraDistribution *	internal_preconditioner_matrix_distribution_pt () const

const LinearAlgebraDistribution *	preconditioner_matrix_distribution_pt () const

Protected Member Functions inherited from oomph::DistributableLinearAlgebraObject
void	clear_distribution ()

Protected Attributes inherited from oomph::BlockPreconditioner< CRDoubleMatrix >
MapMatrix< unsigned, CRDoubleMatrix * >	Replacement_dof_block_pt
	The replacement dof-level blocks. More...

Vector< LinearAlgebraDistribution * >	Block_distribution_pt
	The distribution for the blocks. More...

Vector< Vector< unsigned > >	Block_to_dof_map_coarse

Vector< Vector< unsigned > >	Block_to_dof_map_fine
	Mapping for the block types to the most fine grain dof types. More...

Vector< Vector< unsigned > >	Doftype_coarsen_map_coarse

Vector< Vector< unsigned > >	Doftype_coarsen_map_fine

Vector< LinearAlgebraDistribution * >	Internal_block_distribution_pt
	Storage for the default distribution for each internal block. More...

Vector< LinearAlgebraDistribution * >	Dof_block_distribution_pt

Vector< unsigned >	Allow_multiple_element_type_in_mesh

Vector< const Mesh * >	Mesh_pt

Vector< unsigned >	Ndof_types_in_mesh

unsigned	Internal_nblock_types

unsigned	Internal_ndof_types

Protected Attributes inherited from oomph::Preconditioner
bool	Silent_preconditioner_setup
	Boolean to indicate whether or not the build should be done silently. More...

std::ostream *	Stream_pt
	Pointer to the output stream – defaults to std::cout. More...

Detailed Description

SubBiharmonic Preconditioner - an inexact preconditioner for the 3x3 top left hand corner sub block matrix. Used as part of the BiharmonicPreconditioner<MATRIX>

Constructor & Destructor Documentation

◆ InexactSubBiharmonicPreconditioner() [1/2]

oomph::InexactSubBiharmonicPreconditioner::InexactSubBiharmonicPreconditioner	(	BiharmonicPreconditioner *	master_prec_pt,
		const bool	use_amg
	)

inline

Constructor for the inexact block preconditioner.
This a helper class for BiharmonicPreconditioner and cannot be used as a standalone preconditioner.
master_prec_pt is the pointer to the master BiharmonicPreconditioner.

     {
       // Block mapping for ExactSubBiharmonicPreconditioner
       Vector<unsigned> block_lookup(3);
       block_lookup[0] = 0;
       block_lookup[1] = 1;
       block_lookup[2] = 2;
  
       // set as subsidiary block preconditioner
       this->turn_into_subsidiary_block_preconditioner(master_prec_pt,
                                                       block_lookup);
  
       // zero all the preconditioner pt
       S_00_preconditioner_pt = 0;
       Lumped_J_11_preconditioner_pt = 0;
       Lumped_J_22_preconditioner_pt = 0;
  
       // store the preconditioner type
       Use_amg = use_amg;
     }

References Lumped_J_11_preconditioner_pt, Lumped_J_22_preconditioner_pt, S_00_preconditioner_pt, oomph::BlockPreconditioner< CRDoubleMatrix >::turn_into_subsidiary_block_preconditioner(), and Use_amg.

◆ ~InexactSubBiharmonicPreconditioner()

oomph::InexactSubBiharmonicPreconditioner::~InexactSubBiharmonicPreconditioner ( )

inline

destructor - just calls this->clean_up_memory()

     {
       this->clean_up_memory();
     }

References clean_up_memory().

◆ InexactSubBiharmonicPreconditioner() [2/2]

oomph::InexactSubBiharmonicPreconditioner::InexactSubBiharmonicPreconditioner ( const InexactSubBiharmonicPreconditioner & )

delete

Broken copy constructor.

Member Function Documentation

◆ clean_up_memory()

void oomph::InexactSubBiharmonicPreconditioner::clean_up_memory ( )

inlinevirtual

clean the memory

Reimplemented from oomph::Preconditioner.

     {
       // delete the sub preconditioner pt
       delete S_00_preconditioner_pt;
       S_00_preconditioner_pt = 0;
  
       // delete the lumped preconditioners
       delete Lumped_J_11_preconditioner_pt;
       Lumped_J_11_preconditioner_pt = 0;
       delete Lumped_J_22_preconditioner_pt;
       Lumped_J_22_preconditioner_pt = 0;
  
       // Number of block types
       unsigned n = Matrix_of_block_pointers.nrow();
  
       // delete the block matrices
       for (unsigned i = 0; i < n; i++)
       {
         for (unsigned j = 0; j < n; j++)
         {
           if (Matrix_of_block_pointers(i, j) != 0)
           {
             delete Matrix_of_block_pointers(i, j);
             Matrix_of_block_pointers(i, j) = 0;
           }
         }
       }
     }

References i, j, Lumped_J_11_preconditioner_pt, Lumped_J_22_preconditioner_pt, Matrix_of_block_pointers, n, and S_00_preconditioner_pt.

Referenced by setup(), and ~InexactSubBiharmonicPreconditioner().

◆ compute_inexact_schur_complement()

void oomph::InexactSubBiharmonicPreconditioner::compute_inexact_schur_complement ( )

Computes the inexact schur complement of the block J_00 using lumping as an approximate inverse of block J_11 and J_22 (where possible)

computes the schur complement for the inexact sub biharmonic preconditioner

   {
     // if required get pointers to the vector components of J01 and J10
     int* J_01_row_start = 0;
     int* J_01_column_index = 0;
     double* J_01_value = 0;
     int* J_10_row_start = 0;
     int* J_10_column_index = 0;
  
     // J_01 matrix
     J_01_row_start = Matrix_of_block_pointers(0, 1)->row_start();
     J_01_column_index = Matrix_of_block_pointers(0, 1)->column_index();
     J_01_value = Matrix_of_block_pointers(0, 1)->value();
  
     // J_10 matrix
     J_10_row_start = Matrix_of_block_pointers(1, 0)->row_start();
     J_10_column_index = Matrix_of_block_pointers(1, 0)->column_index();
  
     // if required get pointers to the vector components of J01 and J10
     int* J_02_row_start = 0;
     int* J_02_column_index = 0;
     double* J_02_value = 0;
     int* J_20_row_start = 0;
     int* J_20_column_index = 0;
  
     // J_02 matrix
     J_02_row_start = Matrix_of_block_pointers(0, 2)->row_start();
     J_02_column_index = Matrix_of_block_pointers(0, 2)->column_index();
     J_02_value = Matrix_of_block_pointers(0, 2)->value();
  
     // J_20 matrix
     J_20_row_start = Matrix_of_block_pointers(2, 0)->row_start();
     J_20_column_index = Matrix_of_block_pointers(2, 0)->column_index();
  
     // get the inverse lumped vector of J_11 if required
     double* J_11_lumped_and_inverted = 0;
     J_11_lumped_and_inverted =
       Lumped_J_11_preconditioner_pt->inverse_lumped_vector_pt();
  
     // get the inverse lumped vector of J_22 if required
     double* J_22_lumped_and_inverted = 0;
     J_22_lumped_and_inverted =
       Lumped_J_22_preconditioner_pt->inverse_lumped_vector_pt();
  
     // size of J00 matrix (and S00 matrix)
     unsigned J_00_nrow = Matrix_of_block_pointers(0, 0)->nrow();
  
     // vectors for the schur complement
     Vector<int> S_00_row_start(J_00_nrow + 1);
     Vector<int> S_00_column_index;
     Vector<double> S_00_value;
  
     // number of elements in the x-dimension of the mesh
     unsigned n_element_x =
       dynamic_cast<HermiteQuadMesh<Hijacked<BiharmonicElement<2>>>*>(
         dynamic_cast<BiharmonicPreconditioner*>(
           this->master_block_preconditioner_pt())
           ->bulk_element_mesh_pt())
         ->nelement_in_dim(0);
  
     // nnz in schur complement (initialised to zero)
     unsigned S_00_nnz = 0;
  
     // loop over columns of schur complement matrix
     for (unsigned i = 0; i < J_00_nrow; i++)
     {
       // set column_start
       S_00_row_start[i] = S_00_nnz;
  
       // loop over rows in schur complement matrix
       // the schur complement matrix has 5 well defined bands thus we only
       // perform matrix-matrix multiplication for these bands
       //
       // where the diagonal band is 0 :
       //
       //   band 1 :   -2*n_element_x +/- 5
       //        2 :   -n_element_x +/- 3
       //        3 :   0 +/- 3
       //        4 :   n_element_x +/- 3
       //        5 :   2*n_element_x +/- 5
       //
       // regardless of the type or combination of boundary conditions applied
  
       // Vector for postion of the bands in S_00
       Vector<std::pair<int, int>> band_position(5);
  
       // compute the minimum and maximum positions of each band in terms of
       // row number for column j
       // note : static_cast used because max and min don't work on unsigned
       band_position[0].first =
         std::max(0, static_cast<int>(i - n_element_x * 2 - 5));
       band_position[0].second =
         std::max(0,
                  std::min(static_cast<int>(J_00_nrow - 1),
                           static_cast<int>(i - n_element_x * 2 + 5)));
       band_position[1].first =
         std::max(band_position[0].second + 1,
                  std::max(0, static_cast<int>(i - n_element_x - 3)));
       band_position[1].second =
         std::max(0,
                  std::min(static_cast<int>(J_00_nrow - 1),
                           static_cast<int>(i - n_element_x + 3)));
       band_position[2].first = std::max(band_position[1].second + 1,
                                         std::max(0, static_cast<int>(i - 3)));
       band_position[2].second = std::max(
         0, std::min(static_cast<int>(J_00_nrow - 1), static_cast<int>(i + 3)));
       band_position[3].first =
         std::max(band_position[2].second + 1,
                  std::max(0, static_cast<int>(i + n_element_x - 3)));
       band_position[3].second =
         std::max(0,
                  std::min(static_cast<int>(J_00_nrow - 1),
                           static_cast<int>(i + n_element_x + 3)));
       band_position[4].first =
         std::max(band_position[3].second + 1,
                  std::max(0, static_cast<int>(i + n_element_x * 2 - 5)));
       band_position[4].second =
         std::max(0,
                  std::min(static_cast<int>(J_00_nrow - 1),
                           static_cast<int>(i + n_element_x * 2 + 5)));
  
       // number of bands
       unsigned n_band = 5;
  
       // loop over the bands
       for (unsigned b = 0; b < n_band; b++)
       {
         // loop over the rows in band b
         for (unsigned j = static_cast<unsigned>(band_position[b].first);
              j <= static_cast<unsigned>(band_position[b].second);
              j++)
         {
           ;
  
           // temporary value for the computation of S00(i,j)
           double temp_value = Matrix_of_block_pointers(0, 0)->operator()(i, j);
  
           // iterate through non-zero entries of  column j of A_10
           for (int k = J_01_row_start[i]; k < J_01_row_start[i + 1]; k++)
           {
             if (J_10_column_index[J_10_row_start[J_01_column_index[k]]] <=
                   static_cast<int>(j) &&
                 static_cast<int>(j) <=
                   J_10_column_index[J_10_row_start[J_01_column_index[k] + 1] -
                                     1])
             {
               temp_value -= J_01_value[k] *
                             Matrix_of_block_pointers(1, 0)->operator()(
                               J_01_column_index[k], j) *
                             J_11_lumped_and_inverted[J_01_column_index[k]];
             }
           }
  
           // next compute contribution for A_02*lumped(A_22)'*A_20
  
           // iterate through non-zero entries of  column j of A_10
           for (int k = J_02_row_start[i]; k < J_02_row_start[i + 1]; k++)
           {
             if (J_20_column_index[J_20_row_start[J_02_column_index[k]]] <=
                   static_cast<int>(j) &&
                 static_cast<int>(j) <=
                   J_20_column_index[J_20_row_start[J_02_column_index[k] + 1] -
                                     1])
             {
               temp_value -= J_02_value[k] *
                             Matrix_of_block_pointers(2, 0)->operator()(
                               J_02_column_index[k], j) *
                             J_22_lumped_and_inverted[J_02_column_index[k]];
             }
           }
  
           // add element to schur complement matrix S00
           if (temp_value != 0.0)
           {
             S_00_nnz++;
             S_00_value.push_back(temp_value);
             S_00_column_index.push_back(j);
           }
         }
       }
     }
  
     // last entry of s00 column start
     S_00_row_start[J_00_nrow] = S_00_nnz;
  
     // build the schur complement S00
     S_00_pt = new CRDoubleMatrix(this->block_distribution_pt(0),
                                  J_00_nrow,
                                  S_00_value,
                                  S_00_column_index,
                                  S_00_row_start);
  
     // replace block J01 with J01*lumped(J11)' (if J11 can be lumped)
     unsigned J_01_nnz = Matrix_of_block_pointers(0, 1)->nnz();
     for (unsigned i = 0; i < J_01_nnz; i++)
     {
       J_01_value[i] *= J_11_lumped_and_inverted[J_01_column_index[i]];
     }
  
     // replace block J_02 with J_02*lumped(J_22)' (if J22 can be lumped)
     unsigned J_02_nnz = Matrix_of_block_pointers(0, 2)->nnz();
     for (unsigned i = 0; i < J_02_nnz; i++)
     {
       J_02_value[i] *= J_22_lumped_and_inverted[J_02_column_index[i]];
     }
   }

References b, oomph::BlockPreconditioner< CRDoubleMatrix >::block_distribution_pt(), i, j, k, Lumped_J_11_preconditioner_pt, Lumped_J_22_preconditioner_pt, oomph::BlockPreconditioner< CRDoubleMatrix >::master_block_preconditioner_pt(), Matrix_of_block_pointers, max, min, and S_00_pt.

Referenced by setup().

◆ operator=()

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

delete

Broken assignment operator.

◆ preconditioner_solve()

void oomph::InexactSubBiharmonicPreconditioner::preconditioner_solve	(	const DoubleVector &	r,
		DoubleVector &	z
	)

virtual

Apply preconditioner to r.

preconditioner solve for the inexact sub biharmonic preconditioner

Implements oomph::Preconditioner.

   {
     // get the block vectors
     Vector<DoubleVector> block_r(3);
     get_block_vectors(r, block_r);
  
     // r_0 = r_0 - J_01 * lumped(J_11)'*r_1 - J_02 * lumped(J_22)'*r_2
     // Remember that J_01 has already been premultiplied by lumped(J_11)
     DoubleVector temp;
     Matrix_of_block_pointers(0, 1)->multiply(block_r[1], temp);
     block_r[0] -= temp;
     temp.clear();
     Matrix_of_block_pointers(0, 2)->multiply(block_r[2], temp);
     block_r[0] -= temp;
  
     // apply the inexact preconditioner
     temp.clear();
     S_00_preconditioner_pt->preconditioner_solve(block_r[0], temp);
     return_block_vector(0, temp, z);
  
     // solve: lumped(J_11) x_1 = r_1 - J_10 x_0 for x_1
     // remember temp contains r_0 (...or z_0)
     DoubleVector temp2;
     Matrix_of_block_pointers(1, 0)->multiply(temp, temp2);
     block_r[1] -= temp2;
     DoubleVector z_1;
     Lumped_J_11_preconditioner_pt->preconditioner_solve(block_r[1], z_1);
     return_block_vector(1, z_1, z);
  
     // solve: lumped(J_22) x_2 = r_2 - J_20 x_0 for x_2
     // remember temp contains r_0 (...or z_0)
     temp2.clear();
     Matrix_of_block_pointers(2, 0)->multiply(temp, temp2);
     block_r[2] -= temp2;
     DoubleVector z_2;
     Lumped_J_22_preconditioner_pt->preconditioner_solve(block_r[2], z_2);
     return_block_vector(2, z_2, z);
   }

References oomph::DoubleVector::clear(), oomph::BlockPreconditioner< CRDoubleMatrix >::get_block_vectors(), Lumped_J_11_preconditioner_pt, Lumped_J_22_preconditioner_pt, Matrix_of_block_pointers, oomph::Preconditioner::preconditioner_solve(), UniformPSDSelfTest::r, oomph::BlockPreconditioner< CRDoubleMatrix >::return_block_vector(), and S_00_preconditioner_pt.

◆ setup()

void oomph::InexactSubBiharmonicPreconditioner::setup ( )

virtual

Setup the preconditioner.

setup for the inexact sub biharmonic preconditioner

Implements oomph::Preconditioner.

   {
     // clean up memory first
     this->clean_up_memory();
  
     // setup
     this->block_setup();
  
     // Number of block types
     unsigned n_block_types = this->nblock_types();
  
     // paranoid check for number of blocks
 #ifdef PARANOID
     if (n_block_types != 3)
     {
       std::ostringstream error_message;
       error_message
         << "This preconditioner requires 3 block types.\n"
         << "It is sub preconditioner for the BiharmonicPreconditioner.\n";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
     }
 #endif
  
     // required blocks
     DenseMatrix<bool> required_blocks(n_block_types, n_block_types, false);
     required_blocks(0, 0) = true;
     required_blocks(0, 1) = true;
     required_blocks(1, 0) = true;
     required_blocks(1, 1) = true;
     required_blocks(0, 2) = true;
     required_blocks(2, 0) = true;
     required_blocks(2, 2) = true;
  
     // Matrix of block matrix pointers
     Matrix_of_block_pointers.resize(n_block_types, n_block_types, 0);
  
     // get the blocks
     this->get_blocks(required_blocks, Matrix_of_block_pointers);
  
     // lump the matrix J_11
     Lumped_J_11_preconditioner_pt =
       new MatrixBasedLumpedPreconditioner<CRDoubleMatrix>;
     Lumped_J_11_preconditioner_pt->setup(Matrix_of_block_pointers(1, 1));
  
     delete Matrix_of_block_pointers(1, 1);
     Matrix_of_block_pointers(1, 1) = 0;
  
     // lump the matrix J_22
     Lumped_J_22_preconditioner_pt =
       new MatrixBasedLumpedPreconditioner<CRDoubleMatrix>;
     Lumped_J_22_preconditioner_pt->setup(Matrix_of_block_pointers(2, 2));
     delete Matrix_of_block_pointers(2, 2);
     Matrix_of_block_pointers(2, 2) = 0;
  
     // compute the schur complement
     compute_inexact_schur_complement();
  
     // create the preconditioner for the S00 Schur complement linear system
     if (Use_amg)
     {
 #ifdef OOMPH_HAS_HYPRE
       // Use Hypre Boomer AMG
       S_00_preconditioner_pt = new HyprePreconditioner;
       Biharmonic_schur_complement_Hypre_defaults::set_defaults(
         static_cast<HyprePreconditioner*>(S_00_preconditioner_pt));
 #else
       std::ostringstream error_message;
       error_message << "Request AMG solver but oomph-lib does not have HYPRE";
       throw OomphLibError(
         error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
 #endif
     }
     else
     {
       S_00_preconditioner_pt = new SuperLUPreconditioner;
     }
  
     // setup the preconditioner
     S_00_preconditioner_pt->setup(S_00_pt);
  
     // clean up
     delete S_00_pt;
     S_00_pt = 0;
   }