pseudo_elastic_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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 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 "pseudo_elastic_preconditioner.h"

namespace oomph
{
 //=============================================================================
 /// \short Functions to create instances of optimal subsidiary operators for
 /// the PseudoElasticPreconditioner
 //=============================================================================
 namespace Pseudo_Elastic_Preconditioner_Subsidiary_Operator_Helper
 {
#ifdef OOMPH_HAS_HYPRE

  /// \short AMG w/ GS smoothing for the augmented elastic subsidiary linear
  /// systems
  Preconditioner* get_elastic_preconditioner_hypre()
  {
   HyprePreconditioner* hypre_preconditioner_pt = 
    new HyprePreconditioner("Hypre for diagonal blocks in pseudo-solid");
   hypre_preconditioner_pt->set_amg_iterations(2);
   hypre_preconditioner_pt->amg_using_simple_smoothing();
   if (MPI_Helpers::communicator_pt()->nproc()>1)
    {
     // Jacobi in parallel
     hypre_preconditioner_pt->amg_simple_smoother() = 0;
    }
   else
    {
     // Gauss Seidel in serial (was 3 actually...)
     hypre_preconditioner_pt->amg_simple_smoother() =1;
    }
   hypre_preconditioner_pt->hypre_method() = HyprePreconditioner::BoomerAMG;
   hypre_preconditioner_pt->amg_damping() = 1.0;
   hypre_preconditioner_pt->amg_coarsening() = 6;
   return hypre_preconditioner_pt;
  }


  /// \short AMG w/ GS smoothing for the augmented elastic subsidiary linear
  /// systems -- calls Hypre version to stay consistent with previous default
  Preconditioner* get_elastic_preconditioner()
  {
   return get_elastic_preconditioner_hypre();
  }
  
#endif
  
#ifdef OOMPH_HAS_TRILINOS

  /// \short TrilinosML smoothing for the augmented elastic 
  /// subsidiary linear systems
  Preconditioner* get_elastic_preconditioner_trilinos_ml()
  {
   TrilinosMLPreconditioner* trilinos_prec_pt=
    new TrilinosMLPreconditioner;
   return trilinos_prec_pt;
  }

  /// \short CG with diagonal preconditioner for the lagrange multiplier
  /// subsidiary linear systems.
  Preconditioner* get_lagrange_multiplier_preconditioner()
  {
   InnerIterationPreconditioner
    <TrilinosAztecOOSolver,MatrixBasedDiagPreconditioner>* prec_pt = 
    new InnerIterationPreconditioner
    <TrilinosAztecOOSolver,MatrixBasedDiagPreconditioner>;
   // Note: This makes CG a proper "inner iteration" for
   // which GMRES (may) no longer converge. We should really
   // use FGMRES or GMRESR for this. However, here the solver
   // is so good that it'll converge very quickly anyway
   // so there isn't much to be gained by limiting the number
   // of iterations...
   // prec_pt->max_iter() = 4;
   prec_pt->solver_pt()->solver_type() = TrilinosAztecOOSolver::CG;
   prec_pt->solver_pt()->disable_doc_time();
   return prec_pt;
  }
#endif
 }


 
 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////



 //=============================================================================
 // Setup method for the PseudoElasticPreconditioner.
 //=============================================================================
 void PseudoElasticPreconditioner::setup()
 {
  // clean
  this->clean_up_memory();
  
#ifdef PARANOID
  // paranoid check that meshes have been set
  if (Elastic_mesh_pt==0)
   {
    std::ostringstream error_message;
    error_message << "The elastic mesh must be set.";
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
    
   }
  if (Lagrange_multiplier_mesh_pt==0)
   {
    std::ostringstream error_message;
    error_message << "The Lagrange multiplier mesh must be set.";
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
  // set the meshes
  this->set_nmesh(2);
  this->set_mesh(0,Elastic_mesh_pt);
  this->set_mesh(1,Lagrange_multiplier_mesh_pt);
  
  // Figure out number of dof types
  unsigned n_solid_dof_types = 0;
  unsigned n_dof_types = 0;

  if (this->is_master_block_preconditioner())
   {
    // get the number of solid dof types from the first element
    n_solid_dof_types = this->ndof_types_in_mesh(0);
    
    // get the total number of dof types
    n_dof_types = n_solid_dof_types 
     + this->ndof_types_in_mesh(1);
   }
  else
   {
    n_dof_types = this->ndof_types();
    n_solid_dof_types = n_dof_types - (n_dof_types/3);
   }
#ifdef PARANOID
  if (n_dof_types%3 != 0)
   {
    std::ostringstream error_message;
    error_message << "This preconditioner requires DIM*3 types of DOF";
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
  // determine the dimension
  Dim = n_dof_types/3;
  
#ifdef PARANOID
  // Recast Jacobian matrix to CRDoubleMatrix
  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());

  if (cr_matrix_pt==0)
   {
    std::ostringstream error_message;
    error_message << "FSIPreconditioner only works with"
                  << " CRDoubleMatrix matrices" << std::endl;
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif

  // Setup the dof_list scheme for block_setup. 
  // The natural DOF types order is (in 3D):
  // 0 1 2 3  4  5  6  7  8  <- natural DOF type index
  // x y z xc yc zc lx ly lz <- DOf type
  // 
  // We need to group the directional displacements together:
  // 0 1  2 3  4  5  6  7  8  <- block index
  // x xc y yc xz zc lx ly lz <- DOF type
  //
  // The mapping required is:
  // 0 2 4 1 3 5 6 7 8
  //
  // In general we want:
  //
  //   dof_to_block_map[DOF type] = block type
  Vector<unsigned> dof_list_for_block_setup(n_dof_types);
  for (unsigned dim_i = 0; dim_i < Dim; dim_i++) 
   {
    // bulk solid dof types
    dof_list_for_block_setup[dim_i] = 2*dim_i;

    // constrained solid dof types
    dof_list_for_block_setup[dim_i + Dim] = 2*dim_i + 1;

    // lagr dof types
    dof_list_for_block_setup[dim_i + 2*Dim] = 2*Dim + dim_i;
   }

  // Setup the block ordering. We have the following:
  // Block types [0, 2*Dim) are the solid blocks.
  // Block types [2*Dim, 3*Dim) are the Lagrange multiplier blocks.
  // 
  // For d = 0, 1, 2,
  // Bulk solid doftypes: 2*d
  // Constrained solid doftypes: 2*d+1
  // Lgr doftyoes: 2*Dim + d
  this->block_setup(dof_list_for_block_setup);

  // Create dof_list for subsidiary preconditioner. This will be used later
  // in turn_into_subsidiary_block_preconditioner(...)
  Vector<unsigned> dof_list_for_subsidiary_prec(n_solid_dof_types);

  // The dof types are currently in the order (in 3D):
  // 0 1 2 3  4  5  6  7  8  <- DOF index
  // x y z xc yc zc lx ly lz <- DOF type
  // 
  // We need to group the directional displacements together:
  // x xc y yc xz zc
  //
  // The mapping required is:
  // 0 3 1 4 2 5
  //
  // In general we want:
  //    dof_list[subsidiary DOF type] = master DOF type
  for (unsigned dim_i = 0; dim_i < Dim; dim_i++) 
   {
    // bulk solid dof 
    dof_list_for_subsidiary_prec[2*dim_i] = dim_i;

    // constrained solid dof
    dof_list_for_subsidiary_prec[2*dim_i+1] = dim_i + Dim;
   }

  // Get the solid blocks
  DenseMatrix<CRDoubleMatrix*> solid_matrix_pt(n_solid_dof_types,
                                               n_solid_dof_types,0);

  for (unsigned row_i = 0; row_i < n_solid_dof_types; row_i++) 
   {
    for (unsigned col_i = 0; col_i < n_solid_dof_types; col_i++) 
     {
      solid_matrix_pt(row_i,col_i) = new CRDoubleMatrix;
      this->get_block(row_i,col_i,*solid_matrix_pt(row_i,col_i));
     }
   }

  // compute the scaling (if required)
  if (Use_inf_norm_of_s_scaling)
   {
    Scaling = CRDoubleMatrixHelpers::inf_norm(solid_matrix_pt);
   }
  else
   {
    Scaling = 1.0;
   }

  // Add the scaled identify matrix to the constrained solid blocks.
  for(unsigned d = 0; d < Dim; d++)
   {
    // Where is the constrained block located?
    unsigned block_i = 2*d+1;

    // Data from the constrained block.
    double* s_values = solid_matrix_pt(block_i,block_i)->value();
    int* s_column_index = solid_matrix_pt(block_i,block_i)->column_index();
    int* s_row_start = solid_matrix_pt(block_i,block_i)->row_start();
    int s_nrow_local = solid_matrix_pt(block_i,block_i)->nrow_local();
    int s_first_row = solid_matrix_pt(block_i,block_i)->first_row();

    // Loop through the diagonal entries and add the scaling.
    for (int i = 0; i < s_nrow_local; i++)
     {
      bool found = false;
      for (int j = s_row_start[i]; 
           j < s_row_start[i+1] && !found; j++)
       {
        if (s_column_index[j] == i + s_first_row)
         {
          s_values[j] += Scaling;
          found = true;
         } // if
       } // for row start
      
      // Check if the diagonal entry is found.
      if(!found)
       {
        std::ostringstream error_message;
        error_message << "The diagonal entry for the constained block("
                      << block_i<<","<<block_i<<")\n"
                      << "on local row " << i << " does not exist."
                      << std::endl;
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
       }
     } // for nrow_local
   } // for Dim 


  // setup the solid subsidiary preconditioner ////////////////////////////////
  // this preconditioner uses the full S matrix
  if (E_preconditioner_type == Exact_block_preconditioner)
   {
    ExactBlockPreconditioner<CRDoubleMatrix>* s_prec_pt = 
     new ExactBlockPreconditioner<CRDoubleMatrix>;

    // Add mesh (not actually used since this only acts as a subsidiary
    // preconditioner...
    //s_prec_pt->add_mesh(Elastic_mesh_pt);

    Vector<Vector<unsigned> > doftype_to_doftype_map(n_solid_dof_types,
                                                     Vector<unsigned>(1,0));

    for (unsigned i = 0; i < n_solid_dof_types; i++) 
     {
      doftype_to_doftype_map[i][0] = i;
     }

    s_prec_pt->turn_into_subsidiary_block_preconditioner(
     this,dof_list_for_subsidiary_prec,doftype_to_doftype_map);

    if (Elastic_subsidiary_preconditioner_function_pt != 0)
     {
      s_prec_pt->
       set_subsidiary_preconditioner_function
       (Elastic_subsidiary_preconditioner_function_pt);
     }

    // Set the replacement blocks.
    for (unsigned d = 0; d < Dim; d++) 
    {
      // The dof-block located in the block-ordering.
      unsigned block_i = 2*d+1;

      // The dof-block located in the dof-ordering.
      unsigned dof_block_i = Dim + d;
      this->set_replacement_dof_block(dof_block_i,dof_block_i,
                                      solid_matrix_pt(block_i,block_i));
    }

    s_prec_pt->Preconditioner::setup(matrix_pt());
    Elastic_preconditioner_pt = s_prec_pt;
   }
  // otherwise it is a block based preconditioner
  else 
   {
    GeneralPurposeBlockPreconditioner<CRDoubleMatrix>* s_prec_pt = 0;

    // set the block preconditioning method
    switch (E_preconditioner_type)
     {
     case Block_diagonal_preconditioner:
     {
      s_prec_pt = new BlockDiagonalPreconditioner<CRDoubleMatrix>;
     }
     break;
     case Block_upper_triangular_preconditioner:
     {
      BlockTriangularPreconditioner<CRDoubleMatrix>* block_triangular_prec_pt 
       = new BlockTriangularPreconditioner<CRDoubleMatrix>;
      block_triangular_prec_pt->upper_triangular();

      s_prec_pt = block_triangular_prec_pt;
     }
     break;
     case Block_lower_triangular_preconditioner:
     {
      BlockTriangularPreconditioner<CRDoubleMatrix>* block_triangular_prec_pt 
       = new BlockTriangularPreconditioner<CRDoubleMatrix>;
      block_triangular_prec_pt->lower_triangular();

      s_prec_pt = block_triangular_prec_pt;
     }
     break;
     default:
     {
      std::ostringstream error_msg;
      error_msg << "There is no such block based preconditioner.\n"
      << "Candidates are:\n"
      << "PseudoElasticPreconditioner::Block_diagonal_preconditioner\n"
      << "PseudoElasticPreconditioner::Block_upper_triangular_preconditioner\n"
      << "PseudoElasticPreconditioner::Block_lower_triangular_preconditioner\n";
      throw OomphLibError(error_msg.str(),
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
     }
     break;
     }


    // Add mesh (not actually used since this only acts as a subsidiary
    // preconditioner...
    //s_prec_pt->add_mesh(Elastic_mesh_pt);
    
    // The block to block map
    Vector<Vector<unsigned> > doftype_to_doftype_map(
      Dim,Vector<unsigned>(2,0));
    Vector<unsigned> s_prec_dof_to_block_map(Dim,0);

    unsigned tmp_index = 0;
    for (unsigned d = 0; d < Dim; d++) 
     {
      s_prec_dof_to_block_map[d] = d;
      doftype_to_doftype_map[d][0] = tmp_index++;
      doftype_to_doftype_map[d][1] = tmp_index++;
     }
    
    s_prec_pt->turn_into_subsidiary_block_preconditioner
     (this,dof_list_for_subsidiary_prec,doftype_to_doftype_map);

    if(Elastic_subsidiary_preconditioner_function_pt != 0)
     {
      s_prec_pt->set_subsidiary_preconditioner_function
       (Elastic_subsidiary_preconditioner_function_pt);
     }

    // Set the replacement blocks.
    for (unsigned d = 0; d < Dim; d++) 
    {
      // The dof-block located in the block-ordering.
      unsigned block_i = 2*d+1;

      // Then dof-block located in the dof-ordering.
      unsigned dof_block_i = Dim + d;
      this->set_replacement_dof_block(dof_block_i,dof_block_i,
                                      solid_matrix_pt(block_i,block_i));
    }

    s_prec_pt->set_dof_to_block_map(s_prec_dof_to_block_map);
    s_prec_pt->Preconditioner::setup(matrix_pt());

    Elastic_preconditioner_pt = s_prec_pt;
   }

  // No longer require the solid blocks
  for (unsigned row_i = 0; row_i < n_solid_dof_types; row_i++) 
   {
    for (unsigned col_i = 0; col_i < n_solid_dof_types; col_i++) 
     {
      delete solid_matrix_pt(row_i,col_i); solid_matrix_pt(row_i,col_i) = 0;
     }
   }
 
  // next setup the lagrange multiplier preconditioners
  Lagrange_multiplier_preconditioner_pt.resize(Dim);
  for (unsigned d = 0; d < Dim; d++)
   {
    CRDoubleMatrix* b_pt = new CRDoubleMatrix;
    this->get_block(2*Dim+d,2*d+1,*b_pt);

    // if a non default preconditioner is specified create 
    // the preconditioners
    if (Lagrange_multiplier_subsidiary_preconditioner_function_pt != 0)
     {
      Lagrange_multiplier_preconditioner_pt[d] = 
       (*Lagrange_multiplier_subsidiary_preconditioner_function_pt)();
     }
    // else use default superlu preconditioner
    else
     {
      Lagrange_multiplier_preconditioner_pt[d] = new SuperLUPreconditioner;
     }
    
    // and setup
    Lagrange_multiplier_preconditioner_pt[d]->setup(b_pt);
    delete b_pt; b_pt = 0;
   }

 }
 
 //=============================================================================
 /// \short Apply the elastic subsidiary preconditioner.
 //=============================================================================
 void PseudoElasticPreconditioner::elastic_preconditioner_solve
 (const DoubleVector& r, DoubleVector& z)
 {
  // apply the solid preconditioner
  Elastic_preconditioner_pt->preconditioner_solve(r,z);
 }

 //=============================================================================
 /// \short Apply the lagrange multiplier subsidiary preconditioner.
 //=============================================================================
 void PseudoElasticPreconditioner::
 lagrange_multiplier_preconditioner_solve(const DoubleVector& r,
                                          DoubleVector& z)
 {
  // apply the lagrange multiplier preconditioner
  for (unsigned d = 0; d < Dim; d++)
   {
    DoubleVector x;
    this->get_block_vector(Dim*2+d,r,x);
    DoubleVector y;
    Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(x,y);
    Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(y,x);
    unsigned nrow_local = x.nrow_local();
    double* x_pt = x.values_pt();
    for (unsigned i = 0; i < nrow_local; i++)
     {
      x_pt[i] = x_pt[i] * Scaling;
     }     
    this->return_block_vector(Dim*2+d,x,z);
   }
 }

 //=============================================================================
 /// \short Clears the memory.
 //=============================================================================
 void PseudoElasticPreconditioner::clean_up_memory()
 {
  // clean the block preconditioner base class memory
  this->clear_block_preconditioner_base();
  
  // delete the solid preconditioner
  delete Elastic_preconditioner_pt; Elastic_preconditioner_pt = 0;
  
  // delete the lagrange multiplier preconditioner pt
  unsigned sz = Lagrange_multiplier_preconditioner_pt.size();
  for (unsigned i = 0; i < sz; i++)
   {
    delete Lagrange_multiplier_preconditioner_pt[i];
    Lagrange_multiplier_preconditioner_pt[i] = 0;
   }
 }

 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////

 //=============================================================================
 // Setup method for the PseudoElasticPreconditioner.
 //=============================================================================
 void PseudoElasticPreconditionerOld::setup()
 {
  // clean
  this->clean_up_memory();
  
#ifdef PARANOID
  // paranoid check that meshes have been set
  if (Elastic_mesh_pt==0)
   {
    std::ostringstream error_message;
    error_message << "The elastic mesh must be set.";
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
    
   }
  if (Lagrange_multiplier_mesh_pt==0)
   {
    std::ostringstream error_message;
    error_message << "The Lagrange multiplier mesh must be set.";
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
  // set the mesh
  unsigned n_solid_dof_types = 0;
  unsigned n_dof_types = 0;
  this->set_mesh(0,Elastic_mesh_pt);
  this->set_mesh(1,Lagrange_multiplier_mesh_pt);
  if (this->is_master_block_preconditioner())
   {
    
    // get the number of solid dof types from the first element
    n_solid_dof_types = this->ndof_types_in_mesh(0);
    
    // get the total number of dof types
    n_dof_types = n_solid_dof_types 
     + this->ndof_types_in_mesh(1);
   }
  else
   {
    n_dof_types = this->ndof_types();
    n_solid_dof_types = n_dof_types - (n_dof_types/3);
   }
#ifdef PARANOID
  if (n_dof_types%3 != 0)
   {
    std::ostringstream error_message;
    error_message << "This preconditioner requires DIM*3 types of DOF";
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
  // determine the dimension
  Dim = n_dof_types/3;
  
  // Recast Jacobian matrix to CRDoubleMatrix
  CRDoubleMatrix* cr_matrix_pt = dynamic_cast<CRDoubleMatrix*>(matrix_pt());

#ifdef PARANOID
  if (cr_matrix_pt==0)
   {
    std::ostringstream error_message;
    error_message << "FSIPreconditioner only works with"
                  << " CRDoubleMatrix matrices" << std::endl;
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif

  // Setup the block ordering.
  this->block_setup();

  // Create dof_list for subsidiary preconditioner. This will be used later
  // in turn_into_subsidiary_block_preconditioner(...)
  Vector<unsigned> dof_list_for_subsidiary_prec(n_solid_dof_types);
  for (unsigned i = 0; i < n_solid_dof_types; i++)
   {
    dof_list_for_subsidiary_prec[i] = i;
   }

  // compute the scaling (if required)
  if (Use_inf_norm_of_s_scaling)
   {
    Vector<unsigned> dof_list(n_solid_dof_types);
    for (unsigned i = 0; i < n_solid_dof_types; i++)
     {
      dof_list[i] = i;
     }

    PseudoElasticPreconditionerScalingHelperOld* helper_pt = 
     new PseudoElasticPreconditionerScalingHelperOld(this,
                                                     cr_matrix_pt,
                                                     dof_list,
                                                     Elastic_mesh_pt,
                                                     comm_pt());
    Scaling = helper_pt->s_inf_norm();
    delete helper_pt; helper_pt = 0;
   }
  else
   {
    Scaling = 1.0;
   }
  
 
  // setup the solid subsidiary preconditioner ////////////////////////////////
  // this preconditioner uses the full S matrix
  if (E_preconditioner_type == Exact_block_preconditioner)
   {
    PseudoElasticPreconditionerSubsidiaryPreconditionerOld* s_prec_pt = 
     new PseudoElasticPreconditionerSubsidiaryPreconditionerOld;
    Vector<unsigned> dof_list(n_solid_dof_types);
    for (unsigned i = 0; i < n_solid_dof_types; i++)
     {
      dof_list[i] = i;
     }
    s_prec_pt->turn_into_subsidiary_block_preconditioner(this,dof_list);
    if (Elastic_subsidiary_preconditioner_function_pt != 0)
     {
      s_prec_pt->
       set_subsidiary_preconditioner_function
       (Elastic_subsidiary_preconditioner_function_pt);
     }
    s_prec_pt->scaling() = Scaling;
    s_prec_pt->Preconditioner::setup(matrix_pt());
    Elastic_preconditioner_pt = s_prec_pt;
   }
  
  // otherwise it is a block based preconditioner
  else 
   {
    // create the preconditioner
    PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld* 
     s_prec_pt = 
     new PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld;
    Vector<unsigned> dof_list(n_solid_dof_types);
    for (unsigned i = 0; i < n_solid_dof_types; i++)
     {
      dof_list[i] = i;
     }
    s_prec_pt->turn_into_subsidiary_block_preconditioner(this,dof_list);
    
    // set the subsidiary solve method
    if (Elastic_subsidiary_preconditioner_function_pt != 0)
     {
      s_prec_pt->
       set_subsidiary_preconditioner_function
       (Elastic_subsidiary_preconditioner_function_pt);
     }
    
    // set the scaling
    s_prec_pt->scaling() = Scaling;
    
    // set the block preconditioning method
    switch (E_preconditioner_type)
     {
     case Block_diagonal_preconditioner:
      s_prec_pt->use_block_diagonal_approximation();
      break;
     case Block_upper_triangular_preconditioner:
      s_prec_pt->use_upper_triangular_approximation();
      break;
     case Block_lower_triangular_preconditioner:
      s_prec_pt->use_lower_triangular_approximation();
      break;
     default:
      break;
     }
    
    // setup
    s_prec_pt->Preconditioner::setup(matrix_pt());
    Elastic_preconditioner_pt = s_prec_pt;

   }

  // next setup the lagrange multiplier preconditioners
  Lagrange_multiplier_preconditioner_pt.resize(Dim);
  for (unsigned d = 0; d < Dim; d++)
   {
    CRDoubleMatrix* b_pt = new CRDoubleMatrix;
    this->get_block(2*Dim+d,Dim+d,*b_pt);

    // if a non default preconditioner is specified create 
    // the preconditioners
    if (Lagrange_multiplier_subsidiary_preconditioner_function_pt != 0)
     {
      Lagrange_multiplier_preconditioner_pt[d] = 
       (*Lagrange_multiplier_subsidiary_preconditioner_function_pt)();
     }
    
    // else use default superlu preconditioner
    else
     {
      Lagrange_multiplier_preconditioner_pt[d] = new SuperLUPreconditioner;
     }
    
    // and setup
    Lagrange_multiplier_preconditioner_pt[d]->setup(b_pt);
    delete b_pt; b_pt = 0;
   }

 }
 
 //=============================================================================
 /// \short Apply the elastic subsidiary preconditioner.
 //=============================================================================
 void PseudoElasticPreconditionerOld::elastic_preconditioner_solve
 (const DoubleVector& r, DoubleVector& z)
 {
  // apply the solid preconditioner
  Elastic_preconditioner_pt->preconditioner_solve(r,z);
 }
 

 //=============================================================================
 /// \short Apply the lagrange multiplier subsidiary preconditioner.
 //=============================================================================
 void PseudoElasticPreconditionerOld::
 lagrange_multiplier_preconditioner_solve(const DoubleVector& r,
                                          DoubleVector& z)
 {
  // apply the lagrange multiplier preconditioner
  for (unsigned d = 0; d < Dim; d++)
   {
    DoubleVector x;
    this->get_block_vector(Dim*2+d,r,x);
    DoubleVector y;
    Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(x,y);
    Lagrange_multiplier_preconditioner_pt[d]->preconditioner_solve(y,x);
    unsigned nrow_local = x.nrow_local();
    double* x_pt = x.values_pt();
    for (unsigned i = 0; i < nrow_local; i++)
     {
      x_pt[i] = x_pt[i] * Scaling;
     }     
    this->return_block_vector(Dim*2+d,x,z);
   }
 }

 //=============================================================================
 /// \short Clears the memory.
 //=============================================================================
 void PseudoElasticPreconditionerOld::clean_up_memory()
 {
  // clean the block preconditioner base class memory
  this->clear_block_preconditioner_base();
  
  // delete the solid preconditioner
  delete Elastic_preconditioner_pt;
  Elastic_preconditioner_pt = 0;
  
  // delete the lagrange multiplier preconditioner pt
  unsigned sz = Lagrange_multiplier_preconditioner_pt.size();
  for (unsigned i = 0; i < sz; i++)
   {
    delete Lagrange_multiplier_preconditioner_pt[i];
    Lagrange_multiplier_preconditioner_pt[i] = 0;
   }
 }



 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////



 //=============================================================================
 /// \short Setup the preconditioner
 //=============================================================================
 void PseudoElasticPreconditionerSubsidiaryPreconditionerOld::
 setup()
 {
  // clean memory
  this->clean_up_memory();
     
#ifdef PARANOID
  // paranoid check that this preconditioner has an even number of DOF types
  if (this->ndof_types()%2 != 0)
   {
    std::ostringstream error_message;
    error_message
     << "This SUBSIDIARY preconditioner requires an even number of "
     << "types of DOF";
    throw OomphLibError(
     error_message.str(),
     OOMPH_CURRENT_FUNCTION,
     OOMPH_EXCEPTION_LOCATION);
   }
#endif
     
  // assemble dof_to_block_map
  unsigned ndof_types = this->ndof_types();
  Vector<unsigned> dof_to_block_map(ndof_types,0);
  for (unsigned i = ndof_types/2; i < ndof_types; i++)
   {
    dof_to_block_map[i] = 1;
   }

  this->block_setup(dof_to_block_map);

  // get block 11
  CRDoubleMatrix* s11_pt = new CRDoubleMatrix;
  this->get_block(1,1,*s11_pt);

  // add the scaled identity matrix to block 11
  double* s11_values = s11_pt->value();
  int* s11_column_index = s11_pt->column_index();
  int* s11_row_start = s11_pt->row_start();
  int s11_nrow_local = s11_pt->nrow_local();
  int s11_first_row = s11_pt->first_row();
  for (int i = 0; i < s11_nrow_local; i++)
   {
    bool found = false;
    for (int j = s11_row_start[i]; 
         j < s11_row_start[i+1] && !found; j++)
     {
      if (s11_column_index[j] == i + s11_first_row)
       {
        s11_values[j] += Scaling;
        found = true;
       }
     }
   }

  VectorMatrix<BlockSelector> required_blocks(2,2);
  const bool want_block = true;
  for (unsigned b_i = 0; b_i < 2; b_i++) 
  {
    for (unsigned b_j = 0; b_j < 2; b_j++) 
    {
      required_blocks[b_i][b_j].select_block(b_i,b_j,want_block);
    }
  }

  required_blocks[1][1].set_replacement_block_pt(s11_pt);

  CRDoubleMatrix s_prec_pt = this->get_concatenated_block(required_blocks);

  delete s11_pt; s11_pt = 0;
   
  // setup the preconditioner
  if (Subsidiary_preconditioner_function_pt != 0)
   {
    Preconditioner_pt = (*Subsidiary_preconditioner_function_pt)();
   }
  else
   {
    Preconditioner_pt = new SuperLUPreconditioner;
   }
  Preconditioner_pt->setup(&s_prec_pt);
 }
   
 //=============================================================================
 /// \short Apply the preconditioner.
 //=============================================================================
 void PseudoElasticPreconditionerSubsidiaryPreconditionerOld::
 preconditioner_solve(const DoubleVector& r, DoubleVector& z)
 {
  DoubleVector x;
  this->get_block_ordered_preconditioner_vector(r,x);
  DoubleVector y;
  Preconditioner_pt->preconditioner_solve(x,y);
  this->return_block_ordered_preconditioner_vector(y,z);
 }



 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////


 
 //=============================================================================
 /// clean up the memory
 //=============================================================================
 void PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld::
 clean_up_memory()
 {
  //number of block types
  unsigned n_block = Diagonal_block_preconditioner_pt.size();
   
  //delete diagonal blocks
  for (unsigned i = 0 ; i < n_block; i++)
   {
    delete Diagonal_block_preconditioner_pt[i];
    Diagonal_block_preconditioner_pt[i] = 0;
    if (Method == 1)
     {
      for (unsigned j = i+1; j < n_block; j++)
       {
        delete Off_diagonal_matrix_vector_products(i,j);
        Off_diagonal_matrix_vector_products(i,j) = 0;
       }
     }
    else if (Method == 2)
     {
      for (unsigned j = 0; j < i; j++)
       {
        delete Off_diagonal_matrix_vector_products(i,j);
        Off_diagonal_matrix_vector_products(i,j) = 0;
       }
     }
   }
   
  // clean up the block preconditioner
  this->clear_block_preconditioner_base();
 }
 
 //=============================================================================
 /// \short Setup the preconditioner.
 //=============================================================================
 void PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld::
 setup()
 {
  // clean the memory
  this->clean_up_memory();
  
  // determine the number of DOF types
  unsigned n_dof_types =  this->ndof_types();
  
#ifdef PARANOID
  // must be Dim*2 dof types
  if (n_dof_types%2 != 0)
   {
    std::ostringstream error_message;
    error_message << "This preconditioner requires DIM*3 types of DOF";
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
#endif
  
  // store the dimension of the problem
  unsigned dim = n_dof_types/2;
  
  // assemble the dof to block lookup scheme
  Vector<unsigned> dof_to_block_map(n_dof_types,0);
  for (unsigned d = 0; d < dim; d++)
   {
    dof_to_block_map[d] = d;
    dof_to_block_map[d+dim] = d;
   }
    
  //setup the blocks look up schemes
  this->block_setup(dof_to_block_map);

  // Storage for the diagonal block preconditioners
  Diagonal_block_preconditioner_pt.resize(dim);
  
  // storage for the off diagonal matrix vector products
  Off_diagonal_matrix_vector_products.resize(dim,dim,0);

  // setup the subsidiary preconditioners
  for (unsigned d = 0; d < dim; d++)
   {
    Vector<unsigned> dof_list(2);
    dof_list[0]=d;
    dof_list[1]=d+dim;

    Diagonal_block_preconditioner_pt[d] = new 
     PseudoElasticPreconditionerSubsidiaryPreconditionerOld; 
    Diagonal_block_preconditioner_pt[d]->
     turn_into_subsidiary_block_preconditioner(this,dof_list);
    if (Subsidiary_preconditioner_function_pt != 0)
     {
      Diagonal_block_preconditioner_pt[d]->
       set_subsidiary_preconditioner_function
       (Subsidiary_preconditioner_function_pt);
     }
    Diagonal_block_preconditioner_pt[d]->scaling() = Scaling;

    Diagonal_block_preconditioner_pt[d]->
     Preconditioner::setup(matrix_pt());
    
    // the preconditioning method.\n
    // 0 - block diagonal\n
    // 1 - upper triangular\n
    // 2 - lower triangular\n 
    // next setup the off diagonal mat vec operators if required
    if (Method == 1 || Method == 2)
     {
      unsigned l = d+1;
      unsigned u = dim;
      if (Method==2)
       {
        l = 0;
        u = d;
       }
      for (unsigned j = l; j < u; j++)
       {
        CRDoubleMatrix* block_matrix_pt = new CRDoubleMatrix;
        this->get_block(d,j,*block_matrix_pt);
        Off_diagonal_matrix_vector_products(d,j) 
         = new MatrixVectorProduct();
//        Off_diagonal_matrix_vector_products(d,j)->setup(block_matrix_pt);
        this->setup_matrix_vector_product(Off_diagonal_matrix_vector_products(d,j),
                                          block_matrix_pt,j);

        delete block_matrix_pt; block_matrix_pt = 0;
       }
     }
   } // setup the subsidiary preconditioner.
 } // preconditioner setup
 
 //=============================================================================
 /// Apply preconditioner to r
 //=============================================================================
 void PseudoElasticPreconditionerSubsidiaryBlockPreconditionerOld::
 preconditioner_solve(const DoubleVector &res, DoubleVector &z)
 {
  // copy r
  DoubleVector r(res);

  unsigned n_block;
   
  // Cache umber of block types (also the spatial DIM)
  n_block = this->nblock_types();

  // loop parameters
  int start = n_block-1;
  int end = -1;
  int step = -1;
  if (Method != 1)
   {
    start = 0;
    end = n_block;
    step = 1;
   }

  // the preconditioning method.
  // 0 - block diagonal
  // 1 - upper triangular
  // 2 - lower triangular
  //
  // loop over the DIM
  // 
  // For Method = 0 or 2 (diagonal, lower)
  // start = 2, end = -1, step = -1
  // i = 2,1,0
  // 
  // For Method = 1 (upper)
  // start = 0, end = 3 step = 1
  // i = 0, 1, 2
  for (int i = start; i != end; i+=step)
   {
       
    // solve
    Diagonal_block_preconditioner_pt[i]->preconditioner_solve(r,z);
       
    // if upper or lower triangular
    if (Method != 0)
     {

      // substitute
      //
      for (int j = i + step; j !=end; j+=step)
       {
        DoubleVector x;
        this->get_block_vector(i,z,x);
        DoubleVector y;
        Off_diagonal_matrix_vector_products(j,i)->multiply(x,y);
        x.clear();
        this->get_block_vector(j,r,x);
        x -= y;
        this->return_block_vector(j,x,r);
       } // substitute
     } // if upper or lower
   } // for loop over DIM
 } // Block preconditioner solve
} // namespace oomph