block_preconditioner.cc

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-2025 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 "block_preconditioner.h"
 
namespace oomph
{
  /// Static boolean to allow block_matrix_test(...) to be run.
  /// Defaults to false.
  template<typename MATRIX>
  bool BlockPreconditioner<MATRIX>::Run_block_matrix_test = false;
 
 
  //============================================================================
  /// Determine the size of the matrix blocks and setup the
  /// lookup schemes relating the global degrees of freedom with
  /// their "blocks" and their indices (row/column numbers) in those
  /// blocks.
  /// The distributions of the preconditioner and the blocks are
  /// automatically specified (and assumed to be uniform) at this
  /// stage.
  /// This method should be used if any block contains more than one
  /// type of DOF. The argument vector dof_to_block_map should be of length
  /// ndof. Each element should contain an integer indicating the block number
  /// corresponding to that type of DOF.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::block_setup(
    const Vector<unsigned>& dof_to_block_map_in)
  {
#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
 
    // Create a copy of the vector input so that we can modify it below
    Vector<unsigned> dof_to_block_map = dof_to_block_map_in;
 
    if (is_subsidiary_block_preconditioner())
    {
#ifdef PARANOID
      // Get the size of the Doftype_in_master_preconditioner_coarse.
      unsigned para_doftype_in_master_preconditioner_coarse_size =
        Doftype_in_master_preconditioner_coarse.size();
 
      // Check that the Doftype_in_master_preconditioner_coarse vector is not
      // empty. This must be set (via the function
      // turn_into_subsidiary_block_preconditioner) if this is a
      // subsidiary block preconditioner.
      if (para_doftype_in_master_preconditioner_coarse_size == 0)
      {
        std::ostringstream err_msg;
        err_msg << "The mapping from the dof types of the master "
                << "block preconditioner \n"
                << "to the subsidiary block preconditioner is empty.\n"
                << "Doftype_in_master_preconditioner_coarse.size() == 0 \n"
                << "has turn_into_subsidiary_block_preconditioner(...)\n"
                << "been called with the correct parameters?\n"
                << std::endl;
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
 
      // PARANOID checks for Doftype_coarsen_map_coarse
      // This is also set in the function
      // turn_into_subsidiary_block_preconditioner(...).
      //
      // The Doftype_coarsen_map_coarse vector must satisfy two conditions
      // for it to be valid.
      //
      // 1) The dof type numbers in the dof_coarsen_map vector must be
      //    unique. For example, it does not make sense to have the vector
      //    [[0,1][1,2]] because the first inner vector says
      //    "treat dof types 0 and 1 as dof type 0" and the second inner vector
      //    says "treat dof type 1 and 2 as dof type 1", but dof type 1 is
      //    already being treated as dof type 0.
      //
      // 2) Every SUBSIDIARY dof type must be mapped to a dof type in the
      //    Doftype_coarsen_map_coarse vector.
      //    For example, if there are 5 dof types (passed down from the master
      //    block preconditioner), and this block subsidiary block
      //    preconditioner only deals with 3 dof types, then all 5 dof types
      //    must be mapped to a dof type in the subsidiary preconditioner. For
      //    example if the dof_map is [1,2,3,4,5], then the subsidiary block
      //    preconditioner knows that 5 dof types have been passed down. But if
      //    it only works with three dof types, we MUST have three inner vectors
      //    in the doftype_coarsen_map vector (which corresponds to dof types 0,
      //    1 and 2), the union of the dof types in the three inner vectors must
      //    contain dof types 0, 1, 2, 3 and 4 exactly once. It cannot contain,
      //    say, 0, 1, 5, 7, 9, even though it passes the uniqueness check. We
      //    ensure this by two conditions:
      //
      //    2.1) The Doftype_coarsen_map_coarse vector must contain the same
      //         number of dof types as the dof_map vector.
      //         In other words, recall that Doftype_coarsen_map_coarse is a
      //         2D vector, this must contain the same number of vectors as
      //         there are elements in the dof_to_block_map_in vector.
      //
      //    2.2) The maximum element in the doftype_coarsen_map_coarse vector
      //         is the length of the dof_map vector minus 1.
 
      // A set is deal for checking the above three conditions, we shall insert
      // all the elements in the doftype_coarsen_map_coarse vector into this
      // set.
      std::set<unsigned> doftype_map_set;
 
      // Condition (1): Check for uniqueness by inserting all the values of
      // Doftype_coarsen_map_coarse into a set.
      unsigned para_doftype_coarsen_map_coarse_size =
        Doftype_coarsen_map_coarse.size();
 
      // Loop through the outer vector of Doftype_coarsen_map_coarse
      // then loop through the inner vectors and attempt to insert each
      // element of Doftype_coarsen_map_coarse into doftype_map_set.
      //
      // The inner for loop will throw an error if we cannot insert the
      // element, this means that it is already inserted and thus not unique.
      for (unsigned i = 0; i < para_doftype_coarsen_map_coarse_size; i++)
      {
        // Loop through the inner vector
        unsigned para_doftype_coarsen_map_coarse_i_size =
          Doftype_coarsen_map_coarse[i].size();
        for (unsigned j = 0; j < para_doftype_coarsen_map_coarse_i_size; j++)
        {
          // Attempt to insert all the values of the inner vector into a set.
          std::pair<std::set<unsigned>::iterator, bool> doftype_map_ret =
            doftype_map_set.insert(Doftype_coarsen_map_coarse[i][j]);
 
          if (!doftype_map_ret.second)
          {
            std::ostringstream err_msg;
            err_msg << "Error: the doftype number "
                    << Doftype_coarsen_map_coarse[i][j]
                    << " is already inserted." << std::endl;
            throw OomphLibError(
              err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
          }
        }
      }
 
      // Condition (2.1): Check that the doftype_map_set describes as many
      // values as doftype_in_master_preconditioner_coarse. I.e. if dof_map
      // contains 5 dof types, then the doftype_coarsen_map_coarse vector must
      // also contain 5 dof types.
      if (para_doftype_in_master_preconditioner_coarse_size !=
          doftype_map_set.size())
      {
        std::ostringstream err_msg;
        err_msg << "The size of doftype_in_master_preconditioner_coarse "
                << "must be the same as the total\n"
                << "number of values in the doftype_coarsen_map_coarse vector."
                << std::endl;
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Condition (2.2): Check that the maximum element in the
      // doftype_coarsen_map_coarse vector is the length of the
      // doftype_in_master_preconditioner_coarse minus 1.
      unsigned para_doftype_in_master_preconditioner_coarse_size_minus_one =
        para_doftype_in_master_preconditioner_coarse_size - 1;
      if (para_doftype_in_master_preconditioner_coarse_size_minus_one !=
          *doftype_map_set.rbegin())
      {
        std::ostringstream err_msg;
        err_msg << "The maximum dof type number in the "
                << "doftype_coarsen_map vector must be "
                << para_doftype_in_master_preconditioner_coarse_size_minus_one
                << std::endl;
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Set the mapping from the master preconditioner DOF types to the
      // subsidiary preconditioner DOF types.
      //
      // IMPORTANT: Since DOF types may be coarsened in the master block
      // preconditioner, this may no longer reflect the actual underlying dof
      // types. We must get the actual underlying dof types for the
      // block_setup(...) function to work properly so all the look up schemes
      // for this (subsidiary) block preconditioner is correct and works
      // properly, this is for backwards compatibility purposes and to make sure
      // Richard Muddle's still works at this (subsidiary) level, although it
      // may not be used.
      //
      // If we do not want to make it backwards compatible, we may as well
      // kill the block_setup(...) for subsidiary block preconditioners -
      // but other thing may break. Do it at your own risk (take time to
      // fully understand the whole block preconditioning framework code).
 
      // Create the corresponding Doftype_in_master_preconditioner_fine and
      // Doftype_coarsen_map_fine vectors.
 
      // First resize the vectors.
      Doftype_in_master_preconditioner_fine.resize(0);
      Doftype_coarsen_map_fine.resize(0);
 
      // The Doftype_in_master_preconditioner_fine vector is easy.  We know that
      // the Doftype_coarsen_map_fine in the master preconditioner must be
      // constructed already. So we simply loop through the values in
      // doftype_in_master_preconditioner_coarse, then get the most fine grain
      // dof types from the master preconditioner's Doftype_coarsen_map_fine
      // vector.
      //
      // For example, if the master preconditioner has the vector:
      // Doftype_coarsen_map_fine = [0,1,2,3][4,5,6,7][8,9,10,11][12,13][14,15]
      //
      // and passes the two vectors
      // doftype_in_master_preconditioner_coarse = [1,2,3]
      // doftype_coarsen_map_coarse = [[0][1,2]]
      //
      // Then we want
      // Doftype_in_master_preconditioner_fine = [4,5,6,7,8,9,10,11,12,13]
      //
      // We achieve this by looking up the corresponding fine dof types in the
      // masters' Doftype_coarsen_map_fine vector which corresponds to the
      // values in Doftype_in_master_preconditioner_coarse.
      //
      // That is, the values in Doftype_in_master_preconditioner_coarse gives us
      // the index of sub vector we want in the master's
      // Doftype_coarsen_map_fine vector.
 
#ifdef PARANOID
      // Check that the master block preconditioner's Doftype_coarsen_map_fine
      // is set up. Under the current implementation, this would always be set
      // up properly, but we check it just in case!
      if (master_block_preconditioner_pt()->doftype_coarsen_map_fine().size() ==
          0)
      {
        std::ostringstream err_msg;
        err_msg << "The master block preconditioner's "
                << "Doftype_coarsen_map_fine is not\n"
                << "set up properly.\n"
                << "\n"
                << "This vector is constructed in the function "
                << "block_setup(...).\n"
                << std::endl;
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      unsigned doftype_in_master_preconditioner_coarse_size =
        Doftype_in_master_preconditioner_coarse.size();
      for (unsigned i = 0; i < doftype_in_master_preconditioner_coarse_size;
           i++)
      {
        // The index of the sub vector we want.
        unsigned subvec_index = Doftype_in_master_preconditioner_coarse[i];
 
        // Get the corresponding most fine grain sub vector from the master
        // block preconditioner
        Vector<unsigned> tmp_master_dof_subvec =
          Master_block_preconditioner_pt->get_fine_grain_dof_types_in(
            subvec_index);
 
        Doftype_in_master_preconditioner_fine.insert(
          Doftype_in_master_preconditioner_fine.end(),
          tmp_master_dof_subvec.begin(),
          tmp_master_dof_subvec.end());
      }
 
      // The Doftype_coarsen_map_fine vector is a bit more tricky.
      // The Doftype_coarsen_map_coarse vector describes which coarse dof types
      // of THIS preconditioner are grouped together. We have to translate this
      // into the most fine grain dof types.
      //
      // For example, if
      // Doftype_coarsen_map_coarse = [[0][1,2]]
      // Doftype_in_master_preconditioner_coarse = [1,2,3]
      //
      // and the MASTER preconditioner has:
      // Doftype_coarsen_map_fine= [[0,1,2,3][4,5,6,7][8,9,10,11][12,13][14,15]]
      //
      // Then [[0][1,2]] tell us that the most fine grain DOF types 1 of the
      // master preconditioner most be grouped together, and the most fine
      // grained dof types 2 and 3 of the master preconditioner must be grouped
      // together.
      //
      // This gives the vector [[4,5,6,7] [8,9,10,11,12,13]], translating this
      // into the local DOF types of this preconditioner we have
      // Doftype_coarsen_map_fine = [[0,1,2,3][4,5,6,7,8,9]]. This corresponds
      // with the Doftype_in_master_preconditioner_fine vector we created above:
      // Doftype_in_master_preconditioner_fine = [4,5,6,7,8,9,10,11,12,13]
      //
      // Together, the master block preconditioner says to THIS subsidiary block
      // preconditioner "work on my DOF types [4,5,6,7,8,9,10,11,12,13], but
      // group your DOF type [0,1,2,3] together as DOF type 0 and [4,5,6,7,8,9]
      // together together as DOF type 1".
      //
      // Think of it like this: For each DOF type in Doftype_coarsen_map_coarse
      // we look at how many values this corresponds to in the master
      // preconditioner. In this case, Doftype_coarsen_map_coarse:
      //
      // 1 - corresponds to fine DOF types 0,1,2,3 in this preconditioner,
      // and 4,5,6,7 in the master preconditioner;
      //
      // 2 - corresponds to fine DOF types 4,5,6,7 in this preconditioner,
      // and 8,9,10,11 in the master preconditioner;
      //
      // 3 - corresponds to fine DOF types 8,9 in this preconditioner,
      // and 12,13 in the master preconditioner.
      //
      // Thus Doftype_coarsen_map_fine = [[0,1,2,3][4,5,6,7,8,9]]
      //
      ////////////////////////////////////////////////////////////////////////
      //
      // How to do this: First we create a 2D vector which has the corresponds
      // to the fine dof types in the master preconditioner but starting from
      // 0. For example, take the above example (repeated below):
      //   Passed to this prec by the master prec:
      //   Doftype_coarsen_map_coarse = [[0][1,2]]
      //   Doftype_in_master_preconditioner_coarse = [1,2,3]
      //
      // and the MASTER preconditioner has:
      // Doftype_coarsen_map_fine= [[0,1,2,3][4,5,6,7][8,9,10,11][12,13][14,15]]
      //
      // Step 1:
      // Then, the temp 2D vector we want to create is:
      // master_fine_doftype_translated = [[0 1 2 3], [4,5,6,7], [8,9]]
      // This comes from using Doftype_in_master_preconditioner_coarse
      // then get the number of fine dof types in the master.
      //
      // Step 2:
      // Then:
      //   Loop through the vector Doftype_coarsen_map_coarse,
      //     Loop over the inner vectors in Doftype_coarsen_map_coarse
      //       Each element in the inner vector corresponds to a vector in
      //       master_fine_doftype_translated. We push in the vectors of
      //       master_fine_doftype_translated intp Doftype_coarsen_map_fine
      //
 
      Vector<Vector<unsigned>> master_fine_doftype_translated;
      unsigned dof_type_index = 0;
      for (unsigned i = 0; i < doftype_in_master_preconditioner_coarse_size;
           i++)
      {
        // How many fine DOF types are in the master's
        // Doftype_in_master_preconditioner_coarse[i]?
        unsigned coarse_dof = Doftype_in_master_preconditioner_coarse[i];
 
        unsigned n_master_fine_doftypes =
          Master_block_preconditioner_pt->nfine_grain_dof_types_in(coarse_dof);
 
        Vector<unsigned> tmp_sub_vec;
        for (unsigned j = 0; j < n_master_fine_doftypes; j++)
        {
          tmp_sub_vec.push_back(dof_type_index);
          dof_type_index++;
        }
        master_fine_doftype_translated.push_back(tmp_sub_vec);
      }
 
 
      // master_fine_doftype_translated now contains vectors with values are
      // from 0, 1, 2, ..,
      //
      // Now read out the values of master_fine_doftype_translated and place
      // them in order according to Doftype_coarsen_map_coarse.
      unsigned doftype_coarsen_map_coarse_size =
        Doftype_coarsen_map_coarse.size();
      for (unsigned i = 0; i < doftype_coarsen_map_coarse_size; i++)
      {
        Vector<unsigned> tmp_vec;
        unsigned doftype_coarsen_map_coarse_i_size =
          Doftype_coarsen_map_coarse[i].size();
        for (unsigned j = 0; j < doftype_coarsen_map_coarse_i_size; j++)
        {
          unsigned subvec_i = Doftype_coarsen_map_coarse[i][j];
 
          tmp_vec.insert(tmp_vec.end(),
                         master_fine_doftype_translated[subvec_i].begin(),
                         master_fine_doftype_translated[subvec_i].end());
        }
 
        Doftype_coarsen_map_fine.push_back(tmp_vec);
      }
 
      // Get the number of block types (and DOF types) in this preconditioner
      // from the length of the dof_map vector.
      Internal_ndof_types = Doftype_in_master_preconditioner_fine.size();
 
      // Nblock_types is later updated in block_setup(...)
      Internal_nblock_types = Internal_ndof_types;
 
      // Compute number of rows in this (sub) preconditioner using data from
      // the master.
      Nrow = 0;
      for (unsigned b = 0; b < Internal_ndof_types; b++)
      {
        Nrow += this->internal_dof_block_dimension(b);
      }
 
#ifdef PARANOID
      if (Nrow == 0)
      {
        std::ostringstream error_message;
        error_message
          << "Nrow=0 in subsidiary preconditioner. This seems fishy and\n"
          << "suggests that block_setup() was not called for the \n"
          << "master block preconditioner yet.";
        throw OomphLibWarning(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
 
    // If this is a master block preconditioner, then set the
    // Doftype_coarsen_map_fine and Doftype_coarsen_map_coarse to the
    // identity. Recall that the Doftype_coarsen_map_fine maps the dof types
    // that this preconditioner requires with the most fine grain dof types (the
    // internal dof types) and the Doftype_coarsen_map_coarse maps the dof
    // types that this preconditioner requires with the dof types which this
    // preconditioner is given from a master preconditioner (these dof types may
    // or may not be coarsened). In the case of the master preconditioner, these
    // are the same (since dof types are not coarsened), furthermore the
    // identity mapping is provided to say that dof type 0 maps to dof type 0,
    // dof type 1 maps to dof type 1,
    // dof type 2 maps to dof type 2,
    // etc...
    //
    // If this is not a master block preconditioner, then the vectors
    // Doftype_coarsen_map_fine and Doftype_coarsen_map_coarse is handled
    // by the turn_into_subsidiary_block_preconditioner(...) function.
    if (is_master_block_preconditioner())
    {
      // How many dof types does this preconditioner work with?
      unsigned n_external_dof_types = dof_to_block_map.size();
 
      // Note: at the master level, the n_external_dof_types should be the same
      // as the internal_ndof_types(), since the dof_to_block_map MUST describe
      // the mapping between every dof type (not yet coarsened - so it is the
      // same number as the internal dof types) to the block types. But we
      // distinguish them for clarity. We also check that this is the case.
#ifdef PARANOID
      unsigned n_internal_dof_types = internal_ndof_types();
 
      if (n_internal_dof_types != n_external_dof_types)
      {
        std::ostringstream err_msg;
        err_msg
          << "The internal ndof types and the length of the dof_to_block_map\n"
          << "vector is not the same. Since this is the master block "
          << "preconditioner,\n"
          << "you must describe which block each DOF type belongs to,\n"
          << "no more, no less."
          << "internal_ndof_types = " << n_internal_dof_types << "\n"
          << "dof_to_block_map.size() = " << n_external_dof_types << "\n";
        throw OomphLibWarning(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Clear and reserve space.
      Doftype_coarsen_map_fine.clear();
      Doftype_coarsen_map_coarse.clear();
      Doftype_coarsen_map_fine.reserve(n_external_dof_types);
      Doftype_coarsen_map_coarse.reserve(n_external_dof_types);
 
      // Now push back the identity mapping.
      for (unsigned i = 0; i < n_external_dof_types; i++)
      {
        // Create a vector and push it in.
        Vector<unsigned> tmp_vec(1, i);
        Doftype_coarsen_map_fine.push_back(tmp_vec);
        Doftype_coarsen_map_coarse.push_back(tmp_vec);
      }
    }
    else
    // Else this is a subsidiary block preconditioner.
    {
      // Both the Doftype_coarsen_map_fine and Doftype_coarsen_map_coarse
      // vectors must be already be handled by the
      // turn_into_subsidiary_block_preconditioner(...) function. We check this.
#ifdef PARANOID
      if ((Doftype_coarsen_map_fine.size() == 0) ||
          (Doftype_coarsen_map_coarse.size() == 0))
      {
        std::ostringstream err_msg;
        err_msg << "Either the Doftype_coarsen_map_fine or the \n"
                << "Doftype_coarsen_map_coarse vectors is of size 0.\n"
                << "Did you remember to call the function "
                << "turn_into_subsidiary_block_preconditioner(...)?";
        throw OomphLibWarning(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
    }
 
 
    // Now we create the vector Block_to_dof_map_coarse.
    // Recall that the vector describe which dof types are in which block with
    // the relationship:
    //
    // Block_to_dof_map_coarse[block_number] = Vector[dof types];
    //
    // Note that this is not the internal (underlying) dof type.
    // Nor is this in relation to the parent block preconditioner's dof type.
    // The number of elements in it is the same as dof_to_block_map vector.
    //
    // Since the dof type coarsening feature is added later, we encapsulate this
    // bit of the code so it does not affect things below.
    {
      // Check that the dof_to_block map "makes sense" for the
      // Doftype_coarsen_map_coarse.
      // The Doftype_coarsen_map_coarse describes which doftypes should be
      // considered as a single doftype in THIS preconditioner.
      //
      // For example, if this preconditioner is the LSC block preconditioner
      // applied to a 3D problem, it expects 4 doftypes:
      // 3 velocity, [u, v, w] and 1 pressure [p],
      // giving us the dof type ordering
      // [u v w p].
      //
      // The LSC preconditioner groups the velocity and pressure doftypes
      // separately, thus the dof_to_block_map will be:
      // [0 0 0 1]
      //
      // Then the Doftype_coarsen_map_coarse MUST have length 4, to identify
      // which of the OTHER (possibly coarsened) dof types should be grouped
      // together to be considered as a single dof types of THIS preconditioner.
      //
      // For example, if the preconditioner above this one has the dof type
      // ordering:
      // 0  1  2  3  4  5  6  7  8  9
      // ub vb wb up vp wp ut vt wt p
      // Then we want to tell THIS preconditioner which dof types belongs to
      // u, v, w and p, by providing the optional argument
      // Doftype_coarsen_map_coarse to the
      // turn_into_subsidiary_block_preconditioner(...) function
      // [[0 3 6] [1 4 7] [2 5 8] [9]]
      //
      // So, it is important that the length of dof_to_block_map is the same as
      // the length of Doftype_coarsen_map_coarse. We check this.
      unsigned dof_to_block_map_size = dof_to_block_map.size();
 
#ifdef PARANOID
      if (dof_to_block_map_size != Doftype_coarsen_map_coarse.size())
      {
        std::ostringstream err_msg;
        err_msg
          << "The size of dof_to_block_map and Doftype_coarsen_map_coarse is "
             "not "
          << "the same.\n"
          << "dof_to_block_map.size() = " << dof_to_block_map_size << "\n"
          << "Doftype_coarsen_map_coarse.size() = "
          << Doftype_coarsen_map_coarse.size() << ".\n"
          << "One of the two list is incorrect, please look at the comments\n"
          << "in the source code for more details.";
        throw OomphLibWarning(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // Create the Block_to_dof_map_coarse from
      // the dof_to_block_map and Doftype_coarsen_map_coarse.
 
      // find the maximum block number
      unsigned max_block_number =
        *std::max_element(dof_to_block_map.begin(), dof_to_block_map.end());
 
      // Now we do the following:
      // Lets say the Doftype_coarsen_map_coarse is:
      // [0 3 6]
      // [1 4 7]
      // [2 5 8]
      // [9]
      //
      // (this is the same as the above example)
      //
      // and the dof_to_block_map is [0 0 0 1].
      //
      // Then we need to form the Block_to_dof_map_coarse:
      // [0 3 6 1 4 7 2 5 8]
      // [9]
 
      // Clear anything in the Block_to_dof_map_coarse
      Block_to_dof_map_coarse.clear();
 
      const unsigned tmp_nblock = max_block_number + 1;
 
      Block_to_dof_map_coarse.resize(tmp_nblock);
 
      for (unsigned i = 0; i < dof_to_block_map_size; i++)
      {
        Block_to_dof_map_coarse[dof_to_block_map[i]].push_back(i);
      }
 
      Block_to_dof_map_fine.clear();
      Block_to_dof_map_fine.resize(tmp_nblock);
      for (unsigned block_i = 0; block_i < tmp_nblock; block_i++)
      {
        // get the dof types in this block.
        const unsigned ndof_in_block = Block_to_dof_map_coarse[block_i].size();
        for (unsigned dof_i = 0; dof_i < ndof_in_block; dof_i++)
        {
          const unsigned coarsened_dof_i =
            Block_to_dof_map_coarse[block_i][dof_i];
 
          // Insert the most fine grain dofs which this dof_i corresponds to
          // into block_i
          Vector<unsigned> dof_i_dofs =
            Doftype_coarsen_map_fine[coarsened_dof_i];
 
          Block_to_dof_map_fine[block_i].insert(
            Block_to_dof_map_fine[block_i].end(),
            dof_i_dofs.begin(),
            dof_i_dofs.end());
        }
      }
 
      // Now set the dof_to_block_map to the identify.
      // NOTE: We are now using the internal n dof types. This is because the
      // dof type coarsening feature was built on top of the existing block
      // preconditioning framework which does not handle coarsening of dof
      // types. Hence, under the hood, it still works with the most fine grain
      // dof types and does not do any coarsening.
 
      // Locally cache the internal ndof types (using access function because
      // the Internal_ndof_type variable may not be set up yet if this is a
      // master preconditioner).
      unsigned tmp_internal_ndof_types = internal_ndof_types();
 
      dof_to_block_map.resize(tmp_internal_ndof_types, 0);
 
      for (unsigned i = 0; i < tmp_internal_ndof_types; i++)
      {
        dof_to_block_map[i] = i;
      }
    } // end of Block_to_dof_map_coarse encapsulation
 
#ifdef PARANOID
 
    // Check that the meshes are ok. This only needs to be done in the master
    // because subsidiary preconditioners don't do anything with the meshes
    // here.
    if (is_master_block_preconditioner())
    {
      // This is declared as local_nmesh because there are other variables
      // called nmesh floating about... but this will not exist if PARANOID is
      // switched on.
      unsigned local_nmesh = nmesh();
 
      // Check that some mesh pointers have been assigned.
      if (local_nmesh == 0)
      {
        std::ostringstream error_msg;
        error_msg << "Cannot setup blocks because no meshes have been set.";
        throw OomphLibError(
          error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Each mesh must contain elements with the same number of dof.
      // A stricter check is to ensure that the mesh contains only one type of
      // elements. This is relaxed in same cases.
      for (unsigned mesh_i = 0; mesh_i < local_nmesh; mesh_i++)
      {
        // The number of elements in the current mesh.
        unsigned n_element = mesh_pt(mesh_i)->nelement();
 
        // When the bulk mesh is distributed, there may not be any elements
        // in the surface mesh(es).
        if (n_element > 0)
        {
          // The string of the first element in the current mesh.
          std::string first_element_string =
            typeid(*(mesh_pt(mesh_i)->element_pt(0))).name();
 
          // If there are multiple element types in the current mesh,
          // we can at least make sure that they contain the same types of DOFs.
          if (bool(Allow_multiple_element_type_in_mesh[mesh_i]))
          {
            // The ndof types of the first element.
            unsigned first_element_ndof_type =
              mesh_pt(mesh_i)->element_pt(0)->ndof_types();
 
            // Loop through the meshes and compare the number of types of DOFs.
            for (unsigned el_i = 1; el_i < n_element; el_i++)
            {
              // The ndof type of the current element.
              unsigned current_element_ndof_type =
                mesh_pt(mesh_i)->element_pt(el_i)->ndof_types();
 
              // The string of the current element.
              std::string current_element_string =
                typeid(*(mesh_pt(mesh_i)->element_pt(el_i))).name();
 
              // Compare against the first element.
              if (current_element_ndof_type != first_element_ndof_type)
              {
                std::ostringstream error_message;
                error_message
                  << "Elements in the same mesh MUST have the same number of "
                     "types "
                  << "of DOFs.\n"
                  << "The element in mesh " << mesh_i << ", at position "
                  << el_i << " is: \n"
                  << current_element_string << ", \n"
                  << "with ndof types: " << current_element_ndof_type << ".\n"
                  << "The first element in the same mesh is: \n"
                  << first_element_string << ", \n"
                  << "with ndof types: " << first_element_ndof_type << ".\n";
                throw OomphLibError(error_message.str(),
                                    OOMPH_CURRENT_FUNCTION,
                                    OOMPH_EXCEPTION_LOCATION);
              }
            }
          }
          else
          // There should be only one type of elements in the current mesh.
          // Check that this is the case!
          {
            // Loop through the elements in the current mesh.
            for (unsigned el_i = 1; el_i < n_element; el_i++)
            {
              // The string of the current element.
              std::string current_element_string =
                typeid(*(mesh_pt(mesh_i)->element_pt(el_i))).name();
 
              // Compare against the first element.
              if (current_element_string.compare(first_element_string) != 0)
              {
                std::ostringstream error_message;
                error_message
                  << "By default, a mesh containing block preconditionable "
                  << "elements must contain only one type of element.\n"
                  << "The element in mesh " << mesh_i << ", at position "
                  << el_i << " is: \n"
                  << current_element_string << "\n"
                  << "The first element in the same mesh is: \n"
                  << first_element_string << "\n"
                  << "If this is correct, consider calling the set_mesh(...) "
                     "with\n"
                  << "the optional argument set true to allow multiple "
                     "element\n"
                  << "types in the same mesh.\n"
                  << "Note: A minimal requirement is that the elements in the "
                     "same\n"
                  << "mesh MUST have the same number of DOF types.";
                throw OomphLibError(error_message.str(),
                                    OOMPH_CURRENT_FUNCTION,
                                    OOMPH_EXCEPTION_LOCATION);
              }
            }
          }
        }
      }
    }
 
#endif
    // clear the memory
    this->clear_block_preconditioner_base();
 
    // get my_rank and nproc
#ifdef OOMPH_HAS_MPI
    unsigned my_rank = comm_pt()->my_rank();
    unsigned nproc = comm_pt()->nproc();
#endif
 
 
    /////////////////////////////////////////////////////////////////////////////
    // start of master block preconditioner only operations
    /////////////////////////////////////////////////////////////////////////////
#ifdef OOMPH_HAS_MPI
    unsigned* nreq_sparse = new unsigned[nproc]();
    unsigned* nreq_sparse_for_proc = new unsigned[nproc]();
    unsigned** index_in_dof_block_sparse_send = new unsigned*[nproc]();
    unsigned** dof_number_sparse_send = new unsigned*[nproc]();
    Vector<MPI_Request> send_requests_sparse;
    Vector<MPI_Request> recv_requests_sparse;
#endif
 
    // If this preconditioner is the master preconditioner then we need
    // to assemble the vectors : Dof_number
    //                           Index_in_dof_block
    if (is_master_block_preconditioner())
    {
      // Get the number of dof types in each mesh.
      Ndof_types_in_mesh.assign(nmesh(), 0);
      for (unsigned i = 0; i < nmesh(); i++)
      {
        Ndof_types_in_mesh[i] = mesh_pt(i)->ndof_types();
      }
      // Setup the distribution of this preconditioner, assumed to be the same
      // as the matrix if the matrix is distributable.
      if (dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt()))
      {
        this->build_distribution(
          dynamic_cast<DistributableLinearAlgebraObject*>(matrix_pt())
            ->distribution_pt());
      }
      else
      {
        LinearAlgebraDistribution dist(comm_pt(), matrix_pt()->nrow(), false);
        this->build_distribution(dist);
      }
      Nrow = matrix_pt()->nrow();
 
      // Boolean to indicate whether the matrix is actually distributed,
      // ie distributed and on more than one processor.
      bool matrix_distributed =
        (this->distribution_pt()->distributed() &&
         this->distribution_pt()->communicator_pt()->nproc() > 1);
 
 
      // Matrix must be a CR matrix.
      CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
 
      if (cr_matrix_pt == 0)
      {
        std::ostringstream error_message;
        error_message << "Block setup for distributed matrices only works "
                      << "for CRDoubleMatrices";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
 
 
      // Get distribution.
      unsigned first_row = this->distribution_pt()->first_row();
      unsigned nrow_local = this->distribution_pt()->nrow_local();
      unsigned last_row = first_row + nrow_local - 1;
 
#ifdef OOMPH_HAS_MPI
      // storage for the rows required by each processor in the dense
      // block lookup storage scheme
      // dense_required_rows(p,0) is the minimum global index required by proc p
      //                 ...(p,1) is the maximum global index required by proc p
      DenseMatrix<unsigned> dense_required_rows(nproc, 2);
      for (unsigned p = 0; p < nproc; p++)
      {
        dense_required_rows(p, 0) = this->distribution_pt()->first_row(p);
        dense_required_rows(p, 1) = this->distribution_pt()->first_row(p) +
                                    this->distribution_pt()->nrow_local(p) - 1;
      }
 
      // determine the global rows That are not in the range first_row to
      // first_row+nrow_local for which we should store the
      // Dof_index and Index_in_dof_block for
      // then send the lists to other processors
      std::set<unsigned> sparse_global_rows_for_block_lookup;
      if (matrix_distributed)
      {
        unsigned nnz = cr_matrix_pt->nnz();
        int* column_index = cr_matrix_pt->column_index();
        for (unsigned i = 0; i < nnz; i++)
        {
          unsigned ci = column_index[i];
          if (ci < first_row || ci > last_row)
          {
            sparse_global_rows_for_block_lookup.insert(ci);
          }
        }
      }
 
      int nsparse = sparse_global_rows_for_block_lookup.size();
 
      Global_index_sparse.resize(0);
      std::copy(sparse_global_rows_for_block_lookup.begin(),
                sparse_global_rows_for_block_lookup.end(),
                std::back_inserter(Global_index_sparse));
 
      Index_in_dof_block_sparse.resize(nsparse);
      Dof_number_sparse.resize(nsparse);
      sparse_global_rows_for_block_lookup.clear();
 
      Vector<MPI_Request> recv_requests_sparse_nreq;
      if (matrix_distributed)
      {
        MPI_Aint base_displacement_sparse;
        MPI_Get_address(nreq_sparse, &base_displacement_sparse);
 
        int zero = 0;
        for (unsigned p = 0; p < nproc; p++)
        {
          // determine the global eqn numbers required by this processor
          // that can be classified by processor p
          int begin = 0;
          for (int i = 0; i < nsparse; ++i)
          {
            if (Global_index_sparse[i] < dense_required_rows(p, 0))
            {
              ++begin;
            }
            else
            {
              if (Global_index_sparse[i] <= dense_required_rows(p, 1))
              {
                ++nreq_sparse[p];
              }
              else
              {
                break;
              }
            }
          }
 
          // if this processor has rows to be classified by proc p
          if (nreq_sparse[p] > 0)
          {
            // send the number of global eqn numbers
            MPI_Request req1;
            MPI_Isend(&nreq_sparse[p],
                      1,
                      MPI_UNSIGNED,
                      p,
                      31,
                      comm_pt()->mpi_comm(),
                      &req1);
            send_requests_sparse.push_back(req1);
 
            // send the global eqn numbers
            MPI_Request req2;
            MPI_Isend(&Global_index_sparse[begin],
                      nreq_sparse[p],
                      MPI_UNSIGNED,
                      p,
                      32,
                      comm_pt()->mpi_comm(),
                      &req2);
            send_requests_sparse.push_back(req2);
 
            // post the recvs for the data that will be returned
 
            // the datatypes, displacements, lengths for the two datatypes
            MPI_Datatype types[2];
            MPI_Aint displacements[2];
            int lengths[2];
 
            // index in dof block
            MPI_Type_contiguous(nreq_sparse[p], MPI_UNSIGNED, &types[0]);
            MPI_Type_commit(&types[0]);
            MPI_Get_address(&Index_in_dof_block_sparse[begin],
                            &displacements[0]);
            displacements[0] -= base_displacement_sparse;
            lengths[0] = 1;
 
            // dof number
            MPI_Type_contiguous(nreq_sparse[p], MPI_UNSIGNED, &types[1]);
            MPI_Type_commit(&types[1]);
            MPI_Get_address(&Dof_number_sparse[begin], &displacements[1]);
            displacements[1] -= base_displacement_sparse;
            lengths[1] = 1;
 
            // build the final type
            MPI_Datatype recv_type;
            MPI_Type_create_struct(
              2, lengths, displacements, types, &recv_type);
            MPI_Type_commit(&recv_type);
            MPI_Type_free(&types[0]);
            MPI_Type_free(&types[1]);
 
            // and recv
            MPI_Request req;
            MPI_Irecv(
              nreq_sparse, 1, recv_type, p, 33, comm_pt()->mpi_comm(), &req);
            recv_requests_sparse.push_back(req);
            MPI_Type_free(&recv_type);
          }
 
          // if no communication required, confirm this
          if (nreq_sparse[p] == 0)
          {
            MPI_Request req1;
            MPI_Isend(&zero, 1, MPI_INT, p, 31, comm_pt()->mpi_comm(), &req1);
            send_requests_sparse.push_back(req1);
          }
 
          //
          MPI_Request req;
          MPI_Irecv(&nreq_sparse_for_proc[p],
                    1,
                    MPI_UNSIGNED,
                    p,
                    31,
                    comm_pt()->mpi_comm(),
                    &req);
          recv_requests_sparse_nreq.push_back(req);
        }
      }
#endif
 
      // resize the storage
      Dof_number_dense.resize(nrow_local);
      Index_in_dof_block_dense.resize(nrow_local);
 
      // zero the number of dof types
      Internal_ndof_types = 0;
 
#ifdef PARANOID
      // Vector to keep track of previously assigned block numbers
      // to check consistency between multiple assignments.
      Vector<int> previously_assigned_block_number(nrow_local,
                                                   Data::Is_unclassified);
#endif
 
      // determine whether the problem is distribution
      bool problem_distributed = false;
 
      // the problem method distributed() is only accessible with MPI
#ifdef OOMPH_HAS_MPI
      problem_distributed = any_mesh_distributed();
#endif
 
      // if the problem is not distributed
      if (!problem_distributed)
      {
        // Offset for the block type in the overall system.
        // Different meshes contain different block-preconditionable
        // elements -- their blocks are added one after the other.
        unsigned dof_offset = 0;
 
        // Loop over all meshes.
        for (unsigned m = 0; m < nmesh(); m++)
        {
          // Number of elements in this mesh.
          unsigned n_element = mesh_pt(m)->nelement();
 
          // Find the number of block types that the elements in this mesh
          // are in charge of.
          unsigned ndof_in_element = ndof_types_in_mesh(m);
          Internal_ndof_types += ndof_in_element;
 
          for (unsigned e = 0; e < n_element; e++)
          {
            // List containing pairs of global equation number and
            // dof number for each global dof in an element.
            std::list<std::pair<unsigned long, unsigned>> dof_lookup_list;
 
            // Get list of blocks associated with the element's global unknowns.
            mesh_pt(m)->element_pt(e)->get_dof_numbers_for_unknowns(
              dof_lookup_list);
 
            // Loop over all entries in the list
            // and store the block number.
            typedef std::list<std::pair<unsigned long, unsigned>>::iterator IT;
            for (IT it = dof_lookup_list.begin(); it != dof_lookup_list.end();
                 it++)
            {
              unsigned long global_dof = it->first;
              if (global_dof >= unsigned(first_row) &&
                  global_dof <= unsigned(last_row))
              {
                unsigned dof_number = (it->second) + dof_offset;
                Dof_number_dense[global_dof - first_row] = dof_number;
 
#ifdef PARANOID
                // Check consistency of block numbers if assigned multiple times
                if (previously_assigned_block_number[global_dof - first_row] <
                    0)
                {
                  previously_assigned_block_number[global_dof - first_row] =
                    dof_number;
                }
#endif
              }
            }
          }
 
          // About to do the next mesh which contains block preconditionable
          // elements of a different type; all the dofs that these elements are
          // "in charge of" differ from the ones considered so far.
          // Bump up the block counter to make sure we're not overwriting
          // anything here
          dof_offset += ndof_in_element;
        }
 
#ifdef PARANOID
        // check that every global equation number has been allocated a dof type
        for (unsigned i = 0; i < nrow_local; i++)
        {
          if (previously_assigned_block_number[i] < 0)
          {
            std::ostringstream error_message;
            error_message << "Not all degrees of freedom have had DOF type "
                          << "numbers allocated. Dof number " << i
                          << " is unallocated.";
            throw OomphLibError(error_message.str(),
                                OOMPH_CURRENT_FUNCTION,
                                OOMPH_EXCEPTION_LOCATION);
          }
        }
#endif
      }
      // else the problem is distributed
      else
      {
#ifdef OOMPH_HAS_MPI
        // Offset for the block type in the overall system.
        // Different meshes contain different block-preconditionable
        // elements -- their blocks are added one after the other...
        unsigned dof_offset = 0;
 
        // the set of global degrees of freedom and their corresponding dof
        // number on this processor
        std::map<unsigned long, unsigned> my_dof_map;
 
        // Loop over all meshes
        for (unsigned m = 0; m < nmesh(); m++)
        {
          // Number of elements in this mesh
          unsigned n_element = this->mesh_pt(m)->nelement();
 
          // Find the number of block types that the elements in this mesh
          // are in charge of
          unsigned ndof_in_element = ndof_types_in_mesh(m);
          Internal_ndof_types += ndof_in_element;
 
          // Loop over all elements
          for (unsigned e = 0; e < n_element; e++)
          {
            // if the element is not a halo element
            if (!this->mesh_pt(m)->element_pt(e)->is_halo())
            {
              // List containing pairs of global equation number and
              // dof number for each global dof in an element
              std::list<std::pair<unsigned long, unsigned>> dof_lookup_list;
 
              // Get list of blocks associated with the element's global
              // unknowns
              this->mesh_pt(m)->element_pt(e)->get_dof_numbers_for_unknowns(
                dof_lookup_list);
 
              // update the block numbers and put it in the map.
              typedef std::list<std::pair<unsigned long, unsigned>>::iterator
                IT;
              for (IT it = dof_lookup_list.begin(); it != dof_lookup_list.end();
                   it++)
              {
                it->second = (it->second) + dof_offset;
                my_dof_map[it->first] = it->second;
              }
            }
          }
 
          // About to do the next mesh which contains block preconditionable
          // elements of a different type; all the dofs that these elements are
          // "in charge of" differ from the ones considered so far.
          // Bump up the block counter to make sure we're not overwriting
          // anything here
          dof_offset += ndof_in_element;
        }
 
        // next copy the map of my dofs to two vectors to send
        unsigned my_ndof = my_dof_map.size();
        unsigned long* my_global_dofs = new unsigned long[my_ndof];
        unsigned* my_dof_numbers = new unsigned[my_ndof];
        typedef std::map<unsigned long, unsigned>::iterator IT;
        unsigned pt = 0;
        for (IT it = my_dof_map.begin(); it != my_dof_map.end(); it++)
        {
          my_global_dofs[pt] = it->first;
          my_dof_numbers[pt] = it->second;
          pt++;
        }
 
        // and then clear the map
        my_dof_map.clear();
 
        // count up how many DOFs need to be sent to each processor
        int* first_dof_to_send = new int[nproc];
        int* ndof_to_send = new int[nproc];
        unsigned ptr = 0;
        for (unsigned p = 0; p < nproc; p++)
        {
          first_dof_to_send[p] = 0;
          ndof_to_send[p] = 0;
          while (ptr < my_ndof &&
                 my_global_dofs[ptr] < dense_required_rows(p, 0))
          {
            ptr++;
          }
          first_dof_to_send[p] = ptr;
          while (ptr < my_ndof &&
                 my_global_dofs[ptr] <= dense_required_rows(p, 1))
          {
            ndof_to_send[p]++;
            ptr++;
          }
        }
 
        // next communicate to each processor how many dofs it expects to recv
        int* ndof_to_recv = new int[nproc];
        MPI_Alltoall(ndof_to_send,
                     1,
                     MPI_INT,
                     ndof_to_recv,
                     1,
                     MPI_INT,
                     comm_pt()->mpi_comm());
 
        // the base displacements for the sends
        MPI_Aint base_displacement;
        MPI_Get_address(my_global_dofs, &base_displacement);
 
#ifdef PARANOID
        // storage for paranoid check to ensure that every row as been
        // imported
        std::vector<bool> dof_recv(nrow_local, false);
#endif
 
        // next send and recv
        Vector<MPI_Request> send_requests;
        Vector<MPI_Request> recv_requests;
        Vector<unsigned long*> global_dofs_recv(nproc, 0);
        Vector<unsigned*> dof_numbers_recv(nproc, 0);
        Vector<unsigned> proc;
        for (unsigned p = 0; p < nproc; p++)
        {
          if (p != my_rank)
          {
            // send
            if (ndof_to_send[p] > 0)
            {
              // the datatypes, displacements, lengths for the two datatypes
              MPI_Datatype types[2];
              MPI_Aint displacements[2];
              int lengths[2];
 
              // my global dofs
              MPI_Type_contiguous(
                ndof_to_send[p], MPI_UNSIGNED_LONG, &types[0]);
              MPI_Type_commit(&types[0]);
              MPI_Get_address(my_global_dofs + first_dof_to_send[p],
                              &displacements[0]);
              displacements[0] -= base_displacement;
              lengths[0] = 1;
 
              // my dof numbers
              MPI_Type_contiguous(ndof_to_send[p], MPI_UNSIGNED, &types[1]);
              MPI_Type_commit(&types[1]);
              MPI_Get_address(my_dof_numbers + first_dof_to_send[p],
                              &displacements[1]);
              displacements[1] -= base_displacement;
              lengths[1] = 1;
 
              // build the final type
              MPI_Datatype send_type;
              MPI_Type_create_struct(
                2, lengths, displacements, types, &send_type);
              MPI_Type_commit(&send_type);
              MPI_Type_free(&types[0]);
              MPI_Type_free(&types[1]);
 
              // and send
              MPI_Request req;
              MPI_Isend(my_global_dofs,
                        1,
                        send_type,
                        p,
                        2,
                        comm_pt()->mpi_comm(),
                        &req);
              send_requests.push_back(req);
              MPI_Type_free(&send_type);
            }
 
            // and recv
            if (ndof_to_recv[p] > 0)
            {
              // resize the storage
              global_dofs_recv[p] = new unsigned long[ndof_to_recv[p]];
              dof_numbers_recv[p] = new unsigned[ndof_to_recv[p]];
              proc.push_back(p);
 
              // the datatypes, displacements, lengths for the two datatypes
              MPI_Datatype types[2];
              MPI_Aint displacements[2];
              int lengths[2];
 
              // my global dofs
              MPI_Type_contiguous(
                ndof_to_recv[p], MPI_UNSIGNED_LONG, &types[0]);
              MPI_Type_commit(&types[0]);
              MPI_Get_address(global_dofs_recv[p], &displacements[0]);
              displacements[0] -= base_displacement;
              lengths[0] = 1;
 
              // my dof numbers
              MPI_Type_contiguous(ndof_to_recv[p], MPI_UNSIGNED, &types[1]);
              MPI_Type_commit(&types[1]);
              MPI_Get_address(dof_numbers_recv[p], &displacements[1]);
              displacements[1] -= base_displacement;
              lengths[1] = 1;
 
              // build the final type
              MPI_Datatype recv_type;
              MPI_Type_create_struct(
                2, lengths, displacements, types, &recv_type);
              MPI_Type_commit(&recv_type);
              MPI_Type_free(&types[0]);
              MPI_Type_free(&types[1]);
 
              // and recv
              MPI_Request req;
              MPI_Irecv(my_global_dofs,
                        1,
                        recv_type,
                        p,
                        2,
                        comm_pt()->mpi_comm(),
                        &req);
              recv_requests.push_back(req);
              MPI_Type_free(&recv_type);
            }
          }
          // send to self
          else
          {
            unsigned u = first_dof_to_send[p] + ndof_to_recv[p];
            for (unsigned i = first_dof_to_send[p]; i < u; i++)
            {
#ifdef PARANOID
              // indicate that this dof has ben recv
              dof_recv[my_global_dofs[i] - first_row] = true;
#endif
              Dof_number_dense[my_global_dofs[i] - first_row] =
                my_dof_numbers[i];
            }
          }
        }
 
        // recv and store the data
        unsigned c_recv = recv_requests.size();
        while (c_recv > 0)
        {
          // wait for any communication to finish
          int req_number;
          MPI_Waitany(
            c_recv, &recv_requests[0], &req_number, MPI_STATUS_IGNORE);
          recv_requests.erase(recv_requests.begin() + req_number);
          c_recv--;
 
          // determine the source processor
          unsigned p = proc[req_number];
          proc.erase(proc.begin() + req_number);
 
          // import the data
          for (int i = 0; i < ndof_to_recv[p]; i++)
          {
#ifdef PARANOID
            // indicate that this dof has ben recv
            dof_recv[global_dofs_recv[p][i] - first_row] = true;
#endif
            Dof_number_dense[global_dofs_recv[p][i] - first_row] =
              dof_numbers_recv[p][i];
          }
 
          // delete the data
          delete[] global_dofs_recv[p];
          delete[] dof_numbers_recv[p];
        }
 
        // finally wait for the send requests to complete as we are leaving
        // an MPI block of code
        unsigned csr = send_requests.size();
        if (csr)
        {
          MPI_Waitall(csr, &send_requests[0], MPI_STATUS_IGNORE);
        }
 
        // clean up
        delete[] ndof_to_send;
        delete[] first_dof_to_send;
        delete[] ndof_to_recv;
        delete[] my_global_dofs;
        delete[] my_dof_numbers;
#ifdef PARANOID
        unsigned all_recv = true;
        for (unsigned i = 0; i < nrow_local; i++)
        {
          if (!dof_recv[i])
          {
            all_recv = false;
          }
        }
        if (!all_recv)
        {
          std::ostringstream error_message;
          error_message << "Not all the DOF numbers required were received";
          throw OomphLibError(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
        }
#endif
#else
        std::ostringstream error_message;
        error_message
          << "The problem appears to be distributed, MPI is required";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
#endif
      }
#ifdef OOMPH_HAS_MPI
      Vector<unsigned*> sparse_rows_for_proc(nproc, 0);
      Vector<MPI_Request> sparse_rows_for_proc_requests;
      if (matrix_distributed)
      {
        // wait for number of sparse rows each processor requires
        // post recvs for that data
        if (recv_requests_sparse_nreq.size() > 0)
        {
          MPI_Waitall(recv_requests_sparse_nreq.size(),
                      &recv_requests_sparse_nreq[0],
                      MPI_STATUS_IGNORE);
        }
        for (unsigned p = 0; p < nproc; ++p)
        {
          if (nreq_sparse_for_proc[p] > 0)
          {
            MPI_Request req;
            sparse_rows_for_proc[p] = new unsigned[nreq_sparse_for_proc[p]];
            MPI_Irecv(sparse_rows_for_proc[p],
                      nreq_sparse_for_proc[p],
                      MPI_UNSIGNED,
                      p,
                      32,
                      comm_pt()->mpi_comm(),
                      &req);
            sparse_rows_for_proc_requests.push_back(req);
          }
        }
      }
#endif
 
 
      // for every global degree of freedom required by this processor we now
      // have the corresponding dof number
 
      // clear the Ndof_in_dof_block storage
      Dof_dimension.assign(Internal_ndof_types, 0);
 
      // first consider a non distributed matrix
      if (!matrix_distributed)
      {
        // set the Index_in_dof_block
        unsigned nrow = this->distribution_pt()->nrow();
        Index_in_dof_block_dense.resize(nrow);
        Index_in_dof_block_dense.initialise(0);
        for (unsigned i = 0; i < nrow; i++)
        {
          Index_in_dof_block_dense[i] = Dof_dimension[Dof_number_dense[i]];
          Dof_dimension[Dof_number_dense[i]]++;
        }
      }
 
      // next a distributed matrix
      else
      {
#ifdef OOMPH_HAS_MPI
 
 
        // first compute how many instances of each dof are on this
        // processor
        unsigned* my_nrows_in_dof_block = new unsigned[Internal_ndof_types];
        for (unsigned i = 0; i < Internal_ndof_types; i++)
        {
          my_nrows_in_dof_block[i] = 0;
        }
        for (unsigned i = 0; i < nrow_local; i++)
        {
          my_nrows_in_dof_block[Dof_number_dense[i]]++;
        }
 
        // next share the data
        unsigned* nrow_in_dof_block_recv =
          new unsigned[Internal_ndof_types * nproc];
        MPI_Allgather(my_nrows_in_dof_block,
                      Internal_ndof_types,
                      MPI_UNSIGNED,
                      nrow_in_dof_block_recv,
                      Internal_ndof_types,
                      MPI_UNSIGNED,
                      comm_pt()->mpi_comm());
        delete[] my_nrows_in_dof_block;
 
        // compute my first dof index and Nrows_in_dof_block
        Vector<unsigned> my_first_dof_index(Internal_ndof_types, 0);
        for (unsigned i = 0; i < Internal_ndof_types; i++)
        {
          for (unsigned p = 0; p < my_rank; p++)
          {
            my_first_dof_index[i] +=
              nrow_in_dof_block_recv[p * Internal_ndof_types + i];
          }
          Dof_dimension[i] = my_first_dof_index[i];
          for (unsigned p = my_rank; p < nproc; p++)
          {
            Dof_dimension[i] +=
              nrow_in_dof_block_recv[p * Internal_ndof_types + i];
          }
        }
        delete[] nrow_in_dof_block_recv;
 
        // next compute Index in dof block
        Index_in_dof_block_dense.resize(nrow_local);
        Index_in_dof_block_dense.initialise(0);
        Vector<unsigned> dof_counter(Internal_ndof_types, 0);
        for (unsigned i = 0; i < nrow_local; i++)
        {
          Index_in_dof_block_dense[i] =
            my_first_dof_index[Dof_number_dense[i]] +
            dof_counter[Dof_number_dense[i]];
          dof_counter[Dof_number_dense[i]]++;
        }
 
        // the base displacements for the sends
        if (sparse_rows_for_proc_requests.size() > 0)
        {
          MPI_Waitall(sparse_rows_for_proc_requests.size(),
                      &sparse_rows_for_proc_requests[0],
                      MPI_STATUS_IGNORE);
        }
        MPI_Aint base_displacement;
        MPI_Get_address(dof_number_sparse_send, &base_displacement);
        unsigned first_row = this->distribution_pt()->first_row();
        for (unsigned p = 0; p < nproc; ++p)
        {
          if (nreq_sparse_for_proc[p] > 0)
          {
            // construct the data
            index_in_dof_block_sparse_send[p] =
              new unsigned[nreq_sparse_for_proc[p]];
            dof_number_sparse_send[p] = new unsigned[nreq_sparse_for_proc[p]];
            for (unsigned i = 0; i < nreq_sparse_for_proc[p]; ++i)
            {
              unsigned r = sparse_rows_for_proc[p][i];
              r -= first_row;
              index_in_dof_block_sparse_send[p][i] =
                Index_in_dof_block_dense[r];
              dof_number_sparse_send[p][i] = Dof_number_dense[r];
            }
            delete[] sparse_rows_for_proc[p];
 
            // send the data
            // the datatypes, displacements, lengths for the two datatypes
            MPI_Datatype types[2];
            MPI_Aint displacements[2];
            int lengths[2];
 
            // index in dof block
            MPI_Type_contiguous(
              nreq_sparse_for_proc[p], MPI_UNSIGNED, &types[0]);
            MPI_Type_commit(&types[0]);
            MPI_Get_address(index_in_dof_block_sparse_send[p],
                            &displacements[0]);
            displacements[0] -= base_displacement;
            lengths[0] = 1;
 
            // dof number
            MPI_Type_contiguous(
              nreq_sparse_for_proc[p], MPI_UNSIGNED, &types[1]);
            MPI_Type_commit(&types[1]);
            MPI_Get_address(dof_number_sparse_send[p], &displacements[1]);
            displacements[1] -= base_displacement;
            lengths[1] = 1;
 
            // build the final type
            MPI_Datatype send_type;
            MPI_Type_create_struct(
              2, lengths, displacements, types, &send_type);
            MPI_Type_commit(&send_type);
            MPI_Type_free(&types[0]);
            MPI_Type_free(&types[1]);
 
            // and recv
            MPI_Request req;
            MPI_Isend(dof_number_sparse_send,
                      1,
                      send_type,
                      p,
                      33,
                      comm_pt()->mpi_comm(),
                      &req);
            send_requests_sparse.push_back(req);
            MPI_Type_free(&send_type);
          }
          else
          {
            index_in_dof_block_sparse_send[p] = 0;
            dof_number_sparse_send[p] = 0;
          }
        }
#endif
      }
    }
 
    /////////////////////////////////////////////////////////////////////////////
    // end of master block preconditioner only operations
    /////////////////////////////////////////////////////////////////////////////
 
    // compute the number of rows in each block
 
#ifdef PARANOID
    // check the vector is the correct length
    if (dof_to_block_map.size() != Internal_ndof_types)
    {
      std::ostringstream error_message;
      error_message << "The dof_to_block_map vector (size="
                    << dof_to_block_map.size()
                    << ") must be of size Internal_ndof_types="
                    << Internal_ndof_types;
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // find the maximum block number RAYAY use std::max_element
    unsigned max_block_number = 0;
    for (unsigned i = 0; i < Internal_ndof_types; i++)
    {
      if (dof_to_block_map[i] > max_block_number)
      {
        max_block_number = dof_to_block_map[i];
      }
    }
 
    // resize the storage the the block to dof map
    Block_number_to_dof_number_lookup.clear();
    Block_number_to_dof_number_lookup.resize(max_block_number + 1);
    Ndof_in_block.clear();
    Ndof_in_block.resize(max_block_number + 1);
 
    // resize storage
    Dof_number_to_block_number_lookup.resize(Internal_ndof_types);
 
    // build the storage for the two maps (block to dof) and (dof to block)
    for (unsigned i = 0; i < Internal_ndof_types; i++)
    {
      Dof_number_to_block_number_lookup[i] = dof_to_block_map[i];
      Block_number_to_dof_number_lookup[dof_to_block_map[i]].push_back(i);
      Ndof_in_block[dof_to_block_map[i]]++;
    }
 
#ifdef PARANOID
    // paranoid check that every block number has at least one DOF associated
    // with it
    for (unsigned i = 0; i < max_block_number + 1; i++)
    {
      if (Block_number_to_dof_number_lookup[i].size() == 0)
      {
        std::ostringstream error_message;
        error_message << "block number " << i
                      << " does not have any DOFs associated with it";
        throw OomphLibWarning(error_message.str(),
                              OOMPH_CURRENT_FUNCTION,
                              OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // Update the number of blocks types.
    Internal_nblock_types = max_block_number + 1;
 
    // Distributed or not, depends on if we have more than one processor.
    bool distributed = this->master_distribution_pt()->distributed();
 
    // Create the new block distributions.
    Internal_block_distribution_pt.resize(Internal_nblock_types);
    for (unsigned i = 0; i < Internal_nblock_types; i++)
    {
      unsigned block_dim = 0;
      for (unsigned j = 0; j < Ndof_in_block[i]; j++)
      {
        block_dim +=
          internal_dof_block_dimension(Block_number_to_dof_number_lookup[i][j]);
      }
      Internal_block_distribution_pt[i] =
        new LinearAlgebraDistribution(comm_pt(), block_dim, distributed);
    }
 
    // Work out the distribution of the dof-level blocks.
    // Since several dof types may be coarsened into a single dof type.
    // We get the dof-level block distributions from the parent preconditioner.
 
    // How many dof types are there?
    if (is_subsidiary_block_preconditioner())
    {
      // Delete any pre-existing distributions.
      const unsigned dof_block_distribution_size =
        Dof_block_distribution_pt.size();
      for (unsigned dof_i = 0; dof_i < dof_block_distribution_size; dof_i++)
      {
        delete Dof_block_distribution_pt[dof_i];
      }
      const unsigned ndofs = this->ndof_types();
      Dof_block_distribution_pt.resize(ndofs, 0);
 
      // For each dof type, work out how many parent preconditioner dof types
      // are in it.
      for (unsigned dof_i = 0; dof_i < ndofs; dof_i++)
      {
        // For each external dof, we get the dofs coarsened into it (from the
        // parent preconditioner level, not the most fine grain level).
        const unsigned ncoarsened_dofs_in_dof_i =
          Doftype_coarsen_map_coarse[dof_i].size();
        Vector<LinearAlgebraDistribution*> tmp_dist_pt(ncoarsened_dofs_in_dof_i,
                                                       0);
        for (unsigned parent_dof_i = 0; parent_dof_i < ncoarsened_dofs_in_dof_i;
             parent_dof_i++)
        {
          tmp_dist_pt[parent_dof_i] =
            master_block_preconditioner_pt()->dof_block_distribution_pt(
              Doftype_in_master_preconditioner_coarse
                [Doftype_coarsen_map_coarse[dof_i][parent_dof_i]]);
        }
 
        Dof_block_distribution_pt[dof_i] = new LinearAlgebraDistribution;
 
 
        LinearAlgebraDistributionHelpers::concatenate(
          tmp_dist_pt, *Dof_block_distribution_pt[dof_i]);
      }
    }
 
    // Create Block_distribution_pt
    {
      // Delete any existing distributions in Block_distribution_pt.
      // (This should already be deleted in clear_block_preconditioner_base(...)
      // but we are just being extra safe!).
      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];
      }
 
      Block_distribution_pt.clear();
 
      // Work out the distributions of the concatenated blocks.
      unsigned super_block_size = Block_to_dof_map_coarse.size();
      Block_distribution_pt.resize(super_block_size, 0);
      for (unsigned super_block_i = 0; super_block_i < super_block_size;
           super_block_i++)
      {
        unsigned sub_block_size = Block_to_dof_map_coarse[super_block_i].size();
        Vector<LinearAlgebraDistribution*> tmp_dist_pt(sub_block_size, 0);
 
        for (unsigned sub_block_i = 0; sub_block_i < sub_block_size;
             sub_block_i++)
        {
          tmp_dist_pt[sub_block_i] = dof_block_distribution_pt(
            Block_to_dof_map_coarse[super_block_i][sub_block_i]);
        }
 
        Block_distribution_pt[super_block_i] = new LinearAlgebraDistribution;
 
        LinearAlgebraDistributionHelpers::concatenate(
          tmp_dist_pt, *Block_distribution_pt[super_block_i]);
      }
 
    } // Creating Block_distribution_pt.
 
 
    // Create the distribution of the preconditioner matrix,
    // if this preconditioner is a subsidiary preconditioner then it stored
    // at Distribution_pt;
    // if this preconditioner is a master preconditioner then it is stored
    // at Internal_preconditioner_matrix_distribution_pt.
    LinearAlgebraDistribution dist;
    LinearAlgebraDistributionHelpers::concatenate(
      Internal_block_distribution_pt, dist);
 
    // Build the distribution.
    if (is_subsidiary_block_preconditioner())
    {
      this->build_distribution(dist);
    }
    else
    {
      Internal_preconditioner_matrix_distribution_pt =
        new LinearAlgebraDistribution(dist);
    }
 
    Preconditioner_matrix_distribution_pt = new LinearAlgebraDistribution;
    LinearAlgebraDistributionHelpers::concatenate(
      Block_distribution_pt, *Preconditioner_matrix_distribution_pt);
 
    // Clear all distributions in Auxiliary_block_distribution_pt, except for
    // the one which corresponds to the preconditioner matrix distribution. This
    // is already deleted by clear_block_preconditioner_base(...)
 
    // Create the key which corresponds to
    // preconditioner_matrix_distribution_pt.
    {
      const unsigned nblocks = Block_distribution_pt.size();
      Vector<unsigned> preconditioner_matrix_key(nblocks, 0);
      for (unsigned i = 0; i < nblocks; i++)
      {
        preconditioner_matrix_key[i] = i;
      }
 
      // Now iterate through Auxiliary_block_distribution_pt and delete
      // everything except for the value which corresponds to
      // preconditioner_matrix_key.
      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;
        }
      }
 
      // Clear it just to be safe!
      Auxiliary_block_distribution_pt.clear();
 
      // Insert the preconditioner matrix distribution.
      insert_auxiliary_block_distribution(
        preconditioner_matrix_key, Preconditioner_matrix_distribution_pt);
    } // End of Auxiliary_block_distribution_pt encapsulation.
 
    // Clearing up after comm to assemble sparse lookup schemes.
#ifdef OOMPH_HAS_MPI
    if (send_requests_sparse.size() > 0)
    {
      MPI_Waitall(send_requests_sparse.size(),
                  &send_requests_sparse[0],
                  MPI_STATUS_IGNORE);
    }
    if (recv_requests_sparse.size() > 0)
    {
      MPI_Waitall(recv_requests_sparse.size(),
                  &recv_requests_sparse[0],
                  MPI_STATUS_IGNORE);
    }
    for (unsigned p = 0; p < nproc; p++)
    {
      delete[] index_in_dof_block_sparse_send[p];
      delete[] dof_number_sparse_send[p];
    }
    delete[] index_in_dof_block_sparse_send;
    delete[] dof_number_sparse_send;
    delete[] nreq_sparse;
    delete[] nreq_sparse_for_proc;
#endif
 
    // Next we assemble the lookup schemes for the rows
    // if the matrix is not distributed then we assemble Global_index
    // if the matrix is distributed then Rows_to_send_..., Rows_to_recv_... etc.
    if (!distributed)
    {
      // Resize the storage.
      Global_index.resize(Internal_nblock_types);
      for (unsigned b = 0; b < Internal_nblock_types; b++)
      {
        Global_index[b].resize(Internal_block_distribution_pt[b]->nrow());
      }
 
      // Compute:
      unsigned nrow = this->master_nrow();
      for (unsigned i = 0; i < nrow; i++)
      {
        // the dof type number;
        int dof_number = this->internal_dof_number(i);
        if (dof_number >= 0)
        {
          // the block number;
          unsigned block_number = Dof_number_to_block_number_lookup[dof_number];
 
          // the index in the block.
          unsigned index_in_block = 0;
          unsigned ptr = 0;
          while (int(Block_number_to_dof_number_lookup[block_number][ptr]) !=
                 dof_number)
          {
            index_in_block += internal_dof_block_dimension(
              Block_number_to_dof_number_lookup[block_number][ptr]);
            ptr++;
          }
          index_in_block += internal_index_in_dof(i);
          Global_index[block_number][index_in_block] = i;
        }
      }
    }
    // otherwise the matrix is distributed
    else
    {
#ifdef OOMPH_HAS_MPI
 
      // the pointer to the master distribution
      const LinearAlgebraDistribution* master_distribution_pt =
        this->master_distribution_pt();
 
      // resize the nrows... storage
      Nrows_to_send_for_get_block.resize(Internal_nblock_types, nproc);
      Nrows_to_send_for_get_block.initialise(0);
      Nrows_to_send_for_get_ordered.resize(nproc);
      Nrows_to_send_for_get_ordered.initialise(0);
 
      // loop over my rows
      unsigned nrow_local = master_distribution_pt->nrow_local();
      unsigned first_row = master_distribution_pt->first_row();
      for (unsigned i = 0; i < nrow_local; i++)
      {
        // the block number
        int b = this->internal_block_number(first_row + i);
 
        // check that the DOF i is associated with this preconditioner
        if (b >= 0)
        {
          // the block index
          unsigned j = this->internal_index_in_block(first_row + i);
 
          // the processor this row will be sent to
          unsigned block_p = 0;
          while (!(Internal_block_distribution_pt[b]->first_row(block_p) <= j &&
                   (Internal_block_distribution_pt[b]->first_row(block_p) +
                      Internal_block_distribution_pt[b]->nrow_local(block_p) >
                    j)))
          {
            block_p++;
          }
 
          // and increment the counter
          Nrows_to_send_for_get_block(b, block_p)++;
          Nrows_to_send_for_get_ordered[block_p]++;
        }
      }
 
      // resize the storage for Nrows_to_recv
      Nrows_to_recv_for_get_block.resize(Internal_nblock_types, nproc);
      Nrows_to_recv_for_get_block.initialise(0);
      Nrows_to_recv_for_get_ordered.resize(nproc);
      Nrows_to_recv_for_get_ordered.initialise(0);
 
      // next we send the number of rows that will be sent by this processor
      Vector<unsigned*> nrows_to_send(nproc, 0);
      Vector<unsigned*> nrows_to_recv(nproc, 0);
      Vector<MPI_Request> send_requests_nrow;
      Vector<MPI_Request> recv_requests_nrow;
      Vector<unsigned> proc;
      for (unsigned p = 0; p < nproc; p++)
      {
        if (p != my_rank)
        {
          // send
          proc.push_back(p);
          nrows_to_send[p] = new unsigned[Internal_nblock_types];
          for (unsigned b = 0; b < Internal_nblock_types; b++)
          {
            nrows_to_send[p][b] = Nrows_to_send_for_get_block(b, p);
          }
          MPI_Request s_req;
          MPI_Isend(nrows_to_send[p],
                    Internal_nblock_types,
                    MPI_UNSIGNED,
                    p,
                    3,
                    comm_pt()->mpi_comm(),
                    &s_req);
          send_requests_nrow.push_back(s_req);
 
          // recv
          nrows_to_recv[p] = new unsigned[Internal_nblock_types];
          MPI_Request r_req;
          MPI_Irecv(nrows_to_recv[p],
                    Internal_nblock_types,
                    MPI_UNSIGNED,
                    p,
                    3,
                    comm_pt()->mpi_comm(),
                    &r_req);
          recv_requests_nrow.push_back(r_req);
        }
        // send to self
        else
        {
          for (unsigned b = 0; b < Internal_nblock_types; b++)
          {
            Nrows_to_recv_for_get_block(b, p) =
              Nrows_to_send_for_get_block(b, p);
          }
          Nrows_to_recv_for_get_ordered[p] = Nrows_to_send_for_get_ordered[p];
        }
      }
 
      // create some temporary storage for the global row indices that will
      // be received from another processor.
      DenseMatrix<int*> block_rows_to_send(Internal_nblock_types, nproc, 0);
      Vector<int*> ordered_rows_to_send(nproc, 0);
 
      // resize the rows... storage
      Rows_to_send_for_get_block.resize(Internal_nblock_types, nproc);
      Rows_to_send_for_get_block.initialise(0);
      Rows_to_send_for_get_ordered.resize(nproc);
      Rows_to_send_for_get_ordered.initialise(0);
      Rows_to_recv_for_get_block.resize(Internal_nblock_types, nproc);
      Rows_to_recv_for_get_block.initialise(0);
 
      // resize the storage
      for (unsigned p = 0; p < nproc; p++)
      {
        for (unsigned b = 0; b < Internal_nblock_types; b++)
        {
          Rows_to_send_for_get_block(b, p) =
            new int[Nrows_to_send_for_get_block(b, p)];
          if (p != my_rank)
          {
            block_rows_to_send(b, p) =
              new int[Nrows_to_send_for_get_block(b, p)];
          }
          else
          {
            Rows_to_recv_for_get_block(b, p) =
              new int[Nrows_to_send_for_get_block(b, p)];
          }
        }
        Rows_to_send_for_get_ordered[p] =
          new int[Nrows_to_send_for_get_ordered[p]];
      }
 
 
      // loop over my rows to allocate the nrows
      DenseMatrix<unsigned> ptr_block(Internal_nblock_types, nproc, 0);
      for (unsigned i = 0; i < nrow_local; i++)
      {
        // the block number
        int b = this->internal_block_number(first_row + i);
 
        // check that the DOF i is associated with this preconditioner
        if (b >= 0)
        {
          // the block index
          unsigned j = this->internal_index_in_block(first_row + i);
 
          // the processor this row will be sent to
          unsigned block_p = 0;
          while (!(Internal_block_distribution_pt[b]->first_row(block_p) <= j &&
                   (Internal_block_distribution_pt[b]->first_row(block_p) +
                      Internal_block_distribution_pt[b]->nrow_local(block_p) >
                    j)))
          {
            block_p++;
          }
 
          // and store the row
          Rows_to_send_for_get_block(b, block_p)[ptr_block(b, block_p)] = i;
          if (block_p != my_rank)
          {
            block_rows_to_send(b, block_p)[ptr_block(b, block_p)] =
              j - Internal_block_distribution_pt[b]->first_row(block_p);
          }
          else
          {
            Rows_to_recv_for_get_block(b, block_p)[ptr_block(b, block_p)] =
              j - Internal_block_distribution_pt[b]->first_row(block_p);
          }
          ptr_block(b, block_p)++;
        }
      }
 
      // next block ordered
      for (unsigned p = 0; p < nproc; ++p)
      {
        int pt = 0;
        for (unsigned b = 0; b < Internal_nblock_types; ++b)
        {
          for (unsigned i = 0; i < Nrows_to_send_for_get_block(b, p); ++i)
          {
            Rows_to_send_for_get_ordered[p][pt] =
              Rows_to_send_for_get_block(b, p)[i];
            pt++;
          }
        }
      }
 
      // next process the nrow recvs as they complete
 
      // recv and store the data
      unsigned c = recv_requests_nrow.size();
      while (c > 0)
      {
        // wait for any communication to finish
        int req_number;
        MPI_Waitany(c, &recv_requests_nrow[0], &req_number, MPI_STATUS_IGNORE);
        recv_requests_nrow.erase(recv_requests_nrow.begin() + req_number);
        c--;
 
        // determine the source processor
        unsigned p = proc[req_number];
        proc.erase(proc.begin() + req_number);
 
        // copy the data to its final storage
        Nrows_to_recv_for_get_ordered[p] = 0;
        for (unsigned b = 0; b < Internal_nblock_types; b++)
        {
          Nrows_to_recv_for_get_block(b, p) = nrows_to_recv[p][b];
          Nrows_to_recv_for_get_ordered[p] += nrows_to_recv[p][b];
        }
 
        // and clear
        delete[] nrows_to_recv[p];
      }
 
      // resize the storage for the incoming rows data
      Rows_to_recv_for_get_ordered.resize(nproc, 0);
      for (unsigned p = 0; p < nproc; p++)
      {
        if (p != my_rank)
        {
          for (unsigned b = 0; b < Internal_nblock_types; b++)
          {
            Rows_to_recv_for_get_block(b, p) =
              new int[Nrows_to_recv_for_get_block(b, p)];
          }
        }
      }
 
      // compute the number of sends and recv from this processor
      // to each other processor
      Vector<unsigned> nsend_for_rows(nproc, 0);
      Vector<unsigned> nrecv_for_rows(nproc, 0);
      for (unsigned p = 0; p < nproc; p++)
      {
        if (p != my_rank)
        {
          for (unsigned b = 0; b < Internal_nblock_types; b++)
          {
            if (Nrows_to_send_for_get_block(b, p) > 0)
            {
              nsend_for_rows[p]++;
            }
            if (Nrows_to_recv_for_get_block(b, p) > 0)
            {
              nrecv_for_rows[p]++;
            }
          }
        }
      }
 
      // finally post the sends and recvs
      MPI_Aint base_displacement;
      MPI_Get_address(matrix_pt(), &base_displacement);
      Vector<MPI_Request> req_rows;
      for (unsigned p = 0; p < nproc; p++)
      {
        if (p != my_rank)
        {
          // send
          if (nsend_for_rows[p] > 0)
          {
            MPI_Datatype send_types[nsend_for_rows[p]];
            MPI_Aint send_displacements[nsend_for_rows[p]];
            int send_sz[nsend_for_rows[p]];
            unsigned send_ptr = 0;
            for (unsigned b = 0; b < Internal_nblock_types; b++)
            {
              if (Nrows_to_send_for_get_block(b, p) > 0)
              {
                MPI_Type_contiguous(Nrows_to_send_for_get_block(b, p),
                                    MPI_INT,
                                    &send_types[send_ptr]);
                MPI_Type_commit(&send_types[send_ptr]);
                MPI_Get_address(block_rows_to_send(b, p),
                                &send_displacements[send_ptr]);
                send_displacements[send_ptr] -= base_displacement;
                send_sz[send_ptr] = 1;
                send_ptr++;
              }
            }
            MPI_Datatype final_send_type;
            MPI_Type_create_struct(nsend_for_rows[p],
                                   send_sz,
                                   send_displacements,
                                   send_types,
                                   &final_send_type);
            MPI_Type_commit(&final_send_type);
            for (unsigned i = 0; i < nsend_for_rows[p]; i++)
            {
              MPI_Type_free(&send_types[i]);
            }
            MPI_Request send_req;
            MPI_Isend(matrix_pt(),
                      1,
                      final_send_type,
                      p,
                      4,
                      comm_pt()->mpi_comm(),
                      &send_req);
            req_rows.push_back(send_req);
            MPI_Type_free(&final_send_type);
          }
 
          // recv
          if (nrecv_for_rows[p] > 0)
          {
            MPI_Datatype recv_types[nrecv_for_rows[p]];
            MPI_Aint recv_displacements[nrecv_for_rows[p]];
            int recv_sz[nrecv_for_rows[p]];
            unsigned recv_ptr = 0;
            for (unsigned b = 0; b < Internal_nblock_types; b++)
            {
              if (Nrows_to_recv_for_get_block(b, p) > 0)
              {
                MPI_Type_contiguous(Nrows_to_recv_for_get_block(b, p),
                                    MPI_INT,
                                    &recv_types[recv_ptr]);
                MPI_Type_commit(&recv_types[recv_ptr]);
                MPI_Get_address(Rows_to_recv_for_get_block(b, p),
                                &recv_displacements[recv_ptr]);
                recv_displacements[recv_ptr] -= base_displacement;
                recv_sz[recv_ptr] = 1;
                recv_ptr++;
              }
            }
            MPI_Datatype final_recv_type;
            MPI_Type_create_struct(nrecv_for_rows[p],
                                   recv_sz,
                                   recv_displacements,
                                   recv_types,
                                   &final_recv_type);
            MPI_Type_commit(&final_recv_type);
            for (unsigned i = 0; i < nrecv_for_rows[p]; i++)
            {
              MPI_Type_free(&recv_types[i]);
            }
            MPI_Request recv_req;
            MPI_Irecv(matrix_pt(),
                      1,
                      final_recv_type,
                      p,
                      4,
                      comm_pt()->mpi_comm(),
                      &recv_req);
            req_rows.push_back(recv_req);
            MPI_Type_free(&final_recv_type);
          }
        }
      }
 
      // cleaning up Waitalls
 
 
      // wait for the recv requests so we can compute
      // Nrows_to_recv_for_get_ordered
      unsigned n_req_rows = req_rows.size();
      if (n_req_rows)
      {
        MPI_Waitall(n_req_rows, &req_rows[0], MPI_STATUS_IGNORE);
      }
 
      // resize the storage
      Rows_to_recv_for_get_ordered.resize(nproc);
      Rows_to_recv_for_get_ordered.initialise(0);
 
      // construct block offset
      Vector<int> vec_offset(Internal_nblock_types, 0);
      for (unsigned b = 1; b < Internal_nblock_types; ++b)
      {
        vec_offset[b] = vec_offset[b - 1] +
                        Internal_block_distribution_pt[b - 1]->nrow_local();
      }
 
      //
      for (unsigned p = 0; p < nproc; p++)
      {
        int pt = 0;
        Rows_to_recv_for_get_ordered[p] =
          new int[Nrows_to_recv_for_get_ordered[p]];
        for (unsigned b = 0; b < Internal_nblock_types; b++)
        {
          for (unsigned i = 0; i < Nrows_to_recv_for_get_block(b, p); i++)
          {
            Rows_to_recv_for_get_ordered[p][pt] =
              Rows_to_recv_for_get_block(b, p)[i] + vec_offset[b];
            pt++;
          }
        }
      }
 
      // clean up
      for (unsigned p = 0; p < nproc; p++)
      {
        if (p != my_rank)
        {
          for (unsigned b = 0; b < Internal_nblock_types; b++)
          {
            delete[] block_rows_to_send(b, p);
          }
          if (Nrows_to_send_for_get_ordered[p] > 0)
          {
            delete[] ordered_rows_to_send[p];
          }
        }
      }
 
      // and the send reqs
      unsigned n_req_send_nrow = send_requests_nrow.size();
      if (n_req_send_nrow)
      {
        MPI_Waitall(n_req_send_nrow, &send_requests_nrow[0], MPI_STATUS_IGNORE);
      }
      for (unsigned p = 0; p < nproc; p++)
      {
        delete[] nrows_to_send[p];
      }
#endif
    }
 
    // If we asked for output of blocks to a file then do it.
    if (block_output_on()) output_blocks_to_files(Output_base_filename);
  }
 
  //============================================================================
  //??ds
  /// Function to turn this preconditioner into a
  /// subsidiary preconditioner that operates within a bigger
  /// "master block preconditioner (e.g. a Navier-Stokes 2x2 block
  /// preconditioner dealing with the fluid sub-blocks within a
  /// 3x3 FSI preconditioner. Once this is done the master block
  /// preconditioner deals with the block setup etc.
  /// The vector block_map must specify the dof number in the
  /// master preconditioner that corresponds to a block number in this
  /// preconditioner. ??ds horribly misleading comment!
  /// The length of the vector is used to determine the number of
  /// blocks in this preconditioner therefore it must be correctly sized.
  /// This calls the other turn_into_subsidiary_block_preconditioner(...)
  /// function providing an empty doftype_to_doftype_map vector.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::turn_into_subsidiary_block_preconditioner(
    BlockPreconditioner<MATRIX>* master_block_prec_pt,
    const Vector<unsigned>& doftype_in_master_preconditioner_coarse)
  {
    // Create the identity dof_coarsen_map
    Vector<Vector<unsigned>> doftype_coarsen_map_coarse;
    unsigned doftype_in_master_preconditioner_coarse_size =
      doftype_in_master_preconditioner_coarse.size();
 
    for (unsigned dof_i = 0;
         dof_i < doftype_in_master_preconditioner_coarse_size;
         dof_i++)
    {
      // Create a vector of size 1 and value i,
      // then push it into the dof_coarsen_map vector.
      Vector<unsigned> tmp_vec(1, dof_i);
      doftype_coarsen_map_coarse.push_back(tmp_vec);
    }
 
    // Call the other turn_into_subsidiary_block_preconditioner function.
    turn_into_subsidiary_block_preconditioner(
      master_block_prec_pt,
      doftype_in_master_preconditioner_coarse,
      doftype_coarsen_map_coarse);
  }
 
 
  //============================================================================
  /// Function to turn this block preconditioner into a
  /// subsidiary block preconditioner that operates within a bigger
  /// master block preconditioner (e.g. a Navier-Stokes 2x2 block
  /// preconditioner dealing with the fluid sub-blocks within a
  /// 3x3 FSI preconditioner. Once this is done the master block
  /// preconditioner deals with the block setup etc.
  ///
  /// The vector doftype_map must specify the dof type in the
  /// master preconditioner that corresponds to a dof type in this block
  /// preconditioner.
  ///
  /// In general, we want:
  /// doftype_map[doftype in subsidiary prec] = doftype in master prec.
  ///
  /// It tells this block preconditioner which dof types of the master
  /// block preconditioner it is working with.
  ///
  /// The length of the vector is used to determine the number of
  /// dof types in THIS block preconditioner therefore it must be correctly
  /// sized.
  ///
  /// For example, let the master block preconditioner have 5 dof types in total
  /// and a 1-4 dof type splitting where the block (0,0) corresponds to
  /// dof type 0 and the block (1,1) corresponds to dof types 1, 2, 3 and 4
  /// (i.e. it would have given to block_setup the vector [0,1,1,1,1]).
  /// Furthermore, it solves (1,1) block with subsidiary block preconditioner.
  /// Then the doftype_map passed to this function of the subsidiary block
  /// preconditioner would be [1, 2, 3, 4].
  ///
  /// Dof type coarsening (following on from the example above):
  /// Let the subsidiary block preconditioner (THIS block preconditioner)
  /// only works with two DOF types, then the master block preconditioner must
  /// "coarsen" the dof types by providing the optional argument
  /// doftype_coarsen_map vector.
  ///
  /// The doftype_coarsen_map vector (in this case) might be [[0,1], [2,3]]
  /// telling the subsidiary block preconditioner that the SUBSIDIARY dof types
  /// 0 and 1 should be treated as dof type 0 and the subsidiary dof types 2
  /// and 3 should be treated as subsidiary dof type 1.
  ///
  /// If no doftype_coarsen_map vector is provided, then the identity is
  /// used automatically (see the turn_into_subsidiary_block_preconditioner(...)
  /// function with only two arguments). In the above case, the identity
  /// doftype_coarsen_map vector for the subsidiary block preconditioner
  /// would be the 2D vector [[0], [1], [2], [3]] which means
  /// dof type 0 is treated as dof type 0,
  /// dof type 1 is treated as dof type 1,
  /// dof type 2 is treated as dof type 2, and
  /// dof type 3 is treated as dof type 3.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::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)
  {
    // Set the master block preconditioner pointer
    Master_block_preconditioner_pt = master_block_prec_pt;
 
    // Set the Doftype_coarsen_map_coarse.
    Doftype_coarsen_map_coarse = doftype_coarsen_map_coarse;
 
    Doftype_in_master_preconditioner_coarse =
      doftype_in_master_preconditioner_coarse;
  } // end of turn_into_subsidiary_block_preconditioner(...)
 
 
  //============================================================================
  /// Determine the size of the matrix blocks and setup the
  /// lookup schemes relating the global degrees of freedom with
  /// their "blocks" and their indices (row/column numbers) in those
  /// blocks.
  /// The distributions of the preconditioner and the blocks are
  /// automatically specified (and assumed to be uniform) at this
  /// stage.
  /// This method should be used if each DOF type corresponds to a
  /// unique block type.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::block_setup()
  {
#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
 
    // Get the number of dof types.
    unsigned internal_n_dof_types = ndof_types();
 
    // Build the dof to block map - assume that each type of dof corresponds
    // to a different type of block.
    Vector<unsigned> dof_to_block_lookup(internal_n_dof_types);
    for (unsigned i = 0; i < internal_n_dof_types; i++)
    {
      dof_to_block_lookup[i] = i;
    }
 
    // call the block setup method
    this->block_setup(dof_to_block_lookup);
  }
 
 
  //============================================================================
  /// Get the block matrices required for the block preconditioner. Takes a
  /// pointer to a matrix of bools that indicate if a specified sub-block is
  /// required for the preconditioning operation. Computes the required block
  /// matrices, and stores pointers to them in the matrix block_matrix_pt. If an
  /// entry in block_matrix_pt is equal to NULL that sub-block has not been
  /// requested and is therefore not available.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::get_blocks(
    DenseMatrix<bool>& required_blocks,
    DenseMatrix<MATRIX*>& block_matrix_pt) const
  {
    // Cache number of block types
    const unsigned n_block_types = nblock_types();
 
#ifdef PARANOID
    // If required blocks matrix pointer is not the correct size then abort.
    if ((required_blocks.nrow() != n_block_types) ||
        (required_blocks.ncol() != n_block_types))
    {
      std::ostringstream error_message;
      error_message << "The size of the matrix of bools required_blocks "
                    << "(which indicates which blocks are required) is not the "
                    << "right size, required_blocks is "
                    << required_blocks.ncol() << " x " << required_blocks.nrow()
                    << ", whereas it should "
                    << "be " << n_block_types << " x " << n_block_types;
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // If block matrix pointer is not the correct size then abort.
    if ((block_matrix_pt.nrow() != n_block_types) ||
        (block_matrix_pt.ncol() != n_block_types))
    {
      std::ostringstream error_message;
      error_message << "The size of the block matrix pt is not the "
                    << "right size, block_matrix_pt is "
                    << block_matrix_pt.ncol() << " x " << block_matrix_pt.nrow()
                    << ", whereas it should "
                    << "be " << n_block_types << " x " << n_block_types;
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
#endif
 
    // Loop over the blocks
    for (unsigned i = 0; i < n_block_types; i++)
    {
      for (unsigned j = 0; j < n_block_types; j++)
      {
        // If block(i,j) is required then create a matrix and fill it in.
        if (required_blocks(i, j))
        {
          //??ds might want to remove this use of new as well?
          block_matrix_pt(i, j) = new MATRIX;
          get_block(i, j, *block_matrix_pt(i, j));
        }
 
        // Otherwise set pointer to null.
        else
        {
          block_matrix_pt(i, j) = 0;
        }
      }
    }
  }
 
  //============================================================================
  /// Takes the naturally ordered vector and extracts the blocks
  /// indicated by the block number (the values) in the Vector
  /// block_vec_number all at once, then concatenates them without
  /// communication. Here, the values in block_vec_number is the block number
  /// in the current preconditioner.
  /// This is a non-const function because distributions may be created
  /// and stored in Auxiliary_block_distribution_pt for future use.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::get_concatenated_block_vector(
    const Vector<unsigned>& block_vec_number,
    const DoubleVector& v,
    DoubleVector& w)
  {
#ifdef PARANOID
 
    // Check if v is built.
    if (!v.built())
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // v must have the same distribution as the upper-most master block
    // preconditioner, since the upper-most master block preconditioner
    // should have the same distribution as the matrix pointed to
    // by matrix_pt().
    if (*(v.distribution_pt()) != *(this->master_distribution_pt()))
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must match the "
              << " specified master_distribution_pt(). \n"
              << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check to see if there are more blocks defined in the block_vec_number
    // vector than the number of block types. This is not allowed.
    const unsigned para_nblock_types = nblock_types();
    const unsigned para_block_vec_number_size = block_vec_number.size();
    if (para_block_vec_number_size > para_nblock_types)
    {
      std::ostringstream err_msg;
      err_msg << "You have requested " << para_block_vec_number_size
              << " number of blocks, (block_vec_number.size() is "
              << para_block_vec_number_size << ").\n"
              << "But there are only " << para_nblock_types
              << " nblock_types.\n"
              << "Please make sure that block_vec_number is correctly sized.\n";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check if any block numbers defined in block_vec_number is equal to or
    // greater than the number of block types.
    // E.g. if there are 5 block types, we can only have block numbers:
    //  0, 1, 2, 3 and 4.
    for (unsigned i = 0; i < para_block_vec_number_size; i++)
    {
      const unsigned para_required_block = block_vec_number[i];
      if (para_required_block >= para_nblock_types)
      {
        std::ostringstream err_msg;
        err_msg << "block_vec_number[" << i << "] is " << para_required_block
                << ".\n"
                << "But there are only " << para_nblock_types
                << " nblock_types.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    // Check that no block number is inserted twice.
    std::set<unsigned> para_set;
    for (unsigned b = 0; b < para_block_vec_number_size; b++)
    {
      std::pair<std::set<unsigned>::iterator, bool> para_set_ret;
      para_set_ret = para_set.insert(block_vec_number[b]);
 
      if (!para_set_ret.second)
      {
        std::ostringstream err_msg;
        err_msg << "Error: the block number " << block_vec_number[b]
                << " appears twice.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // Number of blocks to get.
    const unsigned n_block = block_vec_number.size();
 
    // Each block is made of dof types. We get the most fine grain dof types.
    // Most fine grain in the sense that these are the dof types that belongs
    // in this block before any coarsening of dof types has taken place.
    // The ordering of the dof types matters, this is handled properly when
    // creating the Block_to_dof_map_fine vector and must be respected here.
    // I.e. we cannot arbitrarily insert dof types (even if they are correct)
    // in the vector most_fine_grain_dof.
    Vector<unsigned> most_fine_grain_dof;
    for (unsigned b = 0; b < n_block; b++)
    {
      const unsigned mapped_b = block_vec_number[b];
      most_fine_grain_dof.insert(most_fine_grain_dof.end(),
                                 Block_to_dof_map_fine[mapped_b].begin(),
                                 Block_to_dof_map_fine[mapped_b].end());
    }
 
    // Get all the dof level vectors in one go.
    Vector<DoubleVector> dof_block_vector;
    internal_get_block_vectors(most_fine_grain_dof, v, dof_block_vector);
 
    // Next we need to build the output DoubleVector w with the correct
    // distribution: the concatenation of the distributions of all the
    // dof-level vectors. This is the same as the concatenation of the
    // distributions of the blocks within this preconditioner.
    //
    // So we first check if it exists already, if not, we create it and
    // store it for future use. We store it because concatenation of
    // distributions requires communication, so concatenation of
    // distributions on-the-fly should be avoided.
    std::map<Vector<unsigned>, LinearAlgebraDistribution*>::const_iterator iter;
 
    // Attempt to get an iterator pointing to the pair with the value
    // block_vec_number.
    iter = Auxiliary_block_distribution_pt.find(block_vec_number);
 
    if (iter != Auxiliary_block_distribution_pt.end())
    // If it exists, build w with the distribution pointed to
    // by pair::second.
    {
      w.build(iter->second);
    }
    else
    // Else, we need to create the distribution and store it in
    // Auxiliary_block_distribution_pt.
    {
      Vector<LinearAlgebraDistribution*> tmp_vec_dist_pt(n_block, 0);
      for (unsigned b = 0; b < n_block; b++)
      {
        tmp_vec_dist_pt[b] = Block_distribution_pt[block_vec_number[b]];
      }
 
      // Note that the distribution is created with new but not deleted here.
      // This is handled in the clean up functions.
      LinearAlgebraDistribution* tmp_dist_pt = new LinearAlgebraDistribution;
      LinearAlgebraDistributionHelpers::concatenate(tmp_vec_dist_pt,
                                                    *tmp_dist_pt);
 
      // Store the pair of Vector<unsigned> and LinearAlgebraDistribution*
      insert_auxiliary_block_distribution(block_vec_number, tmp_dist_pt);
 
      // Build w.
      w.build(tmp_dist_pt);
    }
 
    // Now concatenate all the dof level vectors into the vector w.
    DoubleVectorHelpers::concatenate_without_communication(dof_block_vector, w);
 
  } // get_concatenated_block_vector(...)
 
  //============================================================================
  /// Takes concatenated block ordered vector, b, and copies its
  // entries to the appropriate entries in the naturally ordered vector, v.
  // Here the values in block_vec_number indicates which blocks the vector
  // b is a concatenation of. The block number are those in the current
  // preconditioner. If the preconditioner is a subsidiary block
  // preconditioner the other entries in v that are not associated with it
  // are left alone.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::return_concatenated_block_vector(
    const Vector<unsigned>& block_vec_number,
    const DoubleVector& w,
    DoubleVector& v) const
  {
#ifdef PARANOID
 
    // Check if v is built.
    if (!v.built())
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // v must have the same distribution as the upper-most master block
    // preconditioner, since the upper-most master block preconditioner
    // should have the same distribution as the matrix pointed to
    // by matrix_pt().
    if (*(v.distribution_pt()) != *(this->master_distribution_pt()))
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must match the "
              << " specified master_distribution_pt(). \n"
              << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check to see if there are more blocks defined in the block_vec_number
    // vector than the number of block types. This is not allowed.
    const unsigned para_block_vec_number_size = block_vec_number.size();
    const unsigned para_n_block = nblock_types();
    if (para_block_vec_number_size > para_n_block)
    {
      std::ostringstream err_msg;
      err_msg << "Trying to return " << para_block_vec_number_size
              << " block vectors.\n"
              << "But there are only " << para_n_block << " block types.\n";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check if any block numbers defined in block_vec_number is equal to or
    // greater than the number of block types.
    // E.g. if there are 5 block types, we can only have block numbers:
    //  0, 1, 2, 3 and 4.
    for (unsigned b = 0; b < para_block_vec_number_size; b++)
    {
      const unsigned para_required_block = block_vec_number[b];
      if (para_required_block > para_n_block)
      {
        std::ostringstream err_msg;
        err_msg << "block_vec_number[" << b << "] is " << para_required_block
                << ".\n"
                << "But there are only " << para_n_block << " block types.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    // Check that no block number is inserted twice.
    std::set<unsigned> para_set;
    for (unsigned b = 0; b < para_block_vec_number_size; b++)
    {
      std::pair<std::set<unsigned>::iterator, bool> para_set_ret;
      para_set_ret = para_set.insert(block_vec_number[b]);
 
      if (!para_set_ret.second)
      {
        std::ostringstream err_msg;
        err_msg << "Error: the block number " << block_vec_number[b]
                << " appears twice.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    // Check that w is built.
    if (!w.built())
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the block vector w must be setup.";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check that the distributions defined by block_vec_number is correct for
    // the distribution from w.
    // Recall that w is the concatenation of the block vectors defined by
    // the values in block_vec_number. We check that this is the case.
    Vector<LinearAlgebraDistribution*> para_vec_dist_pt(
      para_block_vec_number_size, 0);
 
    for (unsigned b = 0; b < para_block_vec_number_size; b++)
    {
      para_vec_dist_pt[b] = Block_distribution_pt[block_vec_number[b]];
    }
 
    LinearAlgebraDistribution para_tmp_dist;
 
    LinearAlgebraDistributionHelpers::concatenate(para_vec_dist_pt,
                                                  para_tmp_dist);
 
    if (*w.distribution_pt() != para_tmp_dist)
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the block vector w does not match \n"
              << "the concatenation of the block distributions defined in \n"
              << "block_vec_number.\n";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Number of blocks to return.
    const unsigned n_block = block_vec_number.size();
 
    // Each block is made of dof types. We get the most fine grain dof types.
    // Most fine grain in the sense that these are the dof types that belongs
    // in this block before any coarsening of dof types has taken place.
    // The ordering of the dof types matters, this is handled properly when
    // creating the Block_to_dof_map_fine vector and must be respected here.
    // I.e. we cannot arbitrarily insert dof types (even if they are correct)
    // in the vector most_fine_grain_dof.
    Vector<unsigned> most_fine_grain_dof;
    for (unsigned b = 0; b < n_block; b++)
    {
      const unsigned mapped_b = block_vec_number[b];
      most_fine_grain_dof.insert(most_fine_grain_dof.end(),
                                 Block_to_dof_map_fine[mapped_b].begin(),
                                 Block_to_dof_map_fine[mapped_b].end());
    }
 
    // The number of most fine grain dof types associated with the blocks
    // defined by block_vec_number.
    const unsigned ndof = most_fine_grain_dof.size();
 
    // Build each dof level vector with the correct distribution.
    Vector<DoubleVector> dof_vector(ndof);
    for (unsigned d = 0; d < ndof; d++)
    {
      dof_vector[d].build(
        internal_block_distribution_pt(most_fine_grain_dof[d]));
    }
 
    // Perform the splitting of w into the most fine grain dof level vectors.
    DoubleVectorHelpers::split_without_communication(w, dof_vector);
 
    // Return all the dof level vectors in one go.
    internal_return_block_vectors(most_fine_grain_dof, dof_vector, v);
  } // return_concatenated_block_vector(...)
 
  //============================================================================
  /// Takes the naturally ordered vector and rearranges it into a
  /// vector of sub vectors corresponding to the blocks, so s[b][i] contains
  /// the i-th entry in the vector associated with block b.
  /// Note: If the preconditioner is a subsidiary preconditioner then only the
  /// sub-vectors associated with the blocks of the subsidiary preconditioner
  /// will be included. Hence the length of v is master_nrow() whereas the
  /// total length of the s vectors is the sum of the lengths of the
  /// individual block vectors defined in block_vec_number.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::get_block_vectors(
    const Vector<unsigned>& block_vec_number,
    const DoubleVector& v,
    Vector<DoubleVector>& s) const
  {
#ifdef PARANOID
 
    // Check if v is built.
    if (!v.built())
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // v must have the same distribution as the upper-most master block
    // preconditioner, since the upper-most master block preconditioner
    // should have the same distribution as the matrix pointed to
    // by matrix_pt().
    if (*(v.distribution_pt()) != *(this->master_distribution_pt()))
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must match the "
              << " specified master_distribution_pt(). \n"
              << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check to see if there are more blocks defined in the block_vec_number
    // vector than the number of block types. This is not allowed.
    const unsigned para_nblock_types = nblock_types();
    const unsigned para_block_vec_number_size = block_vec_number.size();
    if (para_block_vec_number_size > para_nblock_types)
    {
      std::ostringstream err_msg;
      err_msg << "You have requested " << para_block_vec_number_size
              << " number of blocks, (block_vec_number.size() is "
              << para_block_vec_number_size << ").\n"
              << "But there are only " << para_nblock_types
              << " nblock_types.\n"
              << "Please make sure that block_vec_number is correctly sized.\n";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check if any block numbers defined in block_vec_number is equal to or
    // greater than the number of block types.
    // E.g. if there are 5 block types, we can only have block numbers:
    //  0, 1, 2, 3 and 4.
    for (unsigned i = 0; i < para_block_vec_number_size; i++)
    {
      const unsigned para_required_block = block_vec_number[i];
      if (para_required_block > para_nblock_types)
      {
        std::ostringstream err_msg;
        err_msg << "block_vec_number[" << i << "] is " << para_required_block
                << ".\n"
                << "But there are only " << para_nblock_types
                << " nblock_types.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
    // Check that no block number is inserted twice.
    std::set<unsigned> para_set;
    for (unsigned b = 0; b < para_block_vec_number_size; b++)
    {
      std::pair<std::set<unsigned>::iterator, bool> para_set_ret;
      para_set_ret = para_set.insert(block_vec_number[b]);
 
      if (!para_set_ret.second)
      {
        std::ostringstream err_msg;
        err_msg << "Error: the block number " << block_vec_number[b]
                << " appears twice.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // Number of blocks to get.
    const unsigned n_block = block_vec_number.size();
    s.resize(n_block);
 
    // Each block is made of dof types. We get the most fine grain dof types.
    // Most fine grain in the sense that these are the dof types that belongs
    // in this block before any coarsening of dof types has taken place.
    // The ordering of the dof types matters, this is handled properly when
    // creating the Block_to_dof_map_fine vector and must be respected here.
    // I.e. we cannot arbitrarily insert dof types (even if they are correct)
    // in the vector most_fine_grain_dof.
    Vector<unsigned> most_fine_grain_dof;
    for (unsigned b = 0; b < n_block; b++)
    {
      const unsigned mapped_b = block_vec_number[b];
 
      most_fine_grain_dof.insert(most_fine_grain_dof.end(),
                                 Block_to_dof_map_fine[mapped_b].begin(),
                                 Block_to_dof_map_fine[mapped_b].end());
    }
 
    // Get all the dof level vectors in one go.
    Vector<DoubleVector> dof_vector;
    internal_get_block_vectors(most_fine_grain_dof, v, dof_vector);
 
    // For each block vector requested,
    // build the block s[b],
    // concatenate the corresponding dof vector
 
    // Since all the dof vectors are in dof_vector,
    // we need to loop through this.
    // The offset helps us loop through this.
    unsigned offset = 0;
 
    for (unsigned b = 0; b < n_block; b++)
    {
      // The actual block number required.
      const unsigned mapped_b = block_vec_number[b];
 
      // How many most fine grain dofs are in this block?
      const unsigned n_dof = Block_to_dof_map_fine[mapped_b].size();
 
      if (n_dof == 1)
      // No need to concatenate, just copy the DoubleVector.
      {
        s[b] = dof_vector[offset];
      }
      else
      // Concatenate the relevant dof vectors into s[b].
      {
        s[b].build(Block_distribution_pt[mapped_b], 0);
        Vector<DoubleVector*> tmp_vec_pt(n_dof, 0);
        for (unsigned vec_i = 0; vec_i < n_dof; vec_i++)
        {
          tmp_vec_pt[vec_i] = &dof_vector[offset + vec_i];
        }
 
        DoubleVectorHelpers::concatenate_without_communication(tmp_vec_pt,
                                                               s[b]);
      }
 
      // Update the offset.
      offset += n_dof;
    }
  } // get_block_vectors(...)
 
 
  //============================================================================
  /// Takes the naturally ordered vector and rearranges it into a
  /// vector of sub vectors corresponding to the blocks, so s[b][i] contains
  /// the i-th entry in the vector associated with block b.
  /// Note: If the preconditioner is a subsidiary preconditioner then only the
  /// sub-vectors associated with the blocks of the subsidiary preconditioner
  /// will be included. Hence the length of v is master_nrow() whereas the
  /// total length of the s vectors is Nrow.
  /// This is simply a wrapper around the other get_block_vectors(...) function
  /// where the block_vec_number Vector is the identity, i.e.
  /// block_vec_number is [0, 1, ..., nblock_types - 1].
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::get_block_vectors(
    const DoubleVector& v, Vector<DoubleVector>& s) const
  {
    // Get the number of blocks in this block preconditioner.
    const unsigned n_block = nblock_types();
 
    // Create the identity vector.
    Vector<unsigned> required_block(n_block, 0);
    for (unsigned i = 0; i < n_block; i++)
    {
      required_block[i] = i;
    }
 
    // Call the other function which does the work.
    get_block_vectors(required_block, v, s);
  }
 
  //============================================================================
  /// Takes the naturally ordered vector and
  /// rearranges it into a vector of sub vectors corresponding to the blocks,
  /// so s[b][i] contains the i-th entry in the vector associated with block b.
  /// The block_vec_number indicates which blocks we want.
  /// These blocks and vectors are those corresponding to the internal blocks.
  /// Note: If the preconditioner is a subsidiary preconditioner then only the
  /// sub-vectors associated with the blocks of the subsidiary preconditioner
  /// will be included. Hence the length of v is master_nrow() whereas the
  /// total length of the s vectors is the sum of the Nrow of the sub vectors.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::internal_get_block_vectors(
    const Vector<unsigned>& block_vec_number,
    const DoubleVector& v,
    Vector<DoubleVector>& s) const
  {
#ifdef PARANOID
    if (!v.built())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*(v.distribution_pt()) != *(this->master_distribution_pt()))
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must match the "
                    << " specified master_distribution_pt(). \n"
                    << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Number of block types
    // const unsigned nblock = this->internal_nblock_types();
    const unsigned nblock = block_vec_number.size();
 
    // if + only one processor
    //    + more than one processor but matrix_pt is not distributed
    // then use the serial get_block method
    if (this->distribution_pt()->communicator_pt()->nproc() == 1 ||
        !this->distribution_pt()->distributed())
    {
      // Vector of vectors for each section of residual vector
      s.resize(nblock);
 
      // pointer to the data in v
      const double* v_pt = v.values_pt();
 
      // setup the block vector and then insert the data
      for (unsigned b = 0; b < nblock; b++)
      {
        const unsigned required_block = block_vec_number[b];
        s[b].build(Internal_block_distribution_pt[required_block], 0.0);
        double* s_pt = s[b].values_pt();
        unsigned nrow = s[b].nrow();
        for (unsigned i = 0; i < nrow; i++)
        {
          s_pt[i] = v_pt[this->Global_index[required_block][i]];
        }
      }
    }
    // otherwise use mpi
    else
    {
#ifdef OOMPH_HAS_MPI
      // my rank
      unsigned my_rank = this->distribution_pt()->communicator_pt()->my_rank();
 
      // the number of processors
      unsigned nproc = this->distribution_pt()->communicator_pt()->nproc();
 
      // build the vectors
      s.resize(nblock);
      for (unsigned b = 0; b < nblock; b++)
      {
        const unsigned required_block = block_vec_number[b];
        s[b].build(Internal_block_distribution_pt[required_block], 0.0);
      }
 
      // determine the maximum number of rows to be sent or recv
      // and determine the number of blocks each processor will send and recv
      // communication for
      Vector<int> nblock_send(nproc, 0);
      Vector<int> nblock_recv(nproc, 0);
      unsigned max_n_send_or_recv = 0;
      for (unsigned p = 0; p < nproc; p++)
      {
        for (unsigned b = 0; b < nblock; b++)
        {
          const unsigned required_block = block_vec_number[b];
          max_n_send_or_recv = std::max(
            max_n_send_or_recv, Nrows_to_send_for_get_block(required_block, p));
          max_n_send_or_recv = std::max(
            max_n_send_or_recv, Nrows_to_recv_for_get_block(required_block, p));
          if (Nrows_to_send_for_get_block(required_block, p) > 0)
          {
            nblock_send[p]++;
          }
          if (Nrows_to_recv_for_get_block(required_block, p) > 0)
          {
            nblock_recv[p]++;
          }
        }
      }
 
      // create a vectors of 1s the size of the nblock for the mpi indexed
      // data types
      int* block_lengths = new int[max_n_send_or_recv];
      for (unsigned i = 0; i < max_n_send_or_recv; i++)
      {
        block_lengths[i] = 1;
      }
 
      // perform the sends and receives
      Vector<MPI_Request> requests;
      for (unsigned p = 0; p < nproc; p++)
      {
        // send and recv with other processors
        if (p != my_rank)
        {
          // send
          if (nblock_send[p] > 0)
          {
            // create the datatypes vector
            MPI_Datatype block_send_types[nblock_send[p]];
 
            // create the datatypes
            unsigned ptr = 0;
            for (unsigned b = 0; b < nblock; b++)
            {
              const unsigned required_block = block_vec_number[b];
 
              if (Nrows_to_send_for_get_block(required_block, p) > 0)
              {
                MPI_Type_indexed(Nrows_to_send_for_get_block(required_block, p),
                                 block_lengths,
                                 Rows_to_send_for_get_block(required_block, p),
                                 MPI_DOUBLE,
                                 &block_send_types[ptr]);
                MPI_Type_commit(&block_send_types[ptr]);
                ptr++;
              }
            }
 
            // compute the displacements and lengths
            MPI_Aint displacements[nblock_send[p]];
            int lengths[nblock_send[p]];
            for (int i = 0; i < nblock_send[p]; i++)
            {
              lengths[i] = 1;
              displacements[i] = 0;
            }
 
            // build the final datatype
            MPI_Datatype type_send;
            MPI_Type_create_struct(nblock_send[p],
                                   lengths,
                                   displacements,
                                   block_send_types,
                                   &type_send);
            MPI_Type_commit(&type_send);
 
            // send
            MPI_Request send_req;
            MPI_Isend(const_cast<double*>(v.values_pt()),
                      1,
                      type_send,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &send_req);
            MPI_Type_free(&type_send);
            for (int i = 0; i < nblock_send[p]; i++)
            {
              MPI_Type_free(&block_send_types[i]);
            }
            requests.push_back(send_req);
          }
 
          // recv
          if (nblock_recv[p] > 0)
          {
            // create the datatypes vector
            MPI_Datatype block_recv_types[nblock_recv[p]];
 
            // and the displacements
            MPI_Aint displacements[nblock_recv[p]];
 
            // and the lengths
            int lengths[nblock_recv[p]];
 
            // all displacements are computed relative to s[0] values
            MPI_Aint displacements_base;
            MPI_Get_address(s[0].values_pt(), &displacements_base);
 
            // now build
            unsigned ptr = 0;
            for (unsigned b = 0; b < nblock; b++)
            {
              const unsigned required_block = block_vec_number[b];
 
              if (Nrows_to_recv_for_get_block(required_block, p) > 0)
              {
                MPI_Type_indexed(Nrows_to_recv_for_get_block(required_block, p),
                                 block_lengths,
                                 Rows_to_recv_for_get_block(required_block, p),
                                 MPI_DOUBLE,
                                 &block_recv_types[ptr]);
                MPI_Type_commit(&block_recv_types[ptr]);
                MPI_Get_address(s[b].values_pt(), &displacements[ptr]);
                displacements[ptr] -= displacements_base;
                lengths[ptr] = 1;
                ptr++;
              }
            }
 
            // build the final data type
            MPI_Datatype type_recv;
            MPI_Type_create_struct(nblock_recv[p],
                                   lengths,
                                   displacements,
                                   block_recv_types,
                                   &type_recv);
            MPI_Type_commit(&type_recv);
 
            // recv
            MPI_Request recv_req;
            MPI_Irecv(s[0].values_pt(),
                      1,
                      type_recv,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &recv_req);
            MPI_Type_free(&type_recv);
            for (int i = 0; i < nblock_recv[p]; i++)
            {
              MPI_Type_free(&block_recv_types[i]);
            }
            requests.push_back(recv_req);
          }
        }
 
        // communicate with self
        else
        {
          const double* v_values_pt = v.values_pt();
          for (unsigned b = 0; b < nblock; b++)
          {
            const unsigned required_block = block_vec_number[b];
 
            double* w_values_pt = s[b].values_pt();
            for (unsigned i = 0;
                 i < Nrows_to_send_for_get_block(required_block, p);
                 i++)
            {
              w_values_pt[Rows_to_recv_for_get_block(required_block, p)[i]] =
                v_values_pt[Rows_to_send_for_get_block(required_block, p)[i]];
            }
          }
        }
      }
 
      // and then just wait
      unsigned c = requests.size();
      Vector<MPI_Status> stat(c);
      if (c)
      {
        MPI_Waitall(c, &requests[0], &stat[0]);
      }
      delete[] block_lengths;
 
#else
      // throw error
      std::ostringstream error_message;
      error_message << "The preconditioner is distributed and on more than one "
                    << "processor. MPI is required.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
#endif
    }
  }
 
  //============================================================================
  /// A helper function, takes the naturally ordered vector and
  /// rearranges it into a vector of sub vectors corresponding to the blocks,
  /// so s[b][i] contains the i-th entry in the vector associated with block b.
  /// The block_vec_number indicates which blocks we want.
  /// These blocks and vectors are those corresponding to the internal blocks.
  /// Note: If the preconditioner is a subsidiary preconditioner then only the
  /// sub-vectors associated with the blocks of the subsidiary preconditioner
  /// will be included. Hence the length of v is master_nrow() whereas the
  /// total length of the s vectors is the sum of the Nrow of the sub vectors.
  /// This is simply a wrapper around the other internal_get_block_vectors(...)
  /// function with the identity block_vec_number vector.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::internal_get_block_vectors(
    const DoubleVector& v, Vector<DoubleVector>& s) const
  {
    // Number of block types
    const unsigned nblock = this->internal_nblock_types();
    Vector<unsigned> block_vec_number(nblock, 0);
    for (unsigned b = 0; b < nblock; b++)
    {
      block_vec_number[b] = b;
    }
 
    internal_get_block_vectors(block_vec_number, v, s);
  }
 
  //============================================================================
  /// Takes the vector of block vectors, s, and copies its entries into
  /// the naturally ordered vector, v. If this is a subsidiary block
  /// preconditioner only those entries in v that are associated with its
  /// blocks are affected. The block_vec_number indicates which block the
  /// vectors in s came from. The block number corresponds to the block
  /// numbers in this preconditioner.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::return_block_vectors(
    const Vector<unsigned>& block_vec_number,
    const Vector<DoubleVector>& s,
    DoubleVector& v) const
  {
#ifdef PARANOID
 
    // Check if v is built.
    if (!v.built())
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // v must have the same distribution as the upper-most master block
    // preconditioner, since the upper-most master block preconditioner
    // should have the same distribution as the matrix pointed to
    // by matrix_pt().
    if (*(v.distribution_pt()) != *(this->master_distribution_pt()))
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must match the "
              << " specified master_distribution_pt(). \n"
              << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check if the number of vectors in s is the same as the number of block
    // numbers described in block_vec_number.
    const unsigned para_block_vec_number_size = block_vec_number.size();
    const unsigned para_s_size = s.size();
    if (para_block_vec_number_size != para_s_size)
    {
      std::ostringstream err_msg;
      err_msg << "block_vec_number.size() is " << para_block_vec_number_size
              << "\n."
              << "s.size() is " << para_s_size << ".\n"
              << "But they must be the same size!\n";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check to see if there are more blocks defined in the block_vec_number
    // vector than the number of block types. This is not allowed.
    const unsigned para_n_block = nblock_types();
    if (para_block_vec_number_size > para_n_block)
    {
      std::ostringstream err_msg;
      err_msg << "Trying to return " << para_block_vec_number_size
              << " block vectors.\n"
              << "But there are only " << para_n_block << " block types.\n";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check if any block numbers defined in block_vec_number is equal to or
    // greater than the number of block types.
    // E.g. if there are 5 block types, we can only have block numbers:
    //  0, 1, 2, 3 and 4.
    for (unsigned b = 0; b < para_block_vec_number_size; b++)
    {
      const unsigned para_required_block = block_vec_number[b];
      if (para_required_block > para_n_block)
      {
        std::ostringstream err_msg;
        err_msg << "block_vec_number[" << b << "] is " << para_required_block
                << ".\n"
                << "But there are only " << para_n_block << " block types.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    // Check that no block number is inserted twice.
    std::set<unsigned> para_set;
    for (unsigned b = 0; b < para_block_vec_number_size; b++)
    {
      std::pair<std::set<unsigned>::iterator, bool> para_set_ret;
      para_set_ret = para_set.insert(block_vec_number[b]);
 
      if (!para_set_ret.second)
      {
        std::ostringstream err_msg;
        err_msg << "Error: the block number " << block_vec_number[b]
                << " appears twice.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    // Check to see that all the vectors in s are built
    // (since we are trying to return them).
    for (unsigned b = 0; b < para_block_vec_number_size; b++)
    {
      if (!s[b].built())
      {
        std::ostringstream err_msg;
        err_msg << "The distribution of the block vector s[" << b
                << "] must be setup.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
    // Since these are built, we check that the distributions are correct.
    // This are incorrect if the block numbers in block_vec_number and
    // the vectors in s does not match.
    for (unsigned b = 0; b < para_block_vec_number_size; b++)
    {
      if (*(s[b].distribution_pt()) !=
          *(Block_distribution_pt[block_vec_number[b]]))
      {
        std::ostringstream error_message;
        error_message
          << "The distribution of the block vector " << b << " must match the"
          << " specified distribution at "
          << "Block_distribution_pt[" << block_vec_number[b] << "].\n"
          << "The distribution of the Block_distribution_pt is determined by\n"
          << "the vector block_vec_number. Perhaps it is incorrect?\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // Number of blocks to get.
    const unsigned n_block = block_vec_number.size();
 
    // Each block is made of dof types. We get the most fine grain dof types.
    // Most fine grain in the sense that these are the dof types that belongs
    // in this block before any coarsening of dof types has taken place.
    // The ordering of the dof types matters, this is handled properly when
    // creating the Block_to_dof_map_fine vector and must be respected here.
    // I.e. we cannot arbitrarily insert dof types (even if they are correct)
    // in the vector most_fine_grain_dof.
    Vector<unsigned> most_fine_grain_dof;
    for (unsigned b = 0; b < n_block; b++)
    {
      const unsigned mapped_b = block_vec_number[b];
 
      most_fine_grain_dof.insert(most_fine_grain_dof.end(),
                                 Block_to_dof_map_fine[mapped_b].begin(),
                                 Block_to_dof_map_fine[mapped_b].end());
    }
 
    // Split all the blocks into it's most fine grain dof vector.
    Vector<DoubleVector> dof_vector(most_fine_grain_dof.size());
 
    unsigned offset = 0;
 
    // Perform the splitting for each block.
    for (unsigned b = 0; b < n_block; b++)
    {
      // The actual block number.
      const unsigned mapped_b = block_vec_number[b];
 
      // How many most fine grain dof types are associated with this block?
      const unsigned ndof = Block_to_dof_map_fine[mapped_b].size();
 
      if (ndof == 1)
      // No need to split, just copy.
      {
        dof_vector[offset] = s[b];
      }
      else
      // Need to split s[b] into it's most fine grain dof vectors
      {
        // To store pointers to the dof vectors associated with this block.
        Vector<DoubleVector*> tmp_dof_vector_pt(ndof, 0);
 
        for (unsigned d = 0; d < ndof; d++)
        {
          const unsigned offset_plus_d = offset + d;
 
          // build the dof vector.
          dof_vector[offset_plus_d].build(
            Internal_block_distribution_pt[most_fine_grain_dof[offset_plus_d]]);
 
          // Store the pointer.
          tmp_dof_vector_pt[d] = &dof_vector[offset_plus_d];
        }
 
        // Split without communication.
        DoubleVectorHelpers::split_without_communication(s[b],
                                                         tmp_dof_vector_pt);
      }
 
      // Update the offset!
      offset += ndof;
    }
 
    // Return the block vectors all in one go.
    internal_return_block_vectors(most_fine_grain_dof, dof_vector, v);
  } // return_block_vectors(...)
 
 
  //============================================================================
  /// Takes the vector of block vectors, s, and copies its entries into
  /// the naturally ordered vector, v. If this is a subsidiary block
  /// preconditioner only those entries in v that are associated with its
  /// blocks are affected. The block_vec_number indicates which block the
  /// vectors in s came from. The block number corresponds to the block
  /// numbers in this preconditioner.
  /// This is simply a wrapper around the other return_block_vectors(...)
  /// function where the block_vec_number Vector is the identity, i.e.
  /// block_vec_number is [0, 1, ..., nblock_types - 1].
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::return_block_vectors(
    const Vector<DoubleVector>& s, DoubleVector& v) const
  {
    // The number of block types in this preconditioner.
    const unsigned n_block = nblock_types();
 
    // Create the identity vector.
    Vector<unsigned> required_block(n_block, 0);
    for (unsigned i = 0; i < n_block; i++)
    {
      required_block[i] = i;
    }
 
    // Call the other return_block_vectors function which does the work.
    return_block_vectors(required_block, s, v);
  } // return_block_vectors(...)
 
  //============================================================================
  /// Takes the naturally ordered vector and
  /// rearranges it into a vector of sub vectors corresponding to the blocks,
  /// so s[b][i] contains the i-th entry in the vector associated with block b.
  /// The block_vec_number indicates which blocks we want.
  /// These blocks and vectors are those corresponding to the internal blocks.
  /// Note: If the preconditioner is a subsidiary preconditioner then only the
  /// sub-vectors associated with the blocks of the subsidiary preconditioner
  /// will be included. Hence the length of v is master_nrow() whereas the
  /// total length of the s vectors is the sum of the Nrow of the sub vectors.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::internal_return_block_vectors(
    const Vector<unsigned>& block_vec_number,
    const Vector<DoubleVector>& s,
    DoubleVector& v) const
  {
    // the number of blocks
    const unsigned nblock = block_vec_number.size();
 
#ifdef PARANOID
    if (!v.built())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*(v.distribution_pt()) != *(this->master_distribution_pt()))
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must match the "
                    << " specified master_distribution_pt(). \n"
                    << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    for (unsigned b = 0; b < nblock; b++)
    {
      if (!s[b].built())
      {
        std::ostringstream error_message;
        error_message << "The distribution of the block vector " << b
                      << " must be setup.";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      const unsigned required_block = block_vec_number[b];
      if (*(s[b].distribution_pt()) !=
          *(Internal_block_distribution_pt[required_block]))
      {
        std::ostringstream error_message;
        error_message
          << "The distribution of the block vector " << b << " must match the"
          << " specified distribution at Internal_block_distribution_pt[" << b
          << "]";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // if + only one processor
    //    + more than one processor but matrix_pt is not distributed
    // then use the serial get_block method
    if (this->distribution_pt()->communicator_pt()->nproc() == 1 ||
        !this->distribution_pt()->distributed())
    {
      double* v_pt = v.values_pt();
      for (unsigned b = 0; b < nblock; b++)
      {
        const unsigned required_block = block_vec_number[b];
 
        const double* s_pt = s[b].values_pt();
        unsigned nrow = this->internal_block_dimension(required_block);
        for (unsigned i = 0; i < nrow; i++)
        {
          v_pt[this->Global_index[required_block][i]] = s_pt[i];
        }
      }
    }
    // otherwise use mpi
    else
    {
#ifdef OOMPH_HAS_MPI
 
      // my rank
      unsigned my_rank = this->distribution_pt()->communicator_pt()->my_rank();
 
      // the number of processors
      unsigned nproc = this->distribution_pt()->communicator_pt()->nproc();
 
      // determine the maximum number of rows to be sent or recv
      // and determine the number of blocks each processor will send and recv
      // communication for
      Vector<int> nblock_send(nproc, 0);
      Vector<int> nblock_recv(nproc, 0);
      unsigned max_n_send_or_recv = 0;
      for (unsigned p = 0; p < nproc; p++)
      {
        for (unsigned b = 0; b < nblock; b++)
        {
          const unsigned required_block = block_vec_number[b];
 
          max_n_send_or_recv = std::max(
            max_n_send_or_recv, Nrows_to_send_for_get_block(required_block, p));
          max_n_send_or_recv = std::max(
            max_n_send_or_recv, Nrows_to_recv_for_get_block(required_block, p));
          if (Nrows_to_send_for_get_block(required_block, p) > 0)
          {
            nblock_recv[p]++;
          }
          if (Nrows_to_recv_for_get_block(required_block, p) > 0)
          {
            nblock_send[p]++;
          }
        }
      }
 
      // create a vectors of 1s the size of the nblock for the mpi indexed
      // data types
      int* block_lengths = new int[max_n_send_or_recv];
      for (unsigned i = 0; i < max_n_send_or_recv; i++)
      {
        block_lengths[i] = 1;
      }
 
      // perform the sends and receives
      Vector<MPI_Request> requests;
      for (unsigned p = 0; p < nproc; p++)
      {
        // send and recv with other processors
        if (p != my_rank)
        {
          // recv
          if (nblock_recv[p] > 0)
          {
            // create the datatypes vector
            MPI_Datatype block_recv_types[nblock_recv[p]];
 
            // create the datatypes
            unsigned ptr = 0;
            for (unsigned b = 0; b < nblock; b++)
            {
              const unsigned required_block = block_vec_number[b];
 
              if (Nrows_to_send_for_get_block(required_block, p) > 0)
              {
                MPI_Type_indexed(Nrows_to_send_for_get_block(required_block, p),
                                 block_lengths,
                                 Rows_to_send_for_get_block(required_block, p),
                                 MPI_DOUBLE,
                                 &block_recv_types[ptr]);
                MPI_Type_commit(&block_recv_types[ptr]);
                ptr++;
              }
            }
 
            // compute the displacements and lengths
            MPI_Aint displacements[nblock_recv[p]];
            int lengths[nblock_recv[p]];
            for (int i = 0; i < nblock_recv[p]; i++)
            {
              lengths[i] = 1;
              displacements[i] = 0;
            }
 
            // build the final datatype
            MPI_Datatype type_recv;
            MPI_Type_create_struct(nblock_recv[p],
                                   lengths,
                                   displacements,
                                   block_recv_types,
                                   &type_recv);
            MPI_Type_commit(&type_recv);
 
            // recv
            MPI_Request recv_req;
            MPI_Irecv(v.values_pt(),
                      1,
                      type_recv,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &recv_req);
            MPI_Type_free(&type_recv);
            for (int i = 0; i < nblock_recv[p]; i++)
            {
              MPI_Type_free(&block_recv_types[i]);
            }
            requests.push_back(recv_req);
          }
 
          // send
          if (nblock_send[p] > 0)
          {
            // create the datatypes vector
            MPI_Datatype block_send_types[nblock_send[p]];
 
            // and the displacements
            MPI_Aint displacements[nblock_send[p]];
 
            // and the lengths
            int lengths[nblock_send[p]];
 
            // all displacements are computed relative to s[0] values
            MPI_Aint displacements_base;
            MPI_Get_address(const_cast<double*>(s[0].values_pt()),
                            &displacements_base);
 
            // now build
            unsigned ptr = 0;
            for (unsigned b = 0; b < nblock; b++)
            {
              const unsigned required_block = block_vec_number[b];
 
              if (Nrows_to_recv_for_get_block(required_block, p) > 0)
              {
                MPI_Type_indexed(Nrows_to_recv_for_get_block(required_block, p),
                                 block_lengths,
                                 Rows_to_recv_for_get_block(required_block, p),
                                 MPI_DOUBLE,
                                 &block_send_types[ptr]);
                MPI_Type_commit(&block_send_types[ptr]);
                MPI_Get_address(const_cast<double*>(s[b].values_pt()),
                                &displacements[ptr]);
                displacements[ptr] -= displacements_base;
                lengths[ptr] = 1;
                ptr++;
              }
            }
 
            // build the final data type
            MPI_Datatype type_send;
            MPI_Type_create_struct(nblock_send[p],
                                   lengths,
                                   displacements,
                                   block_send_types,
                                   &type_send);
            MPI_Type_commit(&type_send);
 
            // send
            MPI_Request send_req;
            MPI_Isend(const_cast<double*>(s[0].values_pt()),
                      1,
                      type_send,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &send_req);
            MPI_Type_free(&type_send);
            for (int i = 0; i < nblock_send[p]; i++)
            {
              MPI_Type_free(&block_send_types[i]);
            }
            requests.push_back(send_req);
          }
        }
 
        // communicate wih self
        else
        {
          double* v_values_pt = v.values_pt();
          for (unsigned b = 0; b < nblock; b++)
          {
            const unsigned required_block = block_vec_number[b];
 
            const double* w_values_pt = s[b].values_pt();
            for (unsigned i = 0;
                 i < Nrows_to_send_for_get_block(required_block, p);
                 i++)
            {
              v_values_pt[Rows_to_send_for_get_block(required_block, p)[i]] =
                w_values_pt[Rows_to_recv_for_get_block(required_block, p)[i]];
            }
          }
        }
      }
 
      // and then just wait
      unsigned c = requests.size();
      Vector<MPI_Status> stat(c);
      if (c)
      {
        MPI_Waitall(c, &requests[0], &stat[0]);
      }
      delete[] block_lengths;
 
#else
      // throw error
      std::ostringstream error_message;
      error_message << "The preconditioner is distributed and on more than one "
                    << "processor. MPI is required.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
#endif
    }
  }
 
  //============================================================================
  /// A helper function, takes the naturally ordered vector and
  /// rearranges it into a vector of sub vectors corresponding to the blocks,
  /// so s[b][i] contains the i-th entry in the vector associated with block b.
  /// The block_vec_number indicates which blocks we want.
  /// These blocks and vectors are those corresponding to the internal blocks.
  /// Note: If the preconditioner is a subsidiary preconditioner then only the
  /// sub-vectors associated with the blocks of the subsidiary preconditioner
  /// will be included. Hence the length of v is master_nrow() whereas the
  /// total length of the s vectors is the sum of the Nrow of the sub vectors.
  /// This is simply a wrapper around the other internal_get_block_vectors(...)
  /// function with the identity block_vec_number vector.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::internal_return_block_vectors(
    const Vector<DoubleVector>& s, DoubleVector& v) const
  {
    // the number of blocks
    const unsigned nblock = this->internal_nblock_types();
    Vector<unsigned> block_vec_number(nblock, 0);
    for (unsigned b = 0; b < nblock; b++)
    {
      block_vec_number[b] = b;
    }
 
    internal_return_block_vectors(block_vec_number, s, v);
  }
 
  //============================================================================
  /// A helper function, takes the naturally ordered vector, v,
  /// and extracts the n-th block vector, b.
  /// Here n is the block number in the current preconditioner.
  /// NOTE: The ordering of the vector b is the same as the
  /// ordering of the block matrix from internal_get_block(...).
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::internal_get_block_vector(
    const unsigned& b, const DoubleVector& v, DoubleVector& w) const
  {
#ifdef PARANOID
    // the number of blocks
    const unsigned n_blocks = this->internal_nblock_types();
 
    // paranoid check that block i is in this block preconditioner
    if (b >= n_blocks)
    {
      std::ostringstream error_message;
      error_message
        << "Requested block  vector " << b
        << ", however this preconditioner has internal_nblock_types() "
        << "= " << internal_nblock_types() << std::endl;
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (!v.built())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*(v.distribution_pt()) != *(this->master_distribution_pt()))
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must match the "
                    << " specified master_distribution_pt(). \n"
                    << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // rebuild the block vector
    w.build(Internal_block_distribution_pt[b], 0.0);
 
    // if + only one processor
    //    + more than one processor but matrix_pt is not distributed
    // then use the serial get_block method
    if (this->distribution_pt()->communicator_pt()->nproc() == 1 ||
        !this->distribution_pt()->distributed())
    {
      double* w_pt = w.values_pt();
      const double* v_pt = v.values_pt();
      unsigned n_row = w.nrow();
      for (unsigned i = 0; i < n_row; i++)
      {
        w_pt[i] = v_pt[this->Global_index[b][i]];
      }
    }
    // otherwise use mpi
    else
    {
#ifdef OOMPH_HAS_MPI
 
      // my rank
      unsigned my_rank = this->distribution_pt()->communicator_pt()->my_rank();
 
      // the number of processors
      unsigned nproc = this->distribution_pt()->communicator_pt()->nproc();
 
      // determine the maximum number of rows to be sent or recv
      unsigned max_n_send_or_recv = 0;
      for (unsigned p = 0; p < nproc; p++)
      {
        max_n_send_or_recv =
          std::max(max_n_send_or_recv, Nrows_to_send_for_get_block(b, p));
        max_n_send_or_recv =
          std::max(max_n_send_or_recv, Nrows_to_recv_for_get_block(b, p));
      }
 
      // create a vectors of 1s (the size of the nblock for the mpi indexed
      // data types
      int* block_lengths = new int[max_n_send_or_recv];
      for (unsigned i = 0; i < max_n_send_or_recv; i++)
      {
        block_lengths[i] = 1;
      }
 
      // perform the sends and receives
      Vector<MPI_Request> requests;
      for (unsigned p = 0; p < nproc; p++)
      {
        // send and recv with other processors
        if (p != my_rank)
        {
          if (Nrows_to_send_for_get_block(b, p) > 0)
          {
            // create the send datatype
            MPI_Datatype type_send;
            MPI_Type_indexed(Nrows_to_send_for_get_block(b, p),
                             block_lengths,
                             Rows_to_send_for_get_block(b, p),
                             MPI_DOUBLE,
                             &type_send);
            MPI_Type_commit(&type_send);
 
            // send
            MPI_Request send_req;
            MPI_Isend(const_cast<double*>(v.values_pt()),
                      1,
                      type_send,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &send_req);
            MPI_Type_free(&type_send);
            requests.push_back(send_req);
          }
 
          if (Nrows_to_recv_for_get_block(b, p) > 0)
          {
            // create the recv datatype
            MPI_Datatype type_recv;
            MPI_Type_indexed(Nrows_to_recv_for_get_block(b, p),
                             block_lengths,
                             Rows_to_recv_for_get_block(b, p),
                             MPI_DOUBLE,
                             &type_recv);
            MPI_Type_commit(&type_recv);
 
            // recv
            MPI_Request recv_req;
            MPI_Irecv(w.values_pt(),
                      1,
                      type_recv,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &recv_req);
            MPI_Type_free(&type_recv);
            requests.push_back(recv_req);
          }
        }
 
        // communicate with self
        else
        {
          double* w_values_pt = w.values_pt();
          const double* v_values_pt = v.values_pt();
          for (unsigned i = 0; i < Nrows_to_send_for_get_block(b, p); i++)
          {
            w_values_pt[Rows_to_recv_for_get_block(b, p)[i]] =
              v_values_pt[Rows_to_send_for_get_block(b, p)[i]];
          }
        }
      }
 
      // and then just wait
      unsigned c = requests.size();
      Vector<MPI_Status> stat(c);
      if (c)
      {
        MPI_Waitall(c, &requests[0], &stat[0]);
      }
      delete[] block_lengths;
 
#else
      // throw error
      std::ostringstream error_message;
      error_message << "The preconditioner is distributed and on more than one "
                    << "processor. MPI is required.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
#endif
    }
  }
 
  //============================================================================
  /// Takes the naturally ordered vector, v and returns the n-th
  /// block vector, b. Here n is the block number in the current
  /// preconditioner.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::get_block_vector(const unsigned& b,
                                                     const DoubleVector& v,
                                                     DoubleVector& w) const
  {
#ifdef PARANOID
    // the number of blocks
    const unsigned para_n_blocks = nblock_types();
 
    // paranoid check that block i is in this block preconditioner
    if (b >= para_n_blocks)
    {
      std::ostringstream err_msg;
      err_msg << "Requested block vector " << b
              << ", however this preconditioner has only " << para_n_blocks
              << " block types"
              << ".\n";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    if (!v.built())
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*(v.distribution_pt()) != *(this->master_distribution_pt()))
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must match the "
              << " specified master_distribution_pt(). \n"
              << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Recall that, the relationship between the external blocks and the
    // external dof types, as seen by the preconditioner writer is stored in the
    // mapping Block_to_dof_map_coarse.
    //
    // However, each dof type could have been coarsened! The relationship
    // between the dof types of this preconditioner and the parent
    // preconditioner is stored in the mapping Doftype_coarsen_map_coarse. The
    // dof numbers in this map is relative to this preconditioner.
    //
    // Finally, the relationship between the dof types of this preconditioner
    // and the most fine grain dof types is stored in the mapping
    // Doftype_coarsen_map_fine. Again, the dof numbers in this map is relative
    // to this preconditioner.
    //
    // Furthermore, we note that concatenation of vectors without communication
    //  is associative, but not commutative. I.e.
    // (V1+V2)+V3 = V1 + (V2 + V3), where + is concatenation without
    // communication.
    //
    // So all we need is the vectors listed in the correct order.
    //
    // We need only Block_to_dof_map_coarse to tell us which external dof types
    // are in this block, then Doftype_coarsen_map_fine to tell us which most
    // fine grain dofs to concatenate!
    //
    // All the mapping vectors are constructed to respect the ordering of
    // the dof types.
 
    // Get the most fine grain block to dof mapping.
    Vector<unsigned> most_fine_grain_dof = Block_to_dof_map_fine[b];
 
    // How many vectors do we need to concatenate?
    const unsigned n_dof_vec = most_fine_grain_dof.size();
 
    if (n_dof_vec == 1)
    // No need to concatenate, just extract the vector.
    {
      internal_get_block_vector(most_fine_grain_dof[0], v, w);
    }
    else
    // Need to concatenate dof-level vectors.
    {
      Vector<DoubleVector> dof_vector(n_dof_vec);
 
      // Get all the dof-level vectors in one go
      internal_get_block_vectors(most_fine_grain_dof, v, dof_vector);
      // Build w with the correct distribution.
      w.build(Block_distribution_pt[b], 0);
 
      // Concatenate the vectors.
      DoubleVectorHelpers::concatenate_without_communication(dof_vector, w);
 
      dof_vector.clear();
    }
  } // get_block_vector(...)
 
  //============================================================================
  /// Takes the n-th block ordered vector, b,  and copies its entries
  /// to the appropriate entries in the naturally ordered vector, v.
  /// Here n is the block number in the current block preconditioner.
  /// If the preconditioner is a subsidiary block preconditioner
  /// the other entries in v  that are not associated with it
  /// are left alone.
  ///
  /// This version works with the internal block types. This is legacy code
  /// but is kept alive, hence moved to private. Please use the
  /// function "return_block_vector(...)".
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::internal_return_block_vector(
    const unsigned& b, const DoubleVector& w, DoubleVector& v) const
  {
#ifdef PARANOID
    // the number of blocks
    const unsigned n_blocks = this->internal_nblock_types();
 
    // paranoid check that block i is in this block preconditioner
    if (b >= n_blocks)
    {
      std::ostringstream error_message;
      error_message
        << "Requested block  vector " << b
        << ", however this preconditioner has internal_nblock_types() "
        << "= " << internal_nblock_types() << std::endl;
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (!v.built())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*v.distribution_pt() != *this->master_distribution_pt())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must match the "
                    << " specified master_distribution_pt(). \n"
                    << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (!w.built())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the block vector w must be setup.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*w.distribution_pt() != *Internal_block_distribution_pt[b])
    {
      std::ostringstream error_message;
      error_message
        << "The distribution of the block vector w must match the "
        << " specified distribution at Internal_block_distribution_pt[b]";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // if + only one processor
    //    + more than one processor but matrix_pt is not distributed
    // then use the serial get_block method
    if (this->distribution_pt()->communicator_pt()->nproc() == 1 ||
        !this->distribution_pt()->distributed())
    {
      // length of vector
      unsigned n_row = this->internal_block_dimension(b);
 
      // copy back from the block vector to the naturally ordered vector
      double* v_pt = v.values_pt();
      const double* w_pt = w.values_pt();
      for (unsigned i = 0; i < n_row; i++)
      {
        v_pt[this->Global_index[b][i]] = w_pt[i];
      }
    }
    // otherwise use mpi
    else
    {
#ifdef OOMPH_HAS_MPI
 
      // my rank
      unsigned my_rank = this->distribution_pt()->communicator_pt()->my_rank();
 
      // the number of processors
      unsigned nproc = this->distribution_pt()->communicator_pt()->nproc();
 
      // determine the maximum number of rows to be sent or recv
      unsigned max_n_send_or_recv = 0;
      for (unsigned p = 0; p < nproc; p++)
      {
        max_n_send_or_recv =
          std::max(max_n_send_or_recv, Nrows_to_send_for_get_block(b, p));
        max_n_send_or_recv =
          std::max(max_n_send_or_recv, Nrows_to_recv_for_get_block(b, p));
      }
 
      // create a vectors of 1s (the size of the nblock for the mpi indexed
      // data types
      int* block_lengths = new int[max_n_send_or_recv];
      for (unsigned i = 0; i < max_n_send_or_recv; i++)
      {
        block_lengths[i] = 1;
      }
 
      // perform the sends and receives
      Vector<MPI_Request> requests;
      for (unsigned p = 0; p < nproc; p++)
      {
        // send and recv with other processors
        if (p != my_rank)
        {
          if (Nrows_to_recv_for_get_block(b, p) > 0)
          {
            // create the send datatype
            MPI_Datatype type_send;
            MPI_Type_indexed(Nrows_to_recv_for_get_block(b, p),
                             block_lengths,
                             Rows_to_recv_for_get_block(b, p),
                             MPI_DOUBLE,
                             &type_send);
            MPI_Type_commit(&type_send);
 
            // send
            MPI_Request send_req;
            MPI_Isend(const_cast<double*>(w.values_pt()),
                      1,
                      type_send,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &send_req);
            MPI_Type_free(&type_send);
            requests.push_back(send_req);
          }
 
          if (Nrows_to_send_for_get_block(b, p) > 0)
          {
            // create the recv datatype
            MPI_Datatype type_recv;
            MPI_Type_indexed(Nrows_to_send_for_get_block(b, p),
                             block_lengths,
                             Rows_to_send_for_get_block(b, p),
                             MPI_DOUBLE,
                             &type_recv);
            MPI_Type_commit(&type_recv);
 
            // recv
            MPI_Request recv_req;
            MPI_Irecv(v.values_pt(),
                      1,
                      type_recv,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &recv_req);
            MPI_Type_free(&type_recv);
            requests.push_back(recv_req);
          }
        }
 
        // communicate wih self
        else
        {
          const double* w_values_pt = w.values_pt();
          double* v_values_pt = v.values_pt();
          for (unsigned i = 0; i < Nrows_to_send_for_get_block(b, p); i++)
          {
            v_values_pt[Rows_to_send_for_get_block(b, p)[i]] =
              w_values_pt[Rows_to_recv_for_get_block(b, p)[i]];
          }
        }
      }
 
      // and then just wait
      unsigned c = requests.size();
      Vector<MPI_Status> stat(c);
      if (c)
      {
        MPI_Waitall(c, &requests[0], &stat[0]);
      }
      delete[] block_lengths;
 
#else
      // throw error
      std::ostringstream error_message;
      error_message << "The preconditioner is distributed and on more than one "
                    << "processor. MPI is required.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
#endif
    }
  }
 
  //============================================================================
  /// Takes the n-th block ordered vector, b,  and copies its entries
  /// to the appropriate entries in the naturally ordered vector, v.
  /// Here n is the block number in the current block preconditioner.
  /// If the preconditioner is a subsidiary block preconditioner
  /// the other entries in v  that are not associated with it
  /// are left alone.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::return_block_vector(const unsigned& n,
                                                        const DoubleVector& b,
                                                        DoubleVector& v) const
  {
#ifdef PARANOID
    // the number of blocks
    const unsigned para_n_blocks = nblock_types();
 
    // paranoid check that block i is in this block preconditioner
    if (n >= para_n_blocks)
    {
      std::ostringstream err_msg;
      err_msg << "Requested block vector " << b
              << ", however this preconditioner has " << para_n_blocks
              << " block types.\n";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (!v.built())
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*v.distribution_pt() != *this->master_distribution_pt())
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the global vector v must match the "
              << " specified master_distribution_pt(). \n"
              << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (!b.built())
    {
      std::ostringstream err_msg;
      err_msg << "The distribution of the block vector b must be setup.";
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
#endif
 
    // Get the most fine grain dof
    Vector<unsigned> most_fine_grain_dof = Block_to_dof_map_fine[n];
 
    // How many dofs are in this block?
    const unsigned n_dof_vec = Block_to_dof_map_fine[n].size();
 
    if (n_dof_vec == 1)
    // There is only one dof, no need to split.
    {
      internal_return_block_vector(most_fine_grain_dof[0], b, v);
    }
    else
    // Need to split the vector up before we insert them all in one go.
    {
      Vector<DoubleVector> dof_vector(n_dof_vec);
      for (unsigned d = 0; d < n_dof_vec; d++)
      {
        dof_vector[d].build(
          internal_block_distribution_pt(most_fine_grain_dof[d]));
      }
 
      DoubleVectorHelpers::split_without_communication(b, dof_vector);
 
      // return to v
      internal_return_block_vectors(most_fine_grain_dof, dof_vector, v);
    }
  } // return_block_vector(...)
 
  //============================================================================
  /// Given the naturally ordered vector, v, return
  /// the vector rearranged in block order in w. This is a legacy function
  /// from the old block preconditioning framework. Kept alive in case it may
  /// be needed again.
  ///
  /// This uses the variables ending in "get_ordered". We no longer use this
  /// type of method. This function copy values from v and re-order them
  /// in "block order" and place them in w. Block order means that the
  /// values in w are the same as the concatenated block vectors.
  ///
  /// I.e. - v is naturally ordered.
  ///        v -> s_b, v is ordered into blocks vectors
  ///                  (requires communication)
  ///        concatenate_without_communication(s_{0,...,nblocks},w) gives w.
  ///
  /// But this function skips out the concatenation part and builds w directly
  /// from v.
  ///
  /// This is nice but the function is implemented in such a way that it
  /// always use all the (internal) blocks and concatenated with the
  /// identity ordering. I.e. if this preconditioner has 3 block types, then
  /// w will always be:
  /// concatenate_without_communication([s_0, s_1, s_2], w). There is easy
  /// way to change this.
  ///
  /// Furthermore, it does not take into account the new dof type coarsening
  /// feature. So this function will most likely produce the incorrect vector
  /// w from what the user intended. It still works, but w will be the
  /// concatenation of the most fine grain dof block vectors with the
  /// "natural" dof type ordering.
  ///
  /// This has been superseded by the function
  /// get_block_ordered_preconditioner_vector(...) which does the correct
  /// thing.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::
    internal_get_block_ordered_preconditioner_vector(const DoubleVector& v,
                                                     DoubleVector& w) const
  {
#ifdef PARANOID
    if (!v.built())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*v.distribution_pt() != *this->master_distribution_pt())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must match the "
                    << " specified master_distribution_pt(). \n"
                    << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    //  Cleared and resized w for reordered vector
    w.build(this->internal_preconditioner_matrix_distribution_pt(), 0.0);
 
    // if + only one processor
    //    + more than one processor but matrix_pt is not distributed
    // then use the serial get_block method
    if (this->distribution_pt()->communicator_pt()->nproc() == 1 ||
        !this->distribution_pt()->distributed())
    {
      // number of blocks
      unsigned nblock = this->Internal_nblock_types;
 
      // copy to w
      unsigned block_offset = 0;
      double* w_pt = w.values_pt();
      const double* v_pt = v.values_pt();
      for (unsigned b = 0; b < nblock; b++)
      {
        unsigned block_nrow = this->internal_block_dimension(b);
        for (unsigned i = 0; i < block_nrow; i++)
        {
          w_pt[block_offset + i] = v_pt[this->Global_index[b][i]];
        }
        block_offset += block_nrow;
      }
    }
    // otherwise use mpi
    else
    {
#ifdef OOMPH_HAS_MPI
 
      // my rank
      unsigned my_rank = this->distribution_pt()->communicator_pt()->my_rank();
 
      // the number of processors
      unsigned nproc = this->distribution_pt()->communicator_pt()->nproc();
 
      // determine the maximum number of rows to be sent or recv
      unsigned max_n_send_or_recv = 0;
      for (unsigned p = 0; p < nproc; p++)
      {
        max_n_send_or_recv =
          std::max(max_n_send_or_recv, Nrows_to_send_for_get_ordered[p]);
        max_n_send_or_recv =
          std::max(max_n_send_or_recv, Nrows_to_recv_for_get_ordered[p]);
      }
 
      // create a vectors of 1s (the size of the nblock for the mpi indexed
      // data types
      int* block_lengths = new int[max_n_send_or_recv];
      for (unsigned i = 0; i < max_n_send_or_recv; i++)
      {
        block_lengths[i] = 1;
      }
 
      // perform the sends and receives
      Vector<MPI_Request> requests;
      for (unsigned p = 0; p < nproc; p++)
      {
        // send and recv with other processors
        if (p != my_rank)
        {
          if (Nrows_to_send_for_get_ordered[p] > 0)
          {
            // create the send datatype
            MPI_Datatype type_send;
            MPI_Type_indexed(Nrows_to_send_for_get_ordered[p],
                             block_lengths,
                             Rows_to_send_for_get_ordered[p],
                             MPI_DOUBLE,
                             &type_send);
            MPI_Type_commit(&type_send);
 
            // send
            MPI_Request send_req;
            MPI_Isend(const_cast<double*>(v.values_pt()),
                      1,
                      type_send,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &send_req);
            MPI_Type_free(&type_send);
            requests.push_back(send_req);
          }
 
          if (Nrows_to_recv_for_get_ordered[p] > 0)
          {
            // create the recv datatype
            MPI_Datatype type_recv;
            MPI_Type_indexed(Nrows_to_recv_for_get_ordered[p],
                             block_lengths,
                             Rows_to_recv_for_get_ordered[p],
                             MPI_DOUBLE,
                             &type_recv);
            MPI_Type_commit(&type_recv);
 
            // recv
            MPI_Request recv_req;
            MPI_Irecv(w.values_pt(),
                      1,
                      type_recv,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &recv_req);
            MPI_Type_free(&type_recv);
            requests.push_back(recv_req);
          }
        }
 
        // communicate with self
        else
        {
          double* w_values_pt = w.values_pt();
          const double* v_values_pt = v.values_pt();
          for (unsigned i = 0; i < Nrows_to_send_for_get_ordered[p]; i++)
          {
            w_values_pt[Rows_to_recv_for_get_ordered[p][i]] =
              v_values_pt[Rows_to_send_for_get_ordered[p][i]];
          }
        }
      }
 
      // and then just wait
      unsigned c = requests.size();
      Vector<MPI_Status> stat(c);
      if (c)
      {
        MPI_Waitall(c, &requests[0], &stat[0]);
      }
      delete[] block_lengths;
 
#else
      // throw error
      std::ostringstream error_message;
      error_message << "The preconditioner is distributed and on more than one "
                    << "processor. MPI is required.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
#endif
    }
  }
 
  //============================================================================
  /// Given the naturally ordered vector, v, return
  /// the vector rearranged in block order in w. This function calls
  /// get_concatenated_block_vector(...) with the identity block mapping.
  ///
  /// This function has been re-written to work with the new dof type
  /// coarsening feature. The old function is kept alive in
  /// internal_get_block_ordered_preconditioner_vector(...) and is moved to
  /// the private section of the code. The differences between the two are:
  ///
  /// 1) This function extracts all the block vectors (in one go) via the
  ///    function internal_get_block_vectors(...), and concatenates them.
  ///
  /// 2) The old function makes use of the variables ending in "get_ordered",
  ///    thus is slightly more efficient since it does not have to concatenate
  ///    any block vectors.
  ///
  /// 3) The old function no longer respect the new indirections if dof types
  ///    have been coarsened.
  ///
  /// 4) This function extracts the most fine grain dof-level vectors and
  ///    concatenates them. These dof-level vectors respect the re-ordering
  ///    caused by the coarsening of dof types. The overhead associated with
  ///    concatenating DoubleVectors without communication is very small.
  ///
  /// This function should be used.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::get_block_ordered_preconditioner_vector(
    const DoubleVector& v, DoubleVector& w)
  {
#ifdef PARANOID
    if (!v.built())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*v.distribution_pt() != *this->master_distribution_pt())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must match the "
                    << " specified master_distribution_pt(). \n"
                    << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Get the number of blocks.
    unsigned nblocks = this->nblock_types();
 
    // Fill in the identity mapping.
    Vector<unsigned> block_vec_number(nblocks, 0);
    for (unsigned b = 0; b < nblocks; b++)
    {
      block_vec_number[b] = b;
    }
 
    // Do the work.
    get_concatenated_block_vector(block_vec_number, v, w);
  } // get_block_ordered_preconditioner_vector(...)
 
  //============================================================================
  /// Takes the block ordered vector, w, and reorders it in the natural
  /// order. Reordered vector is returned in v. Note: If the preconditioner is
  /// a subsidiary preconditioner then only the components of the vector
  /// associated with the blocks of the subsidiary preconditioner will be
  /// included. Hence the length of v is master_nrow() whereas that of the
  /// vector w is of length this->nrow().
  ///
  /// This is the return function for the function
  /// internal_get_block_ordered_preconditioner_vector(...).
  /// Both internal_get_block_ordered_preconditioner_vector(...) and
  /// internal_return_block_ordered_preconditioner_vector(...) has been
  /// superseded by the functions
  ///
  /// get_block_ordered_preconditioner_vector(...) and
  /// return_block_ordered_preconditioner_vector(...),
  ///
  /// Thus this function is moved to the private section of the code.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::
    internal_return_block_ordered_preconditioner_vector(const DoubleVector& w,
                                                        DoubleVector& v) const
  {
#ifdef PARANOID
    if (!v.built())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*v.distribution_pt() != *this->master_distribution_pt())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must match the "
                    << " specified master_distribution_pt(). \n"
                    << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (!w.built())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the block vector w must be setup.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*w.distribution_pt() !=
        *this->internal_preconditioner_matrix_distribution_pt())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the block vector w must match the "
                    << " specified distribution at Distribution_pt[b]";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
 
    // if + only one processor
    //    + more than one processor but matrix_pt is not distributed
    // then use the serial get_block method
    if (this->distribution_pt()->communicator_pt()->nproc() == 1 ||
        !this->distribution_pt()->distributed())
    {
      // number of blocks
      unsigned nblock = this->Internal_nblock_types;
 
      // copy to w
      unsigned block_offset = 0;
      const double* w_pt = w.values_pt();
      double* v_pt = v.values_pt();
      for (unsigned b = 0; b < nblock; b++)
      {
        unsigned block_nrow = this->internal_block_dimension(b);
        for (unsigned i = 0; i < block_nrow; i++)
        {
          v_pt[this->Global_index[b][i]] = w_pt[block_offset + i];
        }
        block_offset += block_nrow;
      }
    }
    // otherwise use mpi
    else
    {
#ifdef OOMPH_HAS_MPI
 
      // my rank
      unsigned my_rank = this->distribution_pt()->communicator_pt()->my_rank();
 
      // the number of processors
      unsigned nproc = this->distribution_pt()->communicator_pt()->nproc();
 
      // determine the maximum number of rows to be sent or recv
      unsigned max_n_send_or_recv = 0;
      for (unsigned p = 0; p < nproc; p++)
      {
        max_n_send_or_recv =
          std::max(max_n_send_or_recv, Nrows_to_send_for_get_ordered[p]);
        max_n_send_or_recv =
          std::max(max_n_send_or_recv, Nrows_to_recv_for_get_ordered[p]);
      }
 
      // create a vectors of 1s (the size of the nblock for the mpi indexed
      // data types
      int* block_lengths = new int[max_n_send_or_recv];
      for (unsigned i = 0; i < max_n_send_or_recv; i++)
      {
        block_lengths[i] = 1;
      }
 
      // perform the sends and receives
      Vector<MPI_Request> requests;
      for (unsigned p = 0; p < nproc; p++)
      {
        // send and recv with other processors
        if (p != my_rank)
        {
          if (Nrows_to_recv_for_get_ordered[p] > 0)
          {
            // create the send datatype
            MPI_Datatype type_send;
            MPI_Type_indexed(Nrows_to_recv_for_get_ordered[p],
                             block_lengths,
                             Rows_to_recv_for_get_ordered[p],
                             MPI_DOUBLE,
                             &type_send);
            MPI_Type_commit(&type_send);
 
            // send
            MPI_Request send_req;
            MPI_Isend(const_cast<double*>(w.values_pt()),
                      1,
                      type_send,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &send_req);
            MPI_Type_free(&type_send);
            requests.push_back(send_req);
          }
 
          if (Nrows_to_send_for_get_ordered[p] > 0)
          {
            // create the recv datatype
            MPI_Datatype type_recv;
            MPI_Type_indexed(Nrows_to_send_for_get_ordered[p],
                             block_lengths,
                             Rows_to_send_for_get_ordered[p],
                             MPI_DOUBLE,
                             &type_recv);
            MPI_Type_commit(&type_recv);
 
            // recv
            MPI_Request recv_req;
            MPI_Irecv(v.values_pt(),
                      1,
                      type_recv,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &recv_req);
            MPI_Type_free(&type_recv);
            requests.push_back(recv_req);
          }
        }
 
        // communicate wih self
        else
        {
          const double* w_values_pt = w.values_pt();
          double* v_values_pt = v.values_pt();
          for (unsigned i = 0; i < Nrows_to_send_for_get_ordered[p]; i++)
          {
            v_values_pt[Rows_to_send_for_get_ordered[p][i]] =
              w_values_pt[Rows_to_recv_for_get_ordered[p][i]];
          }
        }
      }
 
      // and then just wait
      unsigned c = requests.size();
      Vector<MPI_Status> stat(c);
      if (c)
      {
        MPI_Waitall(c, &requests[0], &stat[0]);
      }
      delete[] block_lengths;
 
#else
      // throw error
      std::ostringstream error_message;
      error_message << "The preconditioner is distributed and on more than one "
                    << "processor. MPI is required.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
#endif
    } // else use mpi
  } // function return_block_ordered_preconditioner_vector
 
 
  //============================================================================
  /// Takes the block ordered vector, w, and reorders it in natural
  /// order. Reordered vector is returned in v. Note: If the preconditioner is
  /// a subsidiary preconditioner then only the components of the vector
  /// associated with the blocks of the subsidiary preconditioner will be
  /// included. Hence the length of v is master_nrow() whereas that of the
  /// vector w is of length this->nrow().
  ///
  /// This is the return function for the function
  /// get_block_ordered_preconditioner_vector(...).
  ///
  /// It calls the function return_concatenated_block_vector(...) with the
  /// identity block number ordering.
  //============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::return_block_ordered_preconditioner_vector(
    const DoubleVector& w, DoubleVector& v) const
  {
#ifdef PARANOID
    if (!v.built())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must be setup.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*v.distribution_pt() != *this->master_distribution_pt())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the global vector v must match the "
                    << " specified master_distribution_pt(). \n"
                    << "i.e. Distribution_pt in the master preconditioner";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (!w.built())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the block vector w must be setup.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
    if (*w.distribution_pt() != *this->preconditioner_matrix_distribution_pt())
    {
      std::ostringstream error_message;
      error_message << "The distribution of the block vector w must match the "
                    << "concatenations of distributions in "
                    << "Block_distribution_pt.\n";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Split into the block vectors.
    const unsigned nblocks = nblock_types();
    Vector<unsigned> block_vec_number(nblocks, 0);
    for (unsigned b = 0; b < nblocks; b++)
    {
      block_vec_number[b] = b;
    }
 
    return_concatenated_block_vector(block_vec_number, w, v);
  } // function return_block_ordered_preconditioner_vector
 
  //=============================================================================
  /// Gets block (i,j) from the matrix pointed to by
  /// Matrix_pt and returns it in output_block. This is associated with the
  /// internal blocks. Please use the other get_block(...) function.
  //=============================================================================
  template<>
  void BlockPreconditioner<CRDoubleMatrix>::internal_get_block(
    const unsigned& block_i,
    const unsigned& block_j,
    CRDoubleMatrix& output_block) const
  {
#ifdef PARANOID
    // the number of blocks
    const unsigned n_blocks = this->internal_nblock_types();
 
    // paranoid check that block i is in this block preconditioner
    if (block_i >= n_blocks || block_j >= n_blocks)
    {
      std::ostringstream error_message;
      error_message
        << "Requested block (" << block_i << "," << block_j
        << "), however this preconditioner has internal_nblock_types() "
        << "= " << internal_nblock_types() << std::endl;
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Check that the matrix is the same as that of the master
    if (is_subsidiary_block_preconditioner())
    {
      if (master_block_preconditioner_pt()->matrix_pt() != matrix_pt())
      {
        std::string err = "Master and subs should have same matrix.";
        throw OomphLibError(
          err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
#endif
 
    // Cast the pointer
    CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());
 
    // if + only one processor
    //    + more than one processor but matrix_pt is not distributed
    // then use the serial get_block method
    if (cr_matrix_pt->distribution_pt()->communicator_pt()->nproc() == 1 ||
        !cr_matrix_pt->distribution_pt()->distributed())
    {
      // pointers for the jacobian matrix is compressed row sparse format
      int* j_row_start;
      int* j_column_index;
      double* j_value;
 
      // sets pointers to jacobian matrix
      j_row_start = cr_matrix_pt->row_start();
      j_column_index = cr_matrix_pt->column_index();
      j_value = cr_matrix_pt->value();
 
      // get the block dimensions
      unsigned block_nrow = this->internal_block_dimension(block_i);
      unsigned block_ncol = this->internal_block_dimension(block_j);
 
      // allocate temporary storage for the component vectors of block (i,j)
      // temp_ptr is used to point to an element in each column - required as
      // cannot assume that order of block's rows in jacobian and the block
      // matrix will be the same
      int* temp_row_start = new int[block_nrow + 1];
      for (unsigned i = 0; i <= block_nrow; i++)
      {
        temp_row_start[i] = 0;
      }
      Vector<int> temp_ptr(block_nrow + 1);
      int block_nnz = 0;
 
      // get number of rows in source matrix
      unsigned master_nrow = this->master_nrow();
 
      // determine how many non zeros there are in the block (i,j)
      // also determines how many non zeros are stored in each row or column -
      // stored in temp_ptr temporarily
      for (unsigned k = 0; k < master_nrow; k++)
      {
        if (internal_block_number(k) == static_cast<int>(block_i))
        {
          for (int l = j_row_start[k]; l < j_row_start[k + 1]; l++)
          {
            if (internal_block_number(j_column_index[l]) ==
                static_cast<int>(block_j))
            {
              block_nnz++;
              temp_ptr[internal_index_in_block(k) + 1]++;
            }
          }
        }
      }
 
      // if the matrix is not empty
      int* temp_column_index = new int[block_nnz];
      double* temp_value = new double[block_nnz];
      if (block_nnz > 0)
      {
        // uses number of elements in each column of block to determine values
        // for the block column start (temp_row_start)
        temp_row_start[0] = 0;
        for (unsigned k = 1; k <= block_nrow; k++)
        {
          temp_row_start[k] = temp_row_start[k - 1] + temp_ptr[k];
          temp_ptr[k] = temp_row_start[k];
        }
 
        // copies the relevant elements of the jacobian to the correct entries
        // of the block matrix
        for (unsigned k = 0; k < master_nrow; k++)
        {
          if (internal_block_number(k) == static_cast<int>(block_i))
          {
            for (int l = j_row_start[k]; l < j_row_start[k + 1]; l++)
            {
              if (internal_block_number(j_column_index[l]) ==
                  static_cast<int>(block_j))
              {
                int kk = temp_ptr[internal_index_in_block(k)]++;
                temp_value[kk] = j_value[l];
                temp_column_index[kk] =
                  internal_index_in_block(j_column_index[l]);
              }
            }
          }
        }
      }
 
 
      // Fill in the compressed row matrix ??ds Note: I kept the calls to
      // build as close as I could to before (had to replace new(dist) with
      // .build(dist) ).
      output_block.build(Internal_block_distribution_pt[block_i]);
      output_block.build_without_copy(
        block_ncol, block_nnz, temp_value, temp_column_index, temp_row_start);
 
#ifdef PARANOID
      // checks to see if block matrix has been set up correctly
      //   block_matrix_test(matrix_pt,block_i,block_j,block_pt);
      if (Run_block_matrix_test)
      {
        // checks to see if block matrix has been set up correctly
        block_matrix_test(block_i, block_j, &output_block);
      }
#endif
    }
 
 
    // otherwise we are dealing with a distributed matrix
    else
    {
#ifdef OOMPH_HAS_MPI
      // number of processors
      unsigned nproc = this->distribution_pt()->communicator_pt()->nproc();
 
      // my rank
      unsigned my_rank = this->distribution_pt()->communicator_pt()->my_rank();
 
      // sets pointers to jacobian matrix
      int* j_row_start = cr_matrix_pt->row_start();
      int* j_column_index = cr_matrix_pt->column_index();
      double* j_value = cr_matrix_pt->value();
 
      // number of non zeros in each row to be sent
      Vector<int*> nnz_send(nproc, 0);
 
      // number of non zeros in each row to be received
      Vector<int*> nnz_recv(nproc, 0);
 
      // storage for data to be sent
      Vector<int*> column_index_for_proc(nproc, 0);
      Vector<double*> values_for_proc(nproc, 0);
 
      // number of non zeros to be sent to each processor
      Vector<unsigned> total_nnz_send(nproc, 0);
 
      // number of rows of the block matrix on this processor
      unsigned nrow_local =
        Internal_block_distribution_pt[block_i]->nrow_local();
 
      // resize the nnz storage and compute nnz_send
      // and send and recv the nnz
      Vector<MPI_Request> send_req;
      Vector<MPI_Request> recv1_req;
      for (unsigned p = 0; p < nproc; p++)
      {
        int nrow_send = Nrows_to_send_for_get_block(block_i, p);
        int nrow_recv = Nrows_to_recv_for_get_block(block_i, p);
 
        // assemble nnz recv
        nnz_recv[p] = new int[nrow_recv];
 
        // assemble the storage to send
        if (nrow_send > 0 && p != my_rank)
        {
          nnz_send[p] = new int[nrow_send];
        }
 
        // compute the number of nnzs in each row and the total number
        // of nnzs
        for (int i = 0; i < nrow_send; i++)
        {
          unsigned row = Rows_to_send_for_get_block(block_i, p)[i];
          int c = 0;
          for (int r = j_row_start[row]; r < j_row_start[row + 1]; r++)
          {
            if (internal_block_number(j_column_index[r]) == int(block_j))
            {
              c++;
            }
          }
          if (p != my_rank)
          {
            nnz_send[p][i] = c;
          }
          else
          {
            nnz_recv[p][i] = c;
          }
          total_nnz_send[p] += c;
        }
 
        // send
        if (p != my_rank)
        {
          if (nrow_send)
          {
            MPI_Request req;
            MPI_Isend(nnz_send[p],
                      nrow_send,
                      MPI_INT,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &req);
            send_req.push_back(req);
          }
 
          // recv
          if (nrow_recv)
          {
            MPI_Request req;
            MPI_Irecv(nnz_recv[p],
                      nrow_recv,
                      MPI_INT,
                      p,
                      0,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &req);
            recv1_req.push_back(req);
          }
        }
      }
 
      // next assemble the values and row_start data to be sent for each
      // processor
      for (unsigned p = 0; p < nproc; p++)
      {
        int nrow_send = Nrows_to_send_for_get_block(block_i, p);
 
        // assemble the storage for the values and column indices to be sent
        if (p != my_rank)
        {
          if (total_nnz_send[p] > 0)
          {
            values_for_proc[p] = new double[total_nnz_send[p]];
            column_index_for_proc[p] = new int[total_nnz_send[p]];
 
            // copy the values and column indices to the storage
            unsigned ptr = 0;
            for (int i = 0; i < nrow_send; i++)
            {
              unsigned row = Rows_to_send_for_get_block(block_i, p)[i];
              for (int r = j_row_start[row]; r < j_row_start[row + 1]; r++)
              {
                if (internal_block_number(j_column_index[r]) == int(block_j))
                {
                  values_for_proc[p][ptr] = j_value[r];
                  column_index_for_proc[p][ptr] =
                    internal_index_in_block(j_column_index[r]);
                  ptr++;
                }
              }
            }
 
            // create the datatypes
            MPI_Datatype types[2];
            MPI_Type_contiguous(total_nnz_send[p], MPI_DOUBLE, &types[0]);
            MPI_Type_commit(&types[0]);
            MPI_Type_contiguous(total_nnz_send[p], MPI_INT, &types[1]);
            MPI_Type_commit(&types[1]);
 
            // get the start address of the vectors
            MPI_Aint displacement[2];
            MPI_Get_address(values_for_proc[p], &displacement[0]);
            MPI_Get_address(column_index_for_proc[p], &displacement[1]);
 
            // compute the displacements
            displacement[1] -= displacement[0];
            displacement[0] -= displacement[0];
 
            // compute the block lengths
            int length[2];
            length[0] = length[1] = 1;
 
            // build the struct data type
            MPI_Datatype final_type;
            MPI_Type_create_struct(2, length, displacement, types, &final_type);
            MPI_Type_commit(&final_type);
            MPI_Type_free(&types[0]);
            MPI_Type_free(&types[1]);
 
            // and send
            MPI_Request req;
            MPI_Isend(values_for_proc[p],
                      1,
                      final_type,
                      p,
                      1,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &req);
            send_req.push_back(req);
            MPI_Type_free(&final_type);
          }
        }
      }
 
      // wait for the recv to complete (the row_start recv which actually
      // contains the number of nnzs in each row)
      int c_recv = recv1_req.size();
      if (c_recv != 0)
      {
        MPI_Waitall(c_recv, &recv1_req[0], MPI_STATUS_IGNORE);
      }
 
      // compute the total number of nnzs to be received
      Vector<int> total_nnz_recv_from_proc(nproc);
      int local_block_nnz = 0;
      for (unsigned p = 0; p < nproc; p++)
      {
        // compute the total nnzs
        for (unsigned i = 0; i < Nrows_to_recv_for_get_block(block_i, p); i++)
        {
          total_nnz_recv_from_proc[p] += nnz_recv[p][i];
        }
        local_block_nnz += total_nnz_recv_from_proc[p];
      }
 
      // compute the offset for each block of nnzs (a matrix row) in the
      // values_recv and column_index_recv vectors
 
      // fisrt determine how many blocks of rows are to be recv
      Vector<int> n_recv_block(nproc, 0);
      for (unsigned p = 0; p < nproc; p++)
      {
        if (Nrows_to_recv_for_get_block(block_i, p) > 0)
        {
          n_recv_block[p] = 1;
        }
        for (unsigned i = 1; i < Nrows_to_recv_for_get_block(block_i, p); i++)
        {
          if (Rows_to_recv_for_get_block(block_i, p)[i] !=
              Rows_to_recv_for_get_block(block_i, p)[i - 1] + 1)
          {
            n_recv_block[p]++;
          }
        }
      }
 
      // next assemble row start recv
      int* row_start_recv = new int[nrow_local + 1];
      for (unsigned i = 0; i <= nrow_local; i++)
      {
        row_start_recv[i] = 0;
      }
      for (unsigned p = 0; p < nproc; p++)
      {
        for (unsigned i = 0; i < Nrows_to_recv_for_get_block(block_i, p); i++)
        {
          row_start_recv[Rows_to_recv_for_get_block(block_i, p)[i]] =
            nnz_recv[p][i];
        }
      }
      int g = row_start_recv[0];
      row_start_recv[0] = 0;
      for (unsigned i = 1; i < nrow_local; i++)
      {
        int temp_g = g;
        g = row_start_recv[i];
        row_start_recv[i] = row_start_recv[i - 1] + temp_g;
      }
      row_start_recv[nrow_local] = row_start_recv[nrow_local - 1] + g;
 
      // next assemble the offset and the number of nzs in each recv block
      Vector<int*> offset_recv_block(nproc, 0);
      Vector<int*> nnz_recv_block(nproc, 0);
      for (unsigned p = 0; p < nproc; p++)
      {
        if (Nrows_to_recv_for_get_block(block_i, p) > 0)
        {
          offset_recv_block[p] = new int[n_recv_block[p]];
          offset_recv_block[p][0] = 0;
          nnz_recv_block[p] = new int[n_recv_block[p]];
          for (int i = 0; i < n_recv_block[p]; i++)
          {
            nnz_recv_block[p][i] = 0;
          }
          unsigned ptr = 0;
          nnz_recv_block[p][ptr] += nnz_recv[p][0];
          offset_recv_block[p][0] =
            row_start_recv[Rows_to_recv_for_get_block(block_i, p)[0]];
          for (unsigned i = 1; i < Nrows_to_recv_for_get_block(block_i, p); i++)
          {
            if (Rows_to_recv_for_get_block(block_i, p)[i] !=
                Rows_to_recv_for_get_block(block_i, p)[i - 1] + 1)
            {
              ptr++;
              offset_recv_block[p][ptr] =
                row_start_recv[Rows_to_recv_for_get_block(block_i, p)[i]];
            }
            nnz_recv_block[p][ptr] += nnz_recv[p][i];
          }
        }
        delete[] nnz_recv[p];
      }
 
      // post the receives
      int* column_index_recv = new int[local_block_nnz];
      double* values_recv = new double[local_block_nnz];
      Vector<MPI_Request> recv2_req;
      for (unsigned p = 0; p < nproc; p++)
      {
        if (p != my_rank)
        {
          if (total_nnz_recv_from_proc[p] != 0)
          {
            // create the datatypes
            MPI_Datatype types[2];
            MPI_Type_indexed(n_recv_block[p],
                             nnz_recv_block[p],
                             offset_recv_block[p],
                             MPI_DOUBLE,
                             &types[0]);
            MPI_Type_commit(&types[0]);
            MPI_Type_indexed(n_recv_block[p],
                             nnz_recv_block[p],
                             offset_recv_block[p],
                             MPI_INT,
                             &types[1]);
            MPI_Type_commit(&types[1]);
 
            // compute the displacements
            MPI_Aint displacements[2];
            MPI_Get_address(values_recv, &displacements[0]);
            MPI_Get_address(column_index_recv, &displacements[1]);
            displacements[1] -= displacements[0];
            displacements[0] -= displacements[0];
 
            // compute the block lengths
            int length[2];
            length[0] = length[1] = 1;
 
            // create the final datatype
            MPI_Datatype final_type;
            MPI_Type_create_struct(
              2, length, displacements, types, &final_type);
            MPI_Type_commit(&final_type);
            MPI_Type_free(&types[0]);
            MPI_Type_free(&types[1]);
 
            // and the recv
            MPI_Request req;
            MPI_Irecv(values_recv,
                      1,
                      final_type,
                      p,
                      1,
                      this->distribution_pt()->communicator_pt()->mpi_comm(),
                      &req);
            recv2_req.push_back(req);
            MPI_Type_free(&final_type);
          }
        }
        else
        {
          // next send the values and column indices to self
          unsigned block_ptr = 0;
          unsigned counter = 0;
          int nrow_send = Nrows_to_send_for_get_block(block_i, my_rank);
          if (nrow_send > 0)
          {
            unsigned offset = offset_recv_block[my_rank][0];
            for (int i = 0; i < nrow_send; i++)
            {
              if (i > 0)
              {
                if (Rows_to_recv_for_get_block(block_i, p)[i] !=
                    Rows_to_recv_for_get_block(block_i, p)[i - 1] + 1)
                {
                  counter = 0;
                  block_ptr++;
                  offset = offset_recv_block[my_rank][block_ptr];
                }
              }
              unsigned row = Rows_to_send_for_get_block(block_i, my_rank)[i];
              for (int r = j_row_start[row]; r < j_row_start[row + 1]; r++)
              {
                if (internal_block_number(j_column_index[r]) == int(block_j))
                {
                  values_recv[offset + counter] = j_value[r];
                  column_index_recv[offset + counter] =
                    internal_index_in_block(j_column_index[r]);
                  counter++;
                }
              }
            }
          }
        }
      }
 
      // wait for the recv to complete (for the column_index and the values_
      c_recv = recv2_req.size();
      if (c_recv != 0)
      {
        MPI_Waitall(c_recv, &recv2_req[0], MPI_STATUS_IGNORE);
      }
 
      // Fill in the compressed row matrix
      output_block.build(Internal_block_distribution_pt[block_i]);
      output_block.build_without_copy(this->internal_block_dimension(block_j),
                                      local_block_nnz,
                                      values_recv,
                                      column_index_recv,
                                      row_start_recv);
 
      // wait for the send to complete (nnz / row_start)
      int c_send = send_req.size();
      if (c_send)
      {
        MPI_Waitall(c_send, &send_req[0], MPI_STATUS_IGNORE);
      }
 
      // delete temp storage used for assembling data for communication
      for (unsigned p = 0; p < nproc; p++)
      {
        delete[] nnz_send[p];
        delete[] column_index_for_proc[p];
        delete[] values_for_proc[p];
        delete[] offset_recv_block[p];
        delete[] nnz_recv_block[p];
      }
#else
      // throw error
      std::ostringstream error_message;
      error_message << "The matrix is distributed and on more than one "
                    << "processor. MPI is required.";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
#endif
    }
  }
 
  //=============================================================================
  /// Gets dof-level block (i,j).
  /// If Replacement_dof_block_pt(i,j) is not null, then the replacement
  /// block is returned via a deep copy.
  ///
  /// Otherwise if this is the uppermost block preconditioner then it calls
  /// internal_get_block(i,j), else if it is a subsidiary
  /// block preconditioner, it will call it's master block preconditioners'
  /// get_dof_level_block function.
  //=============================================================================
  template<>
  void BlockPreconditioner<CRDoubleMatrix>::get_dof_level_block(
    const unsigned& block_i,
    const unsigned& block_j,
    CRDoubleMatrix& output_block,
    const bool& ignore_replacement_block) const
  {
#ifdef PARANOID
    // the number of dof types.
    unsigned para_ndofs = ndof_types();
 
    // paranoid check that block i is in this block preconditioner
    if (block_i >= para_ndofs || block_j >= para_ndofs)
    {
      std::ostringstream err_msg;
      err_msg << "Requested dof block (" << block_i << "," << block_j
              << "), however this preconditioner has ndof_types() "
              << "= " << para_ndofs << std::endl;
      throw OomphLibError(
        err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    CRDoubleMatrix* tmp_block_pt =
      Replacement_dof_block_pt.get(block_i, block_j);
 
    if ((tmp_block_pt == 0) || ignore_replacement_block)
    {
      // Getting the block from parent preconditioner
      const unsigned ndof_in_parent_i =
        Doftype_coarsen_map_coarse[block_i].size();
      const unsigned ndof_in_parent_j =
        Doftype_coarsen_map_coarse[block_j].size();
 
      if (ndof_in_parent_i == 1 && ndof_in_parent_j == 1)
      {
        unsigned parent_dof_i = Doftype_coarsen_map_coarse[block_i][0];
        unsigned parent_dof_j = Doftype_coarsen_map_coarse[block_j][0];
 
        if (is_master_block_preconditioner())
        {
          internal_get_block(parent_dof_i, parent_dof_j, output_block);
        }
        else
        {
          parent_dof_i = Doftype_in_master_preconditioner_coarse[parent_dof_i];
          parent_dof_j = Doftype_in_master_preconditioner_coarse[parent_dof_j];
 
          master_block_preconditioner_pt()->get_dof_level_block(
            parent_dof_i, parent_dof_j, output_block, ignore_replacement_block);
        }
      }
      else
      {
        DenseMatrix<CRDoubleMatrix*> tmp_blocks_pt(
          ndof_in_parent_i, ndof_in_parent_j, 0);
 
        Vector<Vector<unsigned>> new_block(
          ndof_in_parent_i, Vector<unsigned>(ndof_in_parent_j, 0));
 
        for (unsigned dof_i = 0; dof_i < ndof_in_parent_i; dof_i++)
        {
          unsigned parent_dof_i = Doftype_coarsen_map_coarse[block_i][dof_i];
          if (is_subsidiary_block_preconditioner())
          {
            parent_dof_i =
              Doftype_in_master_preconditioner_coarse[parent_dof_i];
          }
 
          for (unsigned dof_j = 0; dof_j < ndof_in_parent_j; dof_j++)
          {
            unsigned parent_dof_j = Doftype_coarsen_map_coarse[block_j][dof_j];
 
            tmp_blocks_pt(dof_i, dof_j) = new CRDoubleMatrix;
 
            new_block[dof_i][dof_j] = 1;
 
            if (is_master_block_preconditioner())
            {
              internal_get_block(
                parent_dof_i, parent_dof_j, *tmp_blocks_pt(dof_i, dof_j));
            }
            else
            {
              parent_dof_j =
                Doftype_in_master_preconditioner_coarse[parent_dof_j];
 
              master_block_preconditioner_pt()->get_dof_level_block(
                parent_dof_i,
                parent_dof_j,
                *tmp_blocks_pt(dof_i, dof_j),
                ignore_replacement_block);
            }
          }
        }
 
        Vector<LinearAlgebraDistribution*> tmp_row_dist_pt(ndof_in_parent_i, 0);
 
        for (unsigned parent_dof_i = 0; parent_dof_i < ndof_in_parent_i;
             parent_dof_i++)
        {
          unsigned mapped_dof_i =
            Doftype_coarsen_map_coarse[block_i][parent_dof_i];
 
          if (is_master_block_preconditioner())
          {
            tmp_row_dist_pt[parent_dof_i] =
              Internal_block_distribution_pt[mapped_dof_i];
          }
          else
          {
            mapped_dof_i =
              Doftype_in_master_preconditioner_coarse[mapped_dof_i];
 
            tmp_row_dist_pt[parent_dof_i] =
              master_block_preconditioner_pt()->dof_block_distribution_pt(
                mapped_dof_i);
          }
        }
 
        Vector<LinearAlgebraDistribution*> tmp_col_dist_pt(ndof_in_parent_j, 0);
 
        for (unsigned parent_dof_j = 0; parent_dof_j < ndof_in_parent_j;
             parent_dof_j++)
        {
          unsigned mapped_dof_j =
            Doftype_coarsen_map_coarse[block_j][parent_dof_j];
 
          if (is_master_block_preconditioner())
          {
            tmp_col_dist_pt[parent_dof_j] =
              Internal_block_distribution_pt[mapped_dof_j];
          }
          else
          {
            mapped_dof_j =
              Doftype_in_master_preconditioner_coarse[mapped_dof_j];
            tmp_col_dist_pt[parent_dof_j] =
              master_block_preconditioner_pt()->dof_block_distribution_pt(
                mapped_dof_j);
          }
        }
 
        CRDoubleMatrixHelpers::concatenate_without_communication(
          tmp_row_dist_pt, tmp_col_dist_pt, tmp_blocks_pt, output_block);
 
        for (unsigned dof_i = 0; dof_i < ndof_in_parent_i; dof_i++)
        {
          for (unsigned dof_j = 0; dof_j < ndof_in_parent_j; dof_j++)
          {
            if (new_block[dof_i][dof_j])
            {
              delete tmp_blocks_pt(dof_i, dof_j);
            }
          }
        }
      }
    }
    else
    {
      CRDoubleMatrixHelpers::deep_copy(tmp_block_pt, output_block);
    }
  }
 
  //=============================================================================
  /// test function to check that every element in the block matrix
  /// (block_i,block_j) matches the corresponding element in the original matrix
  //=============================================================================
  template<typename MATRIX>
  void BlockPreconditioner<MATRIX>::block_matrix_test(
    const unsigned& block_i,
    const unsigned& block_j,
    const MATRIX* block_matrix_pt) const
  {
    // boolean flag to indicate whether test is passed
    bool check = true;
 
    // number of rows in matrix
    unsigned n_row = matrix_pt()->nrow();
 
    // number of columns in matrix
    unsigned n_col = matrix_pt()->ncol();
 
    // loop over rows of original matrix
    for (unsigned i = 0; i < n_row; i++)
    {
      // if this coefficient is associated with a block in this block
      // preconditioner
      if (static_cast<int>(block_i) == this->internal_block_number(i))
      {
        // loop over columns of original matrix
        for (unsigned j = 0; j < n_col; j++)
        {
          // if the coeeficient is associated with a block in this block
          // preconditioner
          if (static_cast<int>(block_j) == this->internal_block_number(j))
          {
            // check whether elements in original matrix and matrix of block
            // pointers match
            if (matrix_pt()->operator()(i, j) !=
                block_matrix_pt->operator()(internal_index_in_block(i),
                                            internal_index_in_block(j)))
            {
              check = false;
            }
          }
        }
      }
    }
 
    // throw error
    if (!check)
    {
      std::ostringstream error_message;
      error_message << "The require elements have not been successfully copied"
                    << " from the original matrix to the block matrices";
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  }
 
 
  template class BlockPreconditioner<CRDoubleMatrix>;
 
} // namespace oomph