block_preconditioner.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Include guards
#ifndef OOMPH_BLOCK_PRECONDITION_HEADER
#define OOMPH_BLOCK_PRECONDITION_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// c++ include
#include <math.h>
#include <typeinfo>
 
// oomph-lib includes
#include "matrices.h"
#include "mesh.h"
#include "vector_matrix.h"
 
// #include "problem.h"
#include "preconditioner.h"
#include "SuperLU_preconditioner.h"
#include "matrix_vector_product.h"
 
namespace oomph
{
  //=============================================================================
  //=============================================================================
  class BlockSelector
  {
  public:
    BlockSelector()
    {
      // Needs to be set to zero because if the build function leaves the
      // Replacement_block_pt alone if replacement_block_pt = 0 (the default
      // argument).
      Replacement_block_pt = 0;
      this->build(0, 0, false);
    }
 
    BlockSelector(const unsigned& row_index,
                  const unsigned& column_index,
                  const bool& wanted,
                  CRDoubleMatrix* replacement_block_pt = 0)
    {
#ifdef PARANOID
      if ((wanted == false) && (replacement_block_pt != 0))
      {
        std::ostringstream err_msg;
        err_msg << "Trying to construct a BlockSelector object with:\n"
                << "replacement_block_pt != 0 and wanted == false"
                << "If you require the block, please set wanted == true.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Needs to be set to zero because if the build function leaves the
      // Replacement_block_pt alone if replacement_block_pt = 0 (the default
      // argument). Thus if it is not set here, it would not be initialised to
      // null.
      Replacement_block_pt = 0;
 
      this->build(row_index, column_index, wanted, replacement_block_pt);
    }
 
    virtual ~BlockSelector()
    {
#ifdef PARANOID
      if (Replacement_block_pt != 0)
      {
        std::ostringstream warning_msg;
        warning_msg << "Warning: BlockSelector destructor is called but...\n"
                    << "replacement_block_pt() is not null.\n"
                    << "Please remember to null this via the function\n"
                    << "BlockSelector::null_replacement_block_pt()\n"
                    << "Row_index: " << Row_index << "\n"
                    << "Column_index: " << Column_index << std::endl;
 
        OomphLibWarning(
          warning_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
 
    void select_block(const unsigned& row_index,
                      const unsigned& column_index,
                      const bool& wanted,
                      CRDoubleMatrix* replacement_block_pt = 0)
    {
#ifdef PARANOID
      if ((wanted == false) && (replacement_block_pt != 0))
      {
        std::ostringstream err_msg;
        err_msg << "Trying to select block with:\n"
                << "replacement_block_pt != 0 and wanted == false"
                << "If you require the block, please set wanted == true.\n"
                << "row_index: " << row_index << "\n"
                << "column_index: " << column_index << "\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      this->build(row_index, column_index, wanted, replacement_block_pt);
    }
 
 
    void want_block()
    {
      Wanted = true;
    }
 
    void do_not_want_block()
    {
#ifdef PARANOID
      if (Replacement_block_pt != 0)
      {
        std::ostringstream warning_msg;
        warning_msg << "Trying to set Wanted = false, but replacement_block_pt "
                       "is not null.\n"
                    << "Please call null_replacement_block_pt()\n"
                    << "(remember to free memory if necessary)\n"
                    << "Row_index: " << Row_index << "\n"
                    << "Column_index: " << Column_index << "\n";
        OomphLibWarning(
          warning_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      null_replacement_block_pt();
 
      Wanted = false;
    }
 
    void null_replacement_block_pt()
    {
      Replacement_block_pt = 0;
    }
 
    void set_replacement_block_pt(CRDoubleMatrix* replacement_block_pt)
    {
#ifdef PARANOID
      if (Wanted == false)
      {
        std::ostringstream err_msg;
        err_msg << "Trying to set replacement_block_pt, but Wanted == false.\n"
                << "Please call want_block()\n"
                << "Row_index: " << Row_index << "\n"
                << "Column_index: " << Column_index << "\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      Replacement_block_pt = replacement_block_pt;
    }
 
    CRDoubleMatrix* replacement_block_pt() const
    {
      return Replacement_block_pt;
    }
 
    void set_row_index(const unsigned& row_index)
    {
      Row_index = row_index;
    }
 
    const unsigned& row_index() const
    {
      return Row_index;
    }
 
    void set_column_index(const unsigned& column_index)
    {
      Column_index = column_index;
    }
 
    const unsigned& column_index() const
    {
      return Column_index;
    }
 
    const bool& wanted() const
    {
      return Wanted;
    }
 
 
    friend std::ostream& operator<<(std::ostream& o_stream,
                                    const BlockSelector& block_selector)
    {
      o_stream << "Row_index = " << block_selector.row_index() << "\n"
               << "Column_index = " << block_selector.column_index() << "\n"
               << "Wanted = " << block_selector.wanted() << "\n"
               << "Replacement_block_pt = ";
      if (block_selector.replacement_block_pt() == 0)
      {
        o_stream << 0;
      }
 
      return o_stream;
    }
 
  private:
    void build(const unsigned& row_index,
               const unsigned& column_index,
               const bool& wanted,
               CRDoubleMatrix* replacement_block_pt = 0)
    {
      Row_index = row_index;
      Column_index = column_index;
      Wanted = wanted;
 
      // Only set the replacement_block_pt if it is wanted. Otherwise we leave
      // it alone. All constructors should set Replacement_block_pt to 0.
      if (replacement_block_pt != 0)
      {
#ifdef PARANOID
        if (Wanted == false)
        {
          std::ostringstream err_msg;
          err_msg
            << "Trying to set replacement_block_pt, but Wanted == false.\n"
            << "Please either not set the replacement_block_pt or call the "
               "function\n"
            << "do_not_want_block()\n"
            << "Row_index: " << Row_index << "\n"
            << "Column_index: " << Column_index << "\n";
          throw OomphLibError(
            err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
#endif
 
        Replacement_block_pt = replacement_block_pt;
      }
    }
 
    unsigned Row_index;
 
    unsigned Column_index;
 
    bool Wanted;
 
    CRDoubleMatrix* Replacement_block_pt;
  };
 
 
  //============================================================================
  //============================================================================
  template<typename MATRIX>
  class BlockPreconditioner : public Preconditioner
  {
  public:
    BlockPreconditioner() : Ndof_types_in_mesh(0)
    {
      // Initially set the master block preconditioner pointer to zero
      // indicating that this is stand-alone preconditioner (i.e. the upper most
      // level preconditioner) that will set up its own block lookup schemes
      // etc.
      Master_block_preconditioner_pt = 0;
 
      // The distribution of the concatenation of the internal block
      // distributions.
      // I.e. LinearAlgebraDistributionHelpers::concatenate
      //        (distributions of internal blocks).
      //
      // The concatenation of the internal block distributions is stored in two
      // places depending on if this is the upper-most master block
      // preconditioner or not.
      //
      // If this is a master block preconditioner
      // (Master_block_preconditioner_pt is null), then it is stored in the
      // variable Internal_preconditioner_matrix_distribution_pt (below). For
      // subsidiary block preconditioners, this remains null.
      //
      // Because BlockPreconditioners are DistributedLinearAlgebraObjects, they
      // have a distribution. For the upper-most master block preconditioner,
      // this is the distribution of the underlying Jacobian.
      //
      // For all subsidiary block preconditioners, this remains null. The
      // concatenation of the distribution of the internal blocks are stored
      // as the distribution of the BlockPreconditioner.
      //
      // This seems inconsistent and I cannot figure out why this is done.
      Internal_preconditioner_matrix_distribution_pt = 0;
 
      // The concatenation of the external block distributions.
      Preconditioner_matrix_distribution_pt = 0;
 
      // Initialise number of rows in this block preconditioner.
      // This information is maintained if used as subsidiary or
      // stand-alone block preconditioner (in the latter case it
      // obviously stores the number of rows within the subsidiary
      // preconditioner.
      Nrow = 0;
 
      // Initialise number of different block types in this preconditioner.
      // This information is maintained if used as subsidiary or
      // stand-alone block preconditioner (in the latter case it
      // obviously stores the number of rows within the subsidiary
      // preconditioner.
      Internal_nblock_types = 0;
 
      // Initialise number of different dof types in this preconditioner.
      // This information is maintained if used as subsidiary or
      // stand-alone block preconditioner (in the latter case it
      // obviously stores the number of rows within the subsidiary
      // preconditioner.
      Internal_ndof_types = 0;
 
      // There are no blocks to start off with.
      Block_distribution_pt.resize(0);
 
      // The distributions of the underlying internal blocks.
      Internal_block_distribution_pt.resize(0);
 
      // The distribution of the dof-level blocks, these are used during the
      // concatenation process to create the distribution of the blocks.
      Dof_block_distribution_pt.resize(0);
 
      // Clear both the Block_to_dof_map_coarse and Block_to_dof_map_fine
      // vectors.
      Block_to_dof_map_coarse.resize(0);
      Block_to_dof_map_fine.resize(0);
 
      // Default the debug flag to false.
      Recursive_debug_flag = false;
 
      // Default the debug flag to false.
      Debug_flag = false;
    } // EOFunc constructor
 
 
    virtual ~BlockPreconditioner()
    {
      this->clear_block_preconditioner_base();
    } // EOFunc destructor
 
    BlockPreconditioner(const BlockPreconditioner&) = delete;
 
    void operator=(const BlockPreconditioner&) = delete;
 
    MATRIX* matrix_pt() const
    {
      if (is_subsidiary_block_preconditioner())
      {
        return master_block_preconditioner_pt()->matrix_pt();
      }
      else
      {
        MATRIX* m_pt = dynamic_cast<MATRIX*>(Preconditioner::matrix_pt());
#ifdef PARANOID
        if (m_pt == 0)
        {
          std::ostringstream error_msg;
          error_msg << "Matrix is not correct type.";
          throw OomphLibError(
            error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
#endif
        return m_pt;
      }
    } // EOFunc matrix_pt()
 
 
    void turn_on_recursive_debug_flag()
    {
      Recursive_debug_flag = true;
      if (is_subsidiary_block_preconditioner())
        this->master_block_preconditioner_pt()->turn_on_recursive_debug_flag();
    }
 
    void turn_off_recursive_debug_flag()
    {
      Recursive_debug_flag = false;
      if (is_subsidiary_block_preconditioner())
        this->master_block_preconditioner_pt()->turn_off_recursive_debug_flag();
    }
 
    void turn_on_debug_flag()
    {
      Debug_flag = true;
    }
 
    void turn_off_debug_flag()
    {
      Debug_flag = false;
    }
 
    void turn_into_subsidiary_block_preconditioner(
      BlockPreconditioner<MATRIX>* master_block_prec_pt,
      const Vector<unsigned>& doftype_in_master_preconditioner_coarse);
 
    void turn_into_subsidiary_block_preconditioner(
      BlockPreconditioner<MATRIX>* 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,
                   MATRIX& output_matrix,
                   const bool& ignore_replacement_block = false) const
    {
#ifdef PARANOID
      // Check the range of i and j, they should not be greater than
      // nblock_types
      unsigned n_block_types = this->nblock_types();
      if ((i > n_block_types) || (j > n_block_types))
      {
        std::ostringstream err_msg;
        err_msg << "Requested block(" << i << "," << j << ")"
                << "\n"
                << "but nblock_types() is " << n_block_types << ".\n"
                << "Maybe your argument to block_setup(...) is incorrect?"
                << std::endl;
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // The logic is this:
      //
      // Block_to_dof_map_coarse tells us which dof types belongs in each block.
      // This is relative to this preconditioner, and describes the external
      // block and dof type mappings (what the preconditioner writer
      // expects/sees).
      //
      // As such, the dof types in Block_to_dof_map_coarse describes the
      // dof-level blocks needed to be concatenated to produce block(i,j). These
      // dofs may have been coarsened.
      //
      // Now, the dof blocks to concatenate comes from:
      // If the dof block exists in Replacement_dof_block_pt, then we make a
      // deep copy of it. Otherwise, if this is the upper-most master block
      // preconditioner then we get it from the original matrix via function
      // internal_get_block(...) otherwise, if this is a subsidiary
      // block preconditioner, we go one level up the hierarchy and repeat the
      // process.
      //
      //
      // A small note about indirections which caused me some headache.
      // Thought I would mention it here in case the below code needs to be
      // re-visited.
      //
      // A small subtlety with indirections:
      //
      // The underlying ordering of the dof-level blocks is STILL AND ALWAYS the
      // `natural' ordering determined by first the elements then the order of
      // the meshes.
      //
      // But during the process of block_setup(...), the external (perceived)
      // block ordering may have changed. So some indirection has to take place,
      // this mapping is maintained in Block_to_dof_map_coarse.
      //
      // For example, take the Lagrangian preconditioner, which expects the
      // natural dof type ordering:
      // 0 1 2 3  4  5
      // u v p uc vc L
      //
      // If the mapping given to block_setup(...) is:
      // dof_to_block_map = [0, 1, 4, 2, 3, 5]
      // saying that: dof type 0 goes to block 0
      //              dof type 1 goes to block 1
      //              dof type 2 goes to block 4
      //              dof type 3 goes to block 2
      //              dof type 4 goes to block 3
      //              dof type 5 goes to block 5
      //
      // The function block_setup(...) will populate the vector
      // Block_to_dof_map_coarse with [[0][1][3][4][2][5]],
      // which says that get_block(0,0) will get the u block
      //                 get_block(1,1) will get the v block
      //                 get_block(2,2) will get the uc block
      //                 get_block(3,3) will get the vc block
      //                 get_block(4,4) will get the p block
      //                 get_block(5,5) will get the L block
      //
      // i.e. the ordering of the block types is a permutation of the dof types,
      // so that we now have the following block ordering:
      // 0 1 2  3  4 5 <- block ordering
      // u v uc vc p L
      //
      // Now, the writer expects to work with the block ordering. I.e. when we
      // replace a dof-level block, say the pressure block, we call
      // set_replacement_dof_block(4,4,Matrix);
      // Similarly, when we want the pressure block, we call
      // get_block(4,4);
      //
      // Now, the below code uses the indirection maintained in
      // Block_to_dof_map_coarse. I.e. When we call get_block(4,4), we use the
      // mapping Block_to_dof_map_coarse[4] = 2, we get the block (2,2)
      // (from Replacement_dof_block_pt or internal_get_block), since the
      // underlying dof_to_block mapping is still the identity, i.e. it has not
      // changed from:
      // 0 1 2 3  4  5
      // u v p uc vc L
      //
      // So, block (4,4) (pressure block) maps to the block (2,2).
 
      // How many external dof types are in this block?
      const unsigned ndofblock_in_row = Block_to_dof_map_coarse[i].size();
      const unsigned ndofblock_in_col = Block_to_dof_map_coarse[j].size();
 
      // If there is only one dof block in this block then there is no need to
      // concatenate.
      if (ndofblock_in_row == 1 && ndofblock_in_col == 1)
      {
        // Get the indirection for the dof we want.
        const unsigned wanted_dof_row = Block_to_dof_map_coarse[i][0];
        const unsigned wanted_dof_col = Block_to_dof_map_coarse[j][0];
 
        // If the block has NOT been replaced or if we want to ignore the
        // replacement, then we call get_dof_level_block(...) which will get the
        // dof-level blocks up the preconditioning hierarchy, dragging down
        // any replacement dof blocks of parent preconditioners if required.
        if ((Replacement_dof_block_pt.get(wanted_dof_row, wanted_dof_col) ==
             0) ||
            ignore_replacement_block)
        {
          get_dof_level_block(wanted_dof_row,
                              wanted_dof_col,
                              output_matrix,
                              ignore_replacement_block);
        }
        else
        // Replacement_dof_block_pt.get(wanted_dof_row,wanted_dof_col) is not
        // null, this means that the block has been replaced. We simply make
        // a deep copy of it.
        {
          CRDoubleMatrixHelpers::deep_copy(
            Replacement_dof_block_pt.get(wanted_dof_row, wanted_dof_col),
            output_matrix);
        }
      }
      else
      // This block contains more than one dof-level block. So we need to
      // concatenate the (external) dof-level blocks.
      {
        // The CRDoubleMatrixHelpers::concatenate_without_communication(...)
        // takes a DenseMatrix of pointers to CRDoubleMatrices to concatenate.
        DenseMatrix<CRDoubleMatrix*> tmp_blocks_pt(
          ndofblock_in_row, ndofblock_in_col, 0);
 
        // Vector of Vectors of unsigns to indicate whether we have created
        // CRDoubleMatrices with new or not... so we can delete it later.
        // 0 - no new CRDoubleMatrix is created.
        // 1 - a new CRDoubleMatrix is created.
        // If we ever use C++11, remove this and use smart pointers.
        Vector<Vector<unsigned>> new_block(
          ndofblock_in_row, Vector<unsigned>(ndofblock_in_col, 0));
 
        // Loop through the number of dof block rows and then the number of dof
        // block columns, either get the pointer from Replacement_dof_block_pt
        // or from get_dof_level_block(...).
        for (unsigned row_dofblock = 0; row_dofblock < ndofblock_in_row;
             row_dofblock++)
        {
          // Indirection for the row, as discuss in the large chunk of text
          // previously.
          const unsigned wanted_dof_row =
            Block_to_dof_map_coarse[i][row_dofblock];
 
          for (unsigned col_dofblock = 0; col_dofblock < ndofblock_in_col;
               col_dofblock++)
          {
            // Indirection for the column as discussed in the large chunk of
            // text above.
            const unsigned wanted_dof_col =
              Block_to_dof_map_coarse[j][col_dofblock];
 
            // Get the pointer from Replacement_dof_block_pt.
            tmp_blocks_pt(row_dofblock, col_dofblock) =
              Replacement_dof_block_pt.get(wanted_dof_row, wanted_dof_col);
 
            // If the pointer from Replacement_dof_block_pt is null, or if
            // we have to ignore replacement blocks, build a new CRDoubleMatrix
            // via get_dof_level_block.
            if ((tmp_blocks_pt(row_dofblock, col_dofblock) == 0) ||
                ignore_replacement_block)
            {
              // We have to create a new CRDoubleMatrix, so put in 1 to indicate
              // that we have to delete it later.
              new_block[row_dofblock][col_dofblock] = 1;
 
              // Create the new CRDoubleMatrix. Note that we do not use the
              // indirection, since the indirection is only used one way.
              tmp_blocks_pt(row_dofblock, col_dofblock) = new CRDoubleMatrix;
 
              // Get the dof-level block, as discussed above.
              get_dof_level_block(wanted_dof_row,
                                  wanted_dof_col,
                                  *tmp_blocks_pt(row_dofblock, col_dofblock),
                                  ignore_replacement_block);
            }
          }
        }
 
        // Concatenation of matrices require the distribution of the individual
        // sub-matrices (for both row and column). This is because concatenation
        // is implemented without communication in such a way that the rows
        // and column values are both permuted, the permutation is defined by
        // the individual distributions of the sub blocks.
        // Without a vector of distributions describing the distributions of
        // of the rows, we do not know how to permute the rows. For the columns,
        // although CRDoubleMatrices does not have a column distribution, the
        // concatenated matrix must have it's columns permuted, so we mirror
        // the diagonal and get the corresponding row distribution.
        //
        // Confused? - Example: Say we have a matrix with dof blocking
        //
        //   | a | b | c | d | e |
        // --|---|---|---|---|---|
        // a |   |   |   |   |   |
        // --|---|---|---|---|---|
        // b |   |   |   |   |   |
        // --|---|---|---|---|---|
        // c |   |   |   |   |   |
        // --|---|---|---|---|---|
        // d |   |   |   |   |   |
        // --|---|---|---|---|---|
        // e |   |   |   |   |   |
        // --|---|---|---|---|---|
        //
        // We wish to concatenate the blocks
        //
        //   | d | e |
        // --|---|---|
        // a |   |   |
        // --|---|---|
        // b |   |   |
        // --|---|---|
        //
        // Then clearly the row distributions required are the distributions for
        // the dof blocks a and b. But block(a,d) does not have a column
        // distribution since it is a CRDoubleMatrix! - We use the distribution
        // mirrored by the diagonal, so the column distributions required to
        // concatenate these blocks is the same as the distributions of the rows
        // for dof block d and e.
 
        // First we do the row distribution.
 
        // Storage for the row distribution pointers.
        Vector<LinearAlgebraDistribution*> tmp_row_dist_pt(ndofblock_in_row, 0);
 
        // Loop through the number of dof blocks in the row. For the upper-most
        // master block preconditioner, the external dof distributions is the
        // same as the internal BLOCK distributions. Recall what we said above
        // about the underlying blocks maintaining it's natural ordering.
        //
        // If this is a subsidiary block preconditioner, then the distributions
        // for the dof blocks are stored in Dof_block_distribution_pt. The
        // reason why this is different for subsidiary block preconditioners is
        // because subsidiary block preconditioners would have it's dof types
        // coarsened. Then a single dof distribution in a subsidiary block
        // preconditioner could be a concatenation of many dof distributions of
        // the master dof distributions.
        for (unsigned row_dof_i = 0; row_dof_i < ndofblock_in_row; row_dof_i++)
        {
          const unsigned mapped_dof_i = Block_to_dof_map_coarse[i][row_dof_i];
          if (is_master_block_preconditioner())
          {
            tmp_row_dist_pt[row_dof_i] =
              Internal_block_distribution_pt[mapped_dof_i];
          }
          else
          {
            tmp_row_dist_pt[row_dof_i] =
              Dof_block_distribution_pt[mapped_dof_i];
          }
        }
 
        // Storage for the column distribution pointers.
        Vector<LinearAlgebraDistribution*> tmp_col_dist_pt(ndofblock_in_col, 0);
 
        // We do the same thing as before.
        for (unsigned col_dof_i = 0; col_dof_i < ndofblock_in_col; col_dof_i++)
        {
          const unsigned mapped_dof_j = Block_to_dof_map_coarse[j][col_dof_i];
          if (is_master_block_preconditioner())
          {
            tmp_col_dist_pt[col_dof_i] =
              Internal_block_distribution_pt[mapped_dof_j];
          }
          else
          {
            tmp_col_dist_pt[col_dof_i] =
              Dof_block_distribution_pt[mapped_dof_j];
          }
        }
 
        // Perform the concatenation.
        CRDoubleMatrixHelpers::concatenate_without_communication(
          tmp_row_dist_pt, tmp_col_dist_pt, tmp_blocks_pt, output_matrix);
 
        // Delete any new CRDoubleMatrices we have created.
        for (unsigned row_i = 0; row_i < ndofblock_in_row; row_i++)
        {
          for (unsigned col_i = 0; col_i < ndofblock_in_col; col_i++)
          {
            if (new_block[row_i][col_i])
            {
              delete tmp_blocks_pt(row_i, col_i);
            }
          }
        }
      } // else need to concatenate
    } // EOFunc get_block(...)
 
 
    MATRIX get_block(const unsigned& i,
                     const unsigned& j,
                     const bool& ignore_replacement_block = false) const
    {
      MATRIX output_matrix;
      get_block(i, j, output_matrix, ignore_replacement_block);
      return output_matrix;
    } // EOFunc get_block(...)
 
    void set_master_matrix_pt(MATRIX* in_matrix_pt)
    {
      if (is_subsidiary_block_preconditioner())
      {
        master_block_preconditioner_pt()->set_master_matrix_pt(in_matrix_pt);
      }
      else
      {
        this->set_matrix_pt(in_matrix_pt);
      }
    }
 
    void get_block_other_matrix(const unsigned& i,
                                const unsigned& j,
                                MATRIX* in_matrix_pt,
                                MATRIX& output_matrix)
    {
      MATRIX* backup_matrix_pt = matrix_pt();
      set_master_matrix_pt(in_matrix_pt);
      get_block(i, j, output_matrix, true);
      set_master_matrix_pt(backup_matrix_pt);
    } // EOFunc get_block_other_matrix(...)
 
    void get_blocks(DenseMatrix<bool>& required_blocks,
                    DenseMatrix<MATRIX*>& block_matrix_pt) const;
 
    void get_dof_level_block(
      const unsigned& i,
      const unsigned& j,
      MATRIX& output_block,
      const bool& ignore_replacement_block = false) const;
 
 
    MATRIX get_concatenated_block(
      const VectorMatrix<BlockSelector>& selected_block)
    {
#ifdef PARANOID
 
      unsigned const para_selected_block_nrow = selected_block.nrow();
      unsigned const para_selected_block_ncol = selected_block.ncol();
      unsigned const para_nblocks = this->nblock_types();
 
      // Check that selected_block size is not 0.
      if (para_selected_block_nrow == 0)
      {
        std::ostringstream error_msg;
        error_msg << "selected_block has nrow = 0.\n";
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Check that the number of blocks is not outside of the range
      // nblock_types(). Since this function builds matrices for block
      // preconditioning, it does not make sense for us to concatenate more
      // blocks than nblock_types().
      if ((para_selected_block_nrow > para_nblocks) ||
          (para_selected_block_ncol > para_nblocks))
      {
        std::ostringstream error_msg;
        error_msg << "Trying to concatenate a " << para_selected_block_nrow
                  << " by " << para_selected_block_ncol << " block matrix,\n"
                  << "but there are only " << para_nblocks << " block types.\n";
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Check that selected_block make sense, i.e. the row indices of each row
      // are the same, and the column indices of each column are the same.
 
      // First check if the row indices are consistent.
      // For each row, loop through the columns, comparing the row index against
      // the first column.
      for (unsigned row_i = 0; row_i < para_selected_block_nrow; row_i++)
      {
        const unsigned col_0_row_index = selected_block[row_i][0].row_index();
 
        for (unsigned col_i = 0; col_i < para_selected_block_ncol; col_i++)
        {
          const unsigned para_b_i = selected_block[row_i][col_i].row_index();
          const unsigned para_b_j = selected_block[row_i][col_i].column_index();
 
          if (col_0_row_index != para_b_i)
          {
            std::ostringstream err_msg;
            err_msg << "Block index for block(" << row_i << "," << col_i << ") "
                    << "contains block indices (" << para_b_i << "," << para_b_j
                    << ").\n"
                    << "But the row index for the first column is "
                    << col_0_row_index << ", they must be the same!\n";
            throw OomphLibError(
              err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Do the same for the column indices, consistency check.
      // For each column, loop through the rows, comparing the column index
      // against the first row.
      for (unsigned col_i = 0; col_i < para_selected_block_ncol; col_i++)
      {
        const unsigned row_0_col_index =
          selected_block[0][col_i].column_index();
 
        for (unsigned row_i = 0; row_i < para_selected_block_nrow; row_i++)
        {
          const unsigned para_b_i = selected_block[row_i][col_i].row_index();
          const unsigned para_b_j = selected_block[row_i][col_i].column_index();
 
          if (row_0_col_index != para_b_j)
          {
            std::ostringstream err_msg;
            err_msg << "Block index for block(" << row_i << "," << col_i << ") "
                    << "contains block indices (" << para_b_i << "," << para_b_j
                    << ").\n"
                    << "But the col index for the first row is "
                    << row_0_col_index << ", they must be the same!\n";
            throw OomphLibError(
              err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Check to see if the values in selected_block is within the range
      // nblock_types()
      //
      // Since we know that the column and row indices are consistent (by the
      // two paranoia checks above), we only need to check the column indices
      // in the first row, and the row indices in the first column.
 
      // Check that the row indices in the first column are within the range
      // nblock_types()
      for (unsigned row_i = 0; row_i < para_selected_block_nrow; row_i++)
      {
        const unsigned para_b_i = selected_block[row_i][0].row_index();
        const unsigned para_b_j = selected_block[row_i][0].column_index();
        if (para_b_i > para_nblocks)
        {
          std::ostringstream err_msg;
          err_msg << "Block index for block(" << row_i << ",0) "
                  << "contains block indices (" << para_b_i << "," << para_b_j
                  << ").\n"
                  << "But there are only " << para_nblocks
                  << " nblock_types().\n";
          throw OomphLibError(
            err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Check that the col indices in the first row are within the range
      // nblock_types()
      for (unsigned col_i = 0; col_i < para_selected_block_ncol; col_i++)
      {
        const unsigned para_b_i = selected_block[0][col_i].row_index();
        const unsigned para_b_j = selected_block[0][col_i].column_index();
        if (para_b_j > para_nblocks)
        {
          std::ostringstream err_msg;
          err_msg << "Block index for block(0," << col_i << ") "
                  << "contains block indices (" << para_b_i << "," << para_b_j
                  << ").\n"
                  << "But there are only " << para_nblocks
                  << " nblock_types().\n";
          throw OomphLibError(
            err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Stricter test - can be removed is required in the future. For the first
      // column, check that the row indices does not repeat.
      std::set<unsigned> row_index_set;
      for (unsigned row_i = 0; row_i < para_selected_block_nrow; row_i++)
      {
        std::pair<std::set<unsigned>::iterator, bool> row_index_set_ret;
 
        const unsigned row_i_index = selected_block[row_i][0].row_index();
 
        row_index_set_ret = row_index_set.insert(row_i_index);
 
        if (!row_index_set_ret.second)
        {
          std::ostringstream err_msg;
          err_msg << "The row " << row_i_index << " is already inserted.\n";
          throw OomphLibError(
            err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Stricter test - can be removed is required in the future. For the first
      // row, check that the column indices does not repeat.
      std::set<unsigned> col_index_set;
      for (unsigned col_i = 0; col_i < para_selected_block_ncol; col_i++)
      {
        std::pair<std::set<unsigned>::iterator, bool> col_index_set_ret;
 
        const unsigned col_i_index = selected_block[0][col_i].column_index();
 
        col_index_set_ret = col_index_set.insert(col_i_index);
 
        if (!col_index_set_ret.second)
        {
          std::ostringstream err_msg;
          err_msg << "The col " << col_i_index << " is already inserted.\n";
          throw OomphLibError(
            err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
      }
 
      // Loop through all the block_pt and check:
      // 1) The non-null pointers point to built matrices.
      // 2) The distribution matches those defined by selected_block within
      //    Block_distribution_pt.
      for (unsigned block_i = 0; block_i < para_selected_block_nrow; block_i++)
      {
        for (unsigned block_j = 0; block_j < para_selected_block_ncol;
             block_j++)
        {
          const CRDoubleMatrix* tmp_block_pt =
            selected_block[block_i][block_j].replacement_block_pt();
 
          if (tmp_block_pt != 0)
          {
            if (!tmp_block_pt->built())
            {
              std::ostringstream err_msg;
              err_msg << "The matrix pointed to by block(" << block_i << ","
                      << block_j << ") is not built.\n";
              throw OomphLibError(err_msg.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
 
            const LinearAlgebraDistribution* const tmp_block_dist_pt =
              tmp_block_pt->distribution_pt();
 
            const unsigned row_selected_block =
              selected_block[block_i][block_j].row_index();
 
            const LinearAlgebraDistribution* const another_tmp_block_dist_pt =
              block_distribution_pt(row_selected_block);
 
            if (*tmp_block_dist_pt != *another_tmp_block_dist_pt)
            {
              std::ostringstream err_msg;
              err_msg << "block_distribution_pt " << row_selected_block << "\n"
                      << "does not match the distribution from the block_pt() "
                      << " in selected_block[" << block_i << "][" << block_j
                      << "].\n";
              throw OomphLibError(err_msg.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
 
      // Attempt a similar check for the column index. This is not as rigorous
      // since a CRDoubleMatrix does not have a distribution for the columns.
      // However, we can check if the ncol of the matrices in block_pt matches
      // those in the block_distribution_pt corresponding to the columns.
      // (I hope this makes sense... both the row and columns are permuted in
      // CRDoubleMatrixHelpers::concatenate_without_communication(...))
      //
      // The test for the row distributions checks if the nrow_local is correct.
      // We do not have the equivalent for columns.
      for (unsigned block_i = 0; block_i < para_selected_block_nrow; block_i++)
      {
        for (unsigned block_j = 0; block_j < para_selected_block_ncol;
             block_j++)
        {
          // Cache the block_pt
          const CRDoubleMatrix* tmp_block_pt =
            selected_block[block_i][block_j].replacement_block_pt();
 
          if (tmp_block_pt != 0)
          {
            const unsigned tmp_block_ncol = tmp_block_pt->ncol();
 
            const unsigned selected_block_col =
              selected_block[block_i][block_j].column_index();
 
            // YES, nrow, this is not incorrect.
            const unsigned another_tmp_block_ncol =
              block_distribution_pt(selected_block_col)->nrow();
 
            if (tmp_block_ncol != another_tmp_block_ncol)
            {
              std::ostringstream err_msg;
              err_msg << "block_pt in selected_block[" << block_i << "]["
                      << block_j << "] "
                      << "has ncol = " << tmp_block_ncol << ".\n"
                      << "But the corresponding block_distribution_pt("
                      << selected_block_col
                      << ") has nrow = " << another_tmp_block_ncol << ").\n";
              throw OomphLibError(err_msg.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
#endif
 
      // The return matrix.
      MATRIX output_matrix;
 
      // How many sub matrices are there in the row and column?
      const unsigned nblock_row = selected_block.nrow();
      const unsigned nblock_col = selected_block.ncol();
 
      // Get the row and col distributions, this is required for concatenation
      // without communication.
      Vector<LinearAlgebraDistribution*> row_dist_pt(nblock_row, 0);
      Vector<LinearAlgebraDistribution*> col_dist_pt(nblock_col, 0);
 
      // For the row distributions, use the first column of selected_block
      // Also, store the index of the block rows.
      Vector<unsigned> block_row_index(nblock_row, 0);
      for (unsigned row_i = 0; row_i < nblock_row; row_i++)
      {
        const unsigned selected_row_index =
          selected_block[row_i][0].row_index();
 
        row_dist_pt[row_i] = Block_distribution_pt[selected_row_index];
        block_row_index[row_i] = selected_row_index;
      }
 
      // For the col distributions, use the first row of selected_block
      Vector<unsigned> block_col_index(nblock_col, 0);
      for (unsigned col_i = 0; col_i < nblock_col; col_i++)
      {
        const unsigned selected_col_index =
          selected_block[0][col_i].column_index();
 
        col_dist_pt[col_i] = Block_distribution_pt[selected_col_index];
        block_col_index[col_i] = selected_col_index;
      }
 
      // Now build the output matrix. The output_matrix needs a distribution,
      // this distribution is a concatenation of the block rows. But because
      // concatenation of distributions requires communication, we try to
      // minimise this process by creating it once, then store a key to the
      // concatenated distribution. First check to see if the block row indices
      // is already a key in Auxiliary_block_distribution_pt, if it is in there,
      // we use the distribution it corresponds to. Otherwise, we create the
      // distribution and store it for possible further use.
      std::map<Vector<unsigned>, LinearAlgebraDistribution*>::const_iterator
        iter;
 
      iter = Auxiliary_block_distribution_pt.find(block_row_index);
 
      if (iter != Auxiliary_block_distribution_pt.end())
      {
        output_matrix.build(iter->second);
      }
      else
      {
        LinearAlgebraDistribution* tmp_dist_pt = new LinearAlgebraDistribution;
        LinearAlgebraDistributionHelpers::concatenate(row_dist_pt,
                                                      *tmp_dist_pt);
        insert_auxiliary_block_distribution(block_row_index, tmp_dist_pt);
        output_matrix.build(tmp_dist_pt);
      }
 
      // Do the same for the column dist, since we might need it for the RHS
      // vector..
      iter = Auxiliary_block_distribution_pt.find(block_col_index);
      if (iter == Auxiliary_block_distribution_pt.end())
      {
        LinearAlgebraDistribution* tmp_dist_pt = new LinearAlgebraDistribution;
        LinearAlgebraDistributionHelpers::concatenate(col_dist_pt,
                                                      *tmp_dist_pt);
        insert_auxiliary_block_distribution(block_col_index, tmp_dist_pt);
      }
 
      // Storage for the pointers to CRDoubleMatrices to concatenate.
      DenseMatrix<CRDoubleMatrix*> block_pt(nblock_row, nblock_col, 0);
 
      // Vector of Vectors of unsigns to indicate whether we have created
      // CRDoubleMatrices with new or not... so we can delete it later.
      // 0 - no new CRDoubleMatrix is created.
      // 1 - a new CRDoubleMatrix is created.
      // If we ever use C++11, remove this and use smart pointers.
      Vector<Vector<unsigned>> new_block(nblock_row,
                                         Vector<unsigned>(nblock_col, 0));
 
      // Get blocks if wanted.
      for (unsigned block_i = 0; block_i < nblock_row; block_i++)
      {
        for (unsigned block_j = 0; block_j < nblock_col; block_j++)
        {
          const bool block_wanted = selected_block[block_i][block_j].wanted();
 
          if (block_wanted)
          {
            CRDoubleMatrix* tmp_block_pt =
              selected_block[block_i][block_j].replacement_block_pt();
 
            if (tmp_block_pt == 0)
            {
              new_block[block_i][block_j] = 1;
 
              block_pt(block_i, block_j) = new CRDoubleMatrix;
 
              // temp variables for readability purposes.
              const unsigned tmp_block_i = block_row_index[block_i];
              const unsigned tmp_block_j = block_col_index[block_j];
 
              // Get the block.
              this->get_block(
                tmp_block_i, tmp_block_j, *block_pt(block_i, block_j));
            }
            else
            {
              block_pt(block_i, block_j) = tmp_block_pt;
            }
          }
        }
      }
 
      // Perform the concatenation.
      CRDoubleMatrixHelpers::concatenate_without_communication(
        row_dist_pt, col_dist_pt, block_pt, output_matrix);
 
      // Delete any new CRDoubleMatrices we created.
      for (unsigned block_i = 0; block_i < nblock_row; block_i++)
      {
        for (unsigned block_j = 0; block_j < nblock_col; block_j++)
        {
          if (new_block[block_i][block_j])
          {
            delete block_pt(block_i, block_j);
          }
        }
      }
 
      return output_matrix;
    } // EOFunc get_concatenated_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;
 
    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
    {
#ifdef PARANOID
      if (Block_to_dof_map_coarse.size() == 0)
      {
        std::ostringstream error_msg;
        error_msg
          << "The Block_to_dof_map_coarse vector is not setup for \n"
          << "this block preconditioner.\n\n"
 
          << "This vector is always set up within block_setup(...).\n"
          << "If block_setup() is already called, then perhaps there is\n"
          << "something wrong with your block preconditionable elements.\n"
          << std::endl;
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Return the number of block types.
      return Block_to_dof_map_coarse.size();
    } // EOFunc nblock_types(...)
 
    unsigned ndof_types() const
    {
#ifdef PARANOID
      // Subsidiary preconditioners don't really need the meshes
      if (this->is_master_block_preconditioner())
      {
        std::ostringstream err_msg;
        unsigned n = nmesh();
        if (n == 0)
        {
          err_msg << "No meshes have been set for this block preconditioner!\n"
                  << "Set one with set_nmesh(...), set_mesh(...)" << std::endl;
          throw OomphLibError(
            err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
          for (unsigned m = 0; m < n; m++)
          {
            if (Mesh_pt[m] == 0)
            {
              err_msg << "The mesh pointer to mesh " << m << " is null!\n"
                      << "Set a non-null one with set_mesh(...)" << std::endl;
              throw OomphLibError(err_msg.str(),
                                  OOMPH_CURRENT_FUNCTION,
                                  OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
      }
#endif
 
      // If this is a subsidiary block preconditioner, then the function
      // turn_into_subsidiary_block_preconditioner(...) should have been called,
      // this function would have set up the look up lists between coarsened
      // dof types and the internal dof types. Of coarse, the user (the writer
      // of the preconditioners) should not care about the internal dof types
      // and what happens under the hood. Thus they should get the number of
      // coarsened dof types (i.e. the number of dof types the preconditioner
      // above (parent preconditioner) decides to give to this preconditioner).
      if (is_subsidiary_block_preconditioner())
      {
#ifdef PARANOID
        if (Doftype_coarsen_map_coarse.size() == 0)
        {
          std::ostringstream error_msg;
          error_msg
            << "The Doftype_coarsen_map_coarse vector is not setup for \n"
            << "this SUBSIDIARY block preconditioner.\n\n"
 
            << "For SUBSIDARY block preconditioners at any level, this\n"
            << "vector is set up in the function \n"
            << "turn_into_subsidiary_block_preconditioner(...).\n\n"
 
            << "Being a SUBSIDIARY block preconditioner means that \n"
            << "(Master_block_preconditioner_pt == 0) is true.\n"
            << "The Master_block_preconditioner_pt MUST be set in the "
            << "function \n"
            << "turn_into_subsidiary_block_preconditioner(...).\n\n"
 
            << "Somewhere, the Master_block_preconditioner_pt pointer is\n"
            << "set but not by the function\n"
            << "turn_into_subsidiary_block_preconditioner(...).\n"
            << std::endl;
          throw OomphLibError(
            error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
#endif
        // return the number of dof types.
        return Doftype_coarsen_map_coarse.size();
      }
      else
      // Otherwise the number of ndof types is the same as the internal number
      // of dof types, since no coarsening of the dof types is done at the
      // top-most master level.
      {
        return internal_ndof_types();
      }
    } // EOFunc ndof_types(...)
 
 
    const Mesh* mesh_pt(const unsigned& i) const
    {
#ifdef PARANOID
      if (is_subsidiary_block_preconditioner())
      {
        std::ostringstream error_msg;
        error_msg << "The mesh_pt() function should not be called on a\n"
                  << "subsidiary block preconditioner." << std::endl;
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      const Mesh* mesh_i_pt = Mesh_pt[i];
 
#ifdef PARANOID
      if (mesh_i_pt == 0)
      {
        std::ostringstream error_msg;
        error_msg << "Mesh pointer " << mesh_i_pt << " is null.";
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      return mesh_i_pt;
    } // EOFunc mesh_pt(...)
 
    unsigned nmesh() const
    {
      return Mesh_pt.size();
    } // EOFunc nmesh()
 
    int block_number(const unsigned& i_dof) const
    {
      int internal_block_number = this->internal_block_number(i_dof);
 
      if (internal_block_number == -1)
      {
        return internal_block_number;
      }
      else
      {
        // Map the internal block to the "external" block number, i.e. what the
        // writer of the preconditioner is expects.
        unsigned block_i = 0;
        while (std::find(Block_to_dof_map_fine[block_i].begin(),
                         Block_to_dof_map_fine[block_i].end(),
                         internal_block_number) ==
               Block_to_dof_map_fine[block_i].end())
        {
          block_i++;
        }
 
        return block_i;
      }
    }
 
    int index_in_block(const unsigned& i_dof) const
    {
      // the dof block number
      int internal_dof_block_number = this->internal_dof_number(i_dof);
 
      if (internal_dof_block_number >= 0)
      {
        // the external block number
        unsigned ex_blk_number = this->block_number(i_dof);
 
        int internal_index_in_dof = this->internal_index_in_dof(i_dof);
 
        // find the processor which this global index in block belongs to.
        unsigned block_proc =
          internal_block_distribution_pt(internal_dof_block_number)
            ->rank_of_global_row(internal_index_in_dof);
 
        // Add up all of the first rows.
        const unsigned ndof_in_block =
          Block_to_dof_map_fine[ex_blk_number].size();
 
        unsigned index = 0;
        for (unsigned dof_i = 0; dof_i < ndof_in_block; dof_i++)
        {
          index += internal_block_distribution_pt(
                     Block_to_dof_map_fine[ex_blk_number][dof_i])
                     ->first_row(block_proc);
        }
 
        // Now add up all the nrow_local up to this dof block.
        unsigned j = 0;
 
        while (int(Block_to_dof_map_fine[ex_blk_number][j]) !=
               internal_dof_block_number)
        {
          index += internal_block_distribution_pt(
                     Block_to_dof_map_fine[ex_blk_number][j])
                     ->nrow_local(block_proc);
          j++;
        }
 
        // Now add the index of this block...
        index += (internal_index_in_dof -
                  internal_block_distribution_pt(internal_dof_block_number)
                    ->first_row(block_proc));
 
        return index;
      }
 
      return -1;
    }
 
    const LinearAlgebraDistribution* block_distribution_pt(
      const unsigned& b) const
    {
#ifdef PARANOID
      if (Block_distribution_pt.size() == 0)
      {
        std::ostringstream error_msg;
        error_msg << "Block distributions are not set up.\n"
                  << "Have you called block_setup(...)?\n"
                  << std::endl;
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
      if (b > nblock_types())
      {
        std::ostringstream error_msg;
        error_msg << "You requested the distribution for the block " << b
                  << ".\n"
                  << "But there are only " << nblock_types()
                  << " block types.\n"
                  << std::endl;
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      return Block_distribution_pt[b];
    } // EOFunc block_distribution_pt(...)
 
    LinearAlgebraDistribution* block_distribution_pt(const unsigned b)
    {
#ifdef PARANOID
      if (Block_distribution_pt.size() == 0)
      {
        std::ostringstream error_msg;
        error_msg << "Block distributions are not set up.\n"
                  << "Have you called block_setup(...)?\n"
                  << std::endl;
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
      if (b > nblock_types())
      {
        std::ostringstream error_msg;
        error_msg << "You requested the distribution for the block " << b
                  << ".\n"
                  << "But there are only " << nblock_types()
                  << " block types.\n"
                  << std::endl;
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      return Block_distribution_pt[b];
    } // EOFunc block_distribution_pt(...)
 
    LinearAlgebraDistribution* dof_block_distribution_pt(const unsigned& b)
    {
#ifdef PARANOID
      if (b > ndof_types())
      {
        std::ostringstream error_msg;
        error_msg << "You requested the distribution for the dof block " << b
                  << ".\n"
                  << "But there are only " << ndof_types() << " DOF types.\n"
                  << std::endl;
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
#endif
 
      if (is_master_block_preconditioner())
      {
#ifdef PARANOID
        if (Internal_block_distribution_pt.size() == 0)
        {
          std::ostringstream error_msg;
          error_msg << "Internal block distributions are not set up.\n"
                    << "Have you called block_setup(...)?\n"
                    << std::endl;
          throw OomphLibError(
            error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
#endif
 
        // The dof block is distribution is the same as the internal
        // block distribution.
        return Internal_block_distribution_pt[b];
      }
      else
      {
#ifdef PARANOID
        if (Dof_block_distribution_pt.size() == 0)
        {
          std::ostringstream error_msg;
          error_msg << "Dof block distributions are not set up.\n"
                    << "Have you called block_setup(...)?\n"
                    << std::endl;
          throw OomphLibError(
            error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
#endif
        return Dof_block_distribution_pt[b];
      }
    } // EOFunc block_distribution_pt(...)
 
 
    const LinearAlgebraDistribution* master_distribution_pt() const
    {
      if (is_master_block_preconditioner())
      {
        return this->distribution_pt();
      }
      else
      {
        return Master_block_preconditioner_pt->master_distribution_pt();
      }
    } // EOFunc master_distribution_pt(...)
 
    unsigned ndof_types_in_mesh(const unsigned& i) const
    {
#ifdef PARANOID
      if (is_subsidiary_block_preconditioner())
      {
        std::ostringstream err_msg;
        err_msg << "A subsidiary block preconditioner should not care about\n"
                << "anything to do with meshes.";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      if (Ndof_types_in_mesh.size() == 0)
      {
        return mesh_pt(i)->ndof_types();
      }
      else
      {
        return Ndof_types_in_mesh[i];
      }
    } // EOFunc ndof_types_in_mesh(...)
 
    bool is_subsidiary_block_preconditioner() const
    {
      return (this->Master_block_preconditioner_pt != 0);
    } // EOFunc is_subsidiary_block_preconditioner()
 
    bool is_master_block_preconditioner() const
    {
      return (this->Master_block_preconditioner_pt == 0);
    } // EOFunc is_master_block_preconditioner()
 
    void set_block_output_to_files(const std::string& basefilename)
    {
      Output_base_filename = basefilename;
    } // EOFunc set_block_output_to_files(...)
 
    void disable_block_output_to_files()
    {
      Output_base_filename.clear();
    } // EOFunc disable_block_output_to_files()
 
    bool block_output_on() const
    {
      return Output_base_filename.size() > 0;
    } // EOFunc block_output_on()
 
    void output_blocks_to_files(const std::string& basefilename,
                                const unsigned& precision = 8) const
    {
      unsigned nblocks = internal_nblock_types();
 
      for (unsigned i = 0; i < nblocks; i++)
      {
        for (unsigned j = 0; j < nblocks; j++)
        {
          // Construct the filename.
          std::string filename(basefilename + "_block_" +
                               StringConversion::to_string(i) + "_" +
                               StringConversion::to_string(j));
 
          // Write out the block.
          get_block(i, j).sparse_indexed_output(filename, precision, true);
        }
      }
    } // EOFunc output_blocks_to_files(...)
 
    void post_block_matrix_assembly_partial_clear()
    {
      if (is_master_block_preconditioner())
      {
        Index_in_dof_block_dense.clear();
        Dof_number_dense.clear();
#ifdef OOMPH_HAS_MPI
        Index_in_dof_block_sparse.clear();
        Dof_number_sparse.clear();
        Global_index_sparse.clear();
        Index_in_dof_block_sparse.clear();
        Dof_number_sparse.clear();
#endif
        Dof_dimension.clear();
      }
      Ndof_in_block.clear();
      Dof_number_to_block_number_lookup.clear();
      Block_number_to_dof_number_lookup.clear();
    } // EOFunc post_block_matrix_assembly_partial_clear()
 
    BlockPreconditioner<MATRIX>* master_block_preconditioner_pt() const
    {
#ifdef PARANOID
      if (is_master_block_preconditioner())
      {
        std::ostringstream error_message;
        error_message << "This block preconditioner does not have "
                      << "a master preconditioner.";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Master_block_preconditioner_pt;
    } // EOFunc master_block_preconditioner_pt()
 
    void clear_block_preconditioner_base()
    {
      Replacement_dof_block_pt.clear();
 
      // clear the Distributions
      this->clear_distribution();
      unsigned nblock = Internal_block_distribution_pt.size();
      for (unsigned b = 0; b < nblock; b++)
      {
        delete Internal_block_distribution_pt[b];
      }
      Internal_block_distribution_pt.resize(0);
 
      // clear the global index
      Global_index.clear();
 
      // call the post block matrix assembly clear
      this->post_block_matrix_assembly_partial_clear();
 
#ifdef OOMPH_HAS_MPI
      // storage if the matrix is distributed
      unsigned nr = Rows_to_send_for_get_block.nrow();
      unsigned nc = Rows_to_send_for_get_block.ncol();
      for (unsigned p = 0; p < nc; p++)
      {
        delete[] Rows_to_send_for_get_ordered[p];
        delete[] Rows_to_recv_for_get_ordered[p];
        for (unsigned b = 0; b < nr; b++)
        {
          delete[] Rows_to_recv_for_get_block(b, p);
          delete[] Rows_to_send_for_get_block(b, p);
        }
      }
      Rows_to_recv_for_get_block.resize(0, 0);
      Nrows_to_recv_for_get_block.resize(0, 0);
      Rows_to_send_for_get_block.resize(0, 0);
      Nrows_to_send_for_get_block.resize(0, 0);
      Rows_to_recv_for_get_ordered.clear();
      Nrows_to_recv_for_get_ordered.clear();
      Rows_to_send_for_get_ordered.clear();
      Nrows_to_send_for_get_ordered.clear();
 
#endif
 
      // zero
      if (is_master_block_preconditioner())
      {
        Nrow = 0;
        Internal_ndof_types = 0;
        Internal_nblock_types = 0;
      }
 
      // delete the prec matrix dist pt
      delete Internal_preconditioner_matrix_distribution_pt;
      Internal_preconditioner_matrix_distribution_pt = 0;
      delete Preconditioner_matrix_distribution_pt;
      Preconditioner_matrix_distribution_pt = 0;
 
      // Delete any existing (external) block distributions.
      const unsigned n_existing_block_dist = Block_distribution_pt.size();
      for (unsigned dist_i = 0; dist_i < n_existing_block_dist; dist_i++)
      {
        delete Block_distribution_pt[dist_i];
      }
 
      // Clear the vector.
      Block_distribution_pt.clear();
 
 
      // Create the identity key.
      Vector<unsigned> preconditioner_matrix_key(n_existing_block_dist, 0);
      for (unsigned i = 0; i < n_existing_block_dist; i++)
      {
        preconditioner_matrix_key[i] = i;
      }
 
      // Now iterate through Auxiliary_block_distribution_pt
      // and delete all distributions, except for the one which corresponds
      // to the identity since this is already deleted.
      std::map<Vector<unsigned>, LinearAlgebraDistribution*>::iterator iter =
        Auxiliary_block_distribution_pt.begin();
 
      while (iter != Auxiliary_block_distribution_pt.end())
      {
        if (iter->first != preconditioner_matrix_key)
        {
          delete iter->second;
          iter++;
        }
        else
        {
          ++iter;
        }
      }
 
      // Now clear it.
      Auxiliary_block_distribution_pt.clear();
 
      // Delete any dof block distributions
      const unsigned ndof_block_dist = Dof_block_distribution_pt.size();
      for (unsigned dof_i = 0; dof_i < ndof_block_dist; dof_i++)
      {
        delete Dof_block_distribution_pt[dof_i];
      }
      Dof_block_distribution_pt.clear();
 
    } // EOFunc clear_block_preconditioner_base()
 
    void document()
    {
      oomph_info << std::endl;
      oomph_info << "===========================================" << std::endl;
      oomph_info << "Block Preconditioner Documentation" << std::endl
                 << std::endl;
      oomph_info << "Number of DOF types: " << internal_ndof_types()
                 << std::endl;
      oomph_info << "Number of block types: " << internal_nblock_types()
                 << std::endl;
      oomph_info << std::endl;
      if (is_subsidiary_block_preconditioner())
      {
        for (unsigned d = 0; d < Internal_ndof_types; d++)
        {
          oomph_info << "Master DOF number " << d << " : "
                     << this->internal_master_dof_number(d) << std::endl;
        }
      }
      oomph_info << std::endl;
      for (unsigned b = 0; b < internal_nblock_types(); b++)
      {
        oomph_info << "Block " << b << " DOF types:";
        for (unsigned i = 0; i < Block_number_to_dof_number_lookup[b].size();
             i++)
        {
          oomph_info << " " << Block_number_to_dof_number_lookup[b][i];
        }
        oomph_info << std::endl;
      }
      oomph_info << std::endl;
      oomph_info << "Master block preconditioner distribution:" << std::endl;
      oomph_info << *master_distribution_pt() << std::endl;
      oomph_info << "Internal preconditioner matrix distribution:" << std::endl;
      oomph_info << *internal_preconditioner_matrix_distribution_pt()
                 << std::endl;
      oomph_info << "Preconditioner matrix distribution:" << std::endl;
      oomph_info << *preconditioner_matrix_distribution_pt() << std::endl;
      for (unsigned b = 0; b < Internal_nblock_types; b++)
      {
        oomph_info << "Internal block " << b << " distribution:" << std::endl;
        oomph_info << *Internal_block_distribution_pt[b] << std::endl;
      }
      for (unsigned b = 0; b < nblock_types(); b++)
      {
        oomph_info << "Block " << b << " distribution:" << std::endl;
        oomph_info << *Block_distribution_pt[b] << std::endl;
      }
 
      // DS: the functions called here no longer exist and this function is
      // never used as far as I can tell, so it should be fine to comment this
      // bit out:
      // if (is_master_block_preconditioner())
      //  {
      //   oomph_info << "First look-up row: " << this->first_lookup_row()
      //              << std::endl;
      //   oomph_info << "Number of look-up rows: "
      //              << this->nlookup_rows() << std::endl;
      //  }
      oomph_info << "===========================================" << std::endl;
      oomph_info << std::endl;
    } // EOFunc document()
 
    Vector<Vector<unsigned>> doftype_coarsen_map_fine() const
    {
      return Doftype_coarsen_map_fine;
    }
 
    Vector<unsigned> get_fine_grain_dof_types_in(const unsigned& i) const
    {
#ifdef PARANOID
      const unsigned n_dof_types = ndof_types();
 
      if (i >= n_dof_types)
      {
        std::ostringstream err_msg;
        err_msg << "Trying to get the most fine grain dof types in dof type "
                << i << ",\nbut there are only " << n_dof_types
                << " number of dof types.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Doftype_coarsen_map_fine[i];
    }
 
    unsigned nfine_grain_dof_types_in(const unsigned& i) const
    {
#ifdef PARANOID
      const unsigned n_dof_types = ndof_types();
 
      if (i >= n_dof_types)
      {
        std::ostringstream err_msg;
        err_msg << "Trying to get the number of most fine grain dof types "
                << "in dof type " << i << ",\n"
                << "but there are only " << n_dof_types
                << " number of dof types.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Doftype_coarsen_map_fine[i].size();
    }
 
    MapMatrix<unsigned, CRDoubleMatrix*> replacement_dof_block_pt() const
    {
      return Replacement_dof_block_pt;
    } // EOFunc replacement_block_pt()
 
    void setup_matrix_vector_product(MatrixVectorProduct* matvec_prod_pt,
                                     CRDoubleMatrix* block_pt,
                                     const Vector<unsigned>& block_col_indices)
    {
      const unsigned nblock = block_col_indices.size();
 
      if (nblock == 1)
      {
        const unsigned col_index = block_col_indices[0];
        matvec_prod_pt->setup(block_pt, Block_distribution_pt[col_index]);
      }
      else
      {
        std::map<Vector<unsigned>, LinearAlgebraDistribution*>::const_iterator
          iter;
 
        iter = Auxiliary_block_distribution_pt.find(block_col_indices);
        if (iter != Auxiliary_block_distribution_pt.end())
        {
          matvec_prod_pt->setup(block_pt, iter->second);
        }
        else
        {
          Vector<LinearAlgebraDistribution*> tmp_vec_dist_pt(nblock, 0);
          for (unsigned b = 0; b < nblock; b++)
          {
            tmp_vec_dist_pt[b] = Block_distribution_pt[block_col_indices[b]];
          }
 
          LinearAlgebraDistribution* tmp_dist_pt =
            new LinearAlgebraDistribution;
          LinearAlgebraDistributionHelpers::concatenate(tmp_vec_dist_pt,
                                                        *tmp_dist_pt);
          insert_auxiliary_block_distribution(block_col_indices, tmp_dist_pt);
          matvec_prod_pt->setup(block_pt, tmp_dist_pt);
        }
      }
    } // EOFunc setup_matrix_vector_product(...)
 
    void setup_matrix_vector_product(MatrixVectorProduct* matvec_prod_pt,
                                     CRDoubleMatrix* block_pt,
                                     const unsigned& block_col_index)
    {
      Vector<unsigned> col_index_vector(1, block_col_index);
      setup_matrix_vector_product(matvec_prod_pt, block_pt, col_index_vector);
    } // EOFunc setup_matrix_vector_product(...)
 
    // private:
 
    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
    {
#ifdef PARANOID
      if (Internal_nblock_types == 0)
      {
        std::ostringstream err_msg;
        err_msg
          << "(Internal_nblock_types == 0) is true. \n"
          << "Did you remember to call the function block_setup(...)?\n\n"
 
          << "This variable is always set up within block_setup(...).\n"
          << "If block_setup() is already called, then perhaps there is\n"
          << "something wrong with your block preconditionable elements.\n"
          << std::endl;
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      if (Internal_nblock_types != internal_ndof_types())
      {
        std::ostringstream err_msg;
        err_msg
          << "The number of internal block types and "
          << "internal dof types does not match... \n\n"
          << "Internally, the number of block types and the number of dof "
          << "types must be the same.\n"
          << std::endl;
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // return the number of internal block types.
      return Internal_nblock_types;
    } // EOFunc internal_nblock_types(...)
 
    unsigned internal_ndof_types() const
    {
      if (is_subsidiary_block_preconditioner())
      // If this is a subsidiary block preconditioner, then the variable
      // Internal_ndof_types must always be set up.
      {
#ifdef PARANOID
        if (Internal_ndof_types == 0)
        {
          std::ostringstream error_msg;
          error_msg
            << "(Internal_ndof_types == 0) is true.\n"
            << "This means that the Master_block_preconditioner_pt pointer is\n"
            << "set but possibly not by the function\n"
            << "turn_into_subsidiary_block_preconditioner(...).\n\n"
 
            << "This goes against the block preconditioning framework "
            << "methodology.\n"
            << "Many machinery relies on the look up lists set up by the \n"
            << "function turn_into_subsidiary_block_preconditioner(...) \n"
            << "between the parent and child block preconditioners.\n"
            << std::endl;
          throw OomphLibError(
            error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
        }
#endif
        return Internal_ndof_types;
      }
      else
      // Else, this is a master block preconditioner, calculate the number of
      // dof types from the meshes.
      {
        unsigned ndof = 0;
        for (unsigned i = 0; i < nmesh(); i++)
        {
          ndof += ndof_types_in_mesh(i);
        }
        return ndof;
      }
    } // EOFunc internal_ndof_types(...)
 
    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,
                            MATRIX& output_block) const;
 
    int internal_block_number(const unsigned& i_dof) const
    {
      int dn = internal_dof_number(i_dof);
      if (dn == -1)
      {
        return dn;
      }
      else
      {
        return Dof_number_to_block_number_lookup[dn];
      }
    } // EOFunc internal_block_number(...)
 
    int internal_index_in_block(const unsigned& i_dof) const
    {
      // the index in the dof block
      unsigned index = internal_index_in_dof(i_dof);
 
      // the dof block number
      int internal_dof_block_number = internal_dof_number(i_dof);
      if (internal_dof_block_number >= 0)
      {
        // the 'actual' block number
        unsigned blk_number = internal_block_number(i_dof);
 
        // compute the index in the block
        unsigned j = 0;
        while (int(Block_number_to_dof_number_lookup[blk_number][j]) !=
               internal_dof_block_number)
        {
          index += internal_dof_block_dimension(
            Block_number_to_dof_number_lookup[blk_number][j]);
          j++;
        }
 
        // and return
        return index;
      }
      return -1;
    } // EOFunc internal_index_in_block(...)
 
    const LinearAlgebraDistribution* internal_block_distribution_pt(
      const unsigned& b) const
    {
#ifdef PARANOID
      if (Internal_block_distribution_pt.size() == 0)
      {
        std::ostringstream error_msg;
        error_msg << "Internal block distributions are not set up.\n"
                  << "Have you called block_setup(...)?\n"
                  << std::endl;
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
      if (b > internal_nblock_types())
      {
        std::ostringstream error_msg;
        error_msg << "You requested the distribution for the internal block "
                  << b << ".\n"
                  << "But there are only " << internal_nblock_types()
                  << " block types.\n"
                  << std::endl;
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Internal_block_distribution_pt[b];
    } // EOFunc internal_block_distribution_pt(...)
 
    void insert_auxiliary_block_distribution(
      const Vector<unsigned>& block_vec_number,
      LinearAlgebraDistribution* dist_pt)
    {
#ifdef PARANOID
      const unsigned max_block_number =
        *std::max_element(block_vec_number.begin(), block_vec_number.end());
 
      const unsigned nblocks = nblock_types();
      if (max_block_number >= nblocks)
      {
        std::ostringstream err_msg;
        err_msg << "Cannot insert into Auxiliary_block_distribution_pt\n"
                << "because " << max_block_number << " is equal to or \n"
                << "greater than " << nblocks << ".\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Now check if the pair already exists in
      // Auxiliary_block_distribution_pt. This is a stricter test and can be
      // removed if required.
 
      // Attempt to get an iterator pointing to the pair with the value
      // block_vec_number.
      std::map<Vector<unsigned>, LinearAlgebraDistribution*>::const_iterator
        iter = Auxiliary_block_distribution_pt.find(block_vec_number);
 
      if (iter != Auxiliary_block_distribution_pt.end())
      // If it exists, we throw an error
      {
        std::ostringstream err_msg;
        err_msg << "Cannot insert into Auxiliary_block_distribution_pt\n"
                << "because the first in the pair already exists.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      Auxiliary_block_distribution_pt.insert(
        std::make_pair(block_vec_number, dist_pt));
    } // insert_auxiliary_block_distribution(...)
 
    void block_matrix_test(const unsigned& i,
                           const unsigned& j,
                           const MATRIX* block_matrix_pt) const;
 
    template<typename myType>
    inline int get_index_of_value(const Vector<myType>& vec,
                                  const myType val,
                                  const bool sorted = false) const
    {
      if (sorted)
      {
        typename Vector<myType>::const_iterator low =
          std::lower_bound(vec.begin(), vec.end(), val);
 
        return (low == vec.end() || *low != val) ? -1 : (low - vec.begin());
      }
      else
      {
        int pos = std::find(vec.begin(), vec.end(), val) - vec.begin();
        return (pos < int(vec.size()) && pos >= 0) ? pos : -1;
      }
    }
 
  private:
  protected:
    void set_nmesh(const unsigned& n)
    {
      Mesh_pt.resize(n, 0);
      Allow_multiple_element_type_in_mesh.resize(n, 0);
    } // EOFunc set_nmesh(...)
 
 
    void set_mesh(const unsigned& i,
                  const Mesh* const mesh_pt,
                  const bool& allow_multiple_element_type_in_mesh = false)
    {
#ifdef PARANOID
      // paranoid check that mesh i can be set
      if (i >= nmesh())
      {
        std::ostringstream err_msg;
        err_msg << "The mesh pointer has space for " << nmesh() << " meshes.\n"
                << "Cannot store a mesh at entry " << i << "\n"
                << "Has set_nmesh(...) been called?";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Check that the mesh pointer is not null.
      if (mesh_pt == 0)
      {
        std::ostringstream err_msg;
        err_msg << "Tried to set the " << i
                << "-th mesh pointer, but it is null.";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // store the mesh pt and n dof types
      Mesh_pt[i] = mesh_pt;
 
      // Does this mesh contain multiple element types?
      Allow_multiple_element_type_in_mesh[i] =
        unsigned(allow_multiple_element_type_in_mesh);
    } // EOFunc set_mesh(...)
 
 
    void set_replacement_dof_block(const unsigned& block_i,
                                   const unsigned& block_j,
                                   CRDoubleMatrix* replacement_dof_block_pt)
    {
#ifdef PARANOID
      // Check if block_setup(...) has been called.
      if (nblock_types() == 0)
      {
        std::ostringstream err_msg;
        err_msg << "nblock_types() is 0, has block_setup(...) been called?\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
 
      // Range checking for replacement dof block.
      unsigned para_ndof_types = this->ndof_types();
 
      if ((block_i >= para_ndof_types) || (block_j >= para_ndof_types))
      {
        std::ostringstream err_msg;
        err_msg << "Replacement dof block (" << block_i << "," << block_j
                << ") is outside of range:\n"
                << para_ndof_types;
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
 
      //    // Check that the most fine grain mapping has been used in
      //    block_setup(...)
      //    // i.e. nblock_types() == ndof_types()
      //    if(ndof_types() != nblock_types())
      //    {
      //      std::ostringstream err_msg;
      //      err_msg << "ndof_types() != nblock_types()\n"
      //              << "Only the dof-level blocks can be replaced.\n"
      //              << "Please re-think your blocking scheme.\n";
      //      throw OomphLibError(err_msg.str(),
      //          OOMPH_CURRENT_FUNCTION,
      //          OOMPH_EXCEPTION_LOCATION);
      //    }
 
      // Check that the replacement block pt is not null
      if (replacement_dof_block_pt == 0)
      {
        std::ostringstream err_msg;
        err_msg << "Replacing block(" << block_i << "," << block_i << ")\n"
                << " but the pointer is NULL." << std::endl;
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Check that the replacement block has been built
      if (!replacement_dof_block_pt->built())
      {
        std::ostringstream err_msg;
        err_msg << "Replacement block(" << block_i << "," << block_i << ")"
                << " is not built." << std::endl;
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Check if the distribution matches. Determine which natural ordering dof
      // this should go to. I.e. we convert from dof-block index to dof index.
      // Luckily, this is stored in Block_to_dof_map_coarse.
      //    const unsigned para_dof_block_i =
      //    Block_to_dof_map_coarse[block_i][0];
      const unsigned para_dof_block_i = block_i;
 
      if (*dof_block_distribution_pt(para_dof_block_i) !=
          *replacement_dof_block_pt->distribution_pt())
      {
        std::ostringstream err_msg;
        err_msg << "The distribution of the replacement dof_block_pt\n"
                << "is different from the Dof_block_distribution_pt["
                << para_dof_block_i << "].\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Now that we know the distribution of the replacement block is
      // correct, we check the number of columns.
      const unsigned para_dof_block_j = block_j;
      unsigned para_replacement_block_ncol = replacement_dof_block_pt->ncol();
      unsigned para_required_ncol =
        dof_block_distribution_pt(para_dof_block_j)->nrow();
      if (para_replacement_block_ncol != para_required_ncol)
      {
        std::ostringstream err_msg;
        err_msg << "Replacement dof block has ncol = "
                << para_replacement_block_ncol << ".\n"
                << "But required ncol is " << para_required_ncol << ".\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Block_to_dof_map_coarse[x][0] sense because we only can use this if
      // nblock_types() == ndof_types(), i.e. each sub-vector is of length 1.
      //
      // We use this indirection so that the placement of the pointer is
      // consistent with internal_get_block(...).
      //    const unsigned dof_block_i = Block_to_dof_map_coarse[block_i][0];
      //    const unsigned dof_block_j = Block_to_dof_map_coarse[block_j][0];
 
      //    Replacement_dof_block_pt(dof_block_i,dof_block_j)
      //      = replacement_dof_block_pt;
 
      Replacement_dof_block_pt(block_i, block_j) = replacement_dof_block_pt;
    }
 
    bool any_mesh_distributed() const
    {
#ifdef OOMPH_HAS_MPI
      // is_mesh_distributed() is only available with MPI
      for (unsigned i = 0, n = nmesh(); i < n; i++)
      {
        if (mesh_pt(i)->is_mesh_distributed())
        {
          return true;
        }
      }
#endif
      return false;
    }
 
    int internal_dof_number(const unsigned& i_dof) const
    {
      if (is_master_block_preconditioner())
      {
#ifdef OOMPH_HAS_MPI
        unsigned first_row = this->distribution_pt()->first_row();
        unsigned nrow_local = this->distribution_pt()->nrow_local();
        unsigned last_row = first_row + nrow_local - 1;
        if (i_dof >= first_row && i_dof <= last_row)
        {
          return static_cast<int>(Dof_number_dense[i_dof - first_row]);
        }
        else
        {
          // int index = this->get_index_of_element(Global_index_sparse,i_dof);
          int index =
            get_index_of_value<unsigned>(Global_index_sparse, i_dof, true);
          if (index >= 0)
          {
            return Dof_number_sparse[index];
          }
        }
        // if we here we couldn't find the i_dof
#ifdef PARANOID
        unsigned my_rank = comm_pt()->my_rank();
        std::ostringstream error_message;
        error_message << "Proc " << my_rank
                      << ": Requested internal_dof_number(...) for global DOF "
                      << i_dof << "\n"
                      << "cannot be found.\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
#endif
#else
        return static_cast<int>(Dof_number_dense[i_dof]);
#endif
      }
      // else this preconditioner  is a subsidiary one, and its Block_number
      // lookup schemes etc haven't been set up
      else
      {
        // Block number in master prec
        unsigned blk_num =
          Master_block_preconditioner_pt->internal_dof_number(i_dof);
 
        // Search through the Block_number_in_master_preconditioner for master
        // block blk_num and return the block number in this preconditioner
        for (unsigned i = 0; i < this->internal_ndof_types(); i++)
        {
          if (Doftype_in_master_preconditioner_fine[i] == blk_num)
          {
            return static_cast<int>(i);
          }
        }
        // if the master block preconditioner number is not found return -1
        return -1;
      }
 
      // Shouldn't get here
      throw OomphLibError(
        "Never get here\n", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      // Dummy return
      return -1;
    }
 
    unsigned internal_index_in_dof(const unsigned& i_dof) const
    {
      if (is_master_block_preconditioner())
      {
#ifdef OOMPH_HAS_MPI
        unsigned first_row = this->distribution_pt()->first_row();
        unsigned nrow_local = this->distribution_pt()->nrow_local();
        unsigned last_row = first_row + nrow_local - 1;
        if (i_dof >= first_row && i_dof <= last_row)
        {
          return static_cast<int>(Index_in_dof_block_dense[i_dof - first_row]);
        }
        else
        {
          // int index = this->get_index_of_element(Global_index_sparse,i_dof);
          int index =
            get_index_of_value<unsigned>(Global_index_sparse, i_dof, true);
          if (index >= 0)
          {
            return Index_in_dof_block_sparse[index];
          }
        }
        // if we here we couldn't find the i_dof
#ifdef PARANOID
        std::ostringstream error_message;
        error_message << "Requested internal_index_in_dof(...) for global DOF "
                      << i_dof << "\n"
                      << "cannot be found.\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
#endif
#else
        return Index_in_dof_block_dense[i_dof];
#endif
      }
      else
      {
        return Master_block_preconditioner_pt->internal_index_in_dof(i_dof);
      }
 
      // Shouldn't get here
      throw OomphLibError(
        "Never get here\n", OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      // Dummy return
      return -1;
    }
 
    unsigned internal_block_dimension(const unsigned& b) const
    {
#ifdef PARANOID
      const unsigned i_nblock_types = internal_nblock_types();
      if (b >= i_nblock_types)
      {
        std::ostringstream err_msg;
        err_msg << "Trying to get internal block dimension for \n"
                << "internal block " << b << ".\n"
                << "But there are only " << i_nblock_types
                << " internal dof types.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      return Internal_block_distribution_pt[b]->nrow();
    }
 
    unsigned internal_dof_block_dimension(const unsigned& i) const
    {
#ifdef PARANOID
      const unsigned i_n_dof_types = internal_ndof_types();
      if (i >= i_n_dof_types)
      {
        std::ostringstream err_msg;
        err_msg << "Trying to get internal dof block dimension for \n"
                << "internal dof block " << i << ".\n"
                << "But there are only " << i_n_dof_types
                << " internal dof types.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // I don't understand the difference between this function and
      // block_dimension(...) but I'm not going to mess with it... David
 
      if (is_master_block_preconditioner())
      {
        return Dof_dimension[i];
      }
      else
      {
        unsigned master_i = internal_master_dof_number(i);
        return Master_block_preconditioner_pt->internal_dof_block_dimension(
          master_i);
      }
    }
 
    unsigned master_nrow() const
    {
      if (is_master_block_preconditioner())
      {
        return Nrow;
      }
      else
      {
        return (this->Master_block_preconditioner_pt->master_nrow());
      }
    }
 
    unsigned internal_master_dof_number(const unsigned& b) const
    {
      if (is_master_block_preconditioner()) return b;
      else
        return Doftype_in_master_preconditioner_fine[b];
    }
 
    const LinearAlgebraDistribution* internal_preconditioner_matrix_distribution_pt()
      const
    {
      if (is_master_block_preconditioner())
        return Internal_preconditioner_matrix_distribution_pt;
      else
        return this->distribution_pt();
    }
 
    const LinearAlgebraDistribution* preconditioner_matrix_distribution_pt()
      const
    {
      return Preconditioner_matrix_distribution_pt;
    }
 
    MapMatrix<unsigned, CRDoubleMatrix*> Replacement_dof_block_pt;
 
    Vector<LinearAlgebraDistribution*> Block_distribution_pt;
 
    Vector<Vector<unsigned>> Block_to_dof_map_coarse;
 
    Vector<Vector<unsigned>> Block_to_dof_map_fine;
 
    Vector<Vector<unsigned>> Doftype_coarsen_map_coarse;
 
    Vector<Vector<unsigned>> Doftype_coarsen_map_fine;
 
    Vector<LinearAlgebraDistribution*> Internal_block_distribution_pt;
 
    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;
 
  private:
    bool Recursive_debug_flag;
 
    bool Debug_flag;
 
    std::map<Vector<unsigned>, LinearAlgebraDistribution*>
      Auxiliary_block_distribution_pt;
 
    unsigned Nrow;
 
    BlockPreconditioner<MATRIX>* Master_block_preconditioner_pt;
 
    Vector<unsigned> Doftype_in_master_preconditioner_fine;
 
    Vector<unsigned> Doftype_in_master_preconditioner_coarse;
 
    Vector<unsigned> Index_in_dof_block_dense;
 
    Vector<unsigned> Dof_number_dense;
 
#ifdef OOMPH_HAS_MPI
 
    // The following three vectors store data on the matrix rows/matrix
    // columns/dofs (the three are equivalent) that are not on this processor.
 
    Vector<unsigned> Global_index_sparse;
 
    Vector<unsigned> Index_in_dof_block_sparse;
 
    Vector<unsigned> Dof_number_sparse;
#endif
 
    Vector<unsigned> Dof_dimension;
 
    Vector<Vector<unsigned>> Global_index;
 
    Vector<Vector<unsigned>> Block_number_to_dof_number_lookup;
 
    Vector<unsigned> Dof_number_to_block_number_lookup;
 
    Vector<unsigned> Ndof_in_block;
 
 
#ifdef OOMPH_HAS_MPI
    DenseMatrix<int*> Rows_to_send_for_get_block;
 
    DenseMatrix<unsigned> Nrows_to_send_for_get_block;
 
    DenseMatrix<int*> Rows_to_recv_for_get_block;
 
    DenseMatrix<unsigned> Nrows_to_recv_for_get_block;
 
    Vector<int*> Rows_to_send_for_get_ordered;
 
    Vector<unsigned> Nrows_to_send_for_get_ordered;
 
    Vector<int*> Rows_to_recv_for_get_ordered;
 
    Vector<unsigned> Nrows_to_recv_for_get_ordered;
#endif
 
    LinearAlgebraDistribution* Internal_preconditioner_matrix_distribution_pt;
 
    LinearAlgebraDistribution* Preconditioner_matrix_distribution_pt;
 
    static bool Run_block_matrix_test;
 
    std::string Output_base_filename;
  };
 
 
} // namespace oomph
#endif