matrices.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// This header file contains classes and inline function definitions for
// matrices and their derived types
 
// Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_MATRICES_HEADER
#define OOMPH_MATRICES_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif
 
 
// Needed for g++ in some cases
#include <iomanip>
 
// oomph-lib headers
#include "Vector.h"
#include "oomph_utilities.h"
#include "linear_algebra_distribution.h"
#include "double_vector.h"
 
 
#ifdef OOMPH_HAS_TRILINOS
#include "trilinos_helpers.h"
#endif
 
namespace oomph
{
// Initialise dense pointer-based matrices/tensors?
#define OOMPH_INITIALISE_DENSE_MATRICES
#undef OOMPH_INITIALISE_DENSE_MATRICES
 
  //=================================================================
  //=================================================================
  template<class T, class MATRIX_TYPE>
  class Matrix
  {
  protected:
    void range_check(const unsigned long& i, const unsigned long& j) const
    {
      if (i >= nrow())
      {
        std::ostringstream error_message;
        error_message << "Range Error: i=" << i << " is not in the range (0,"
                      << nrow() - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else if (j >= ncol())
      {
        std::ostringstream error_message;
        error_message << "Range Error: j=" << j << " is not in the range (0,"
                      << ncol() - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
 
 
  public:
    Matrix() {}
 
    Matrix(const Matrix& matrix) = delete;
 
    void operator=(const Matrix&) = delete;
 
    virtual ~Matrix() {}
 
    virtual unsigned long nrow() const = 0;
 
    virtual unsigned long ncol() const = 0;
 
    inline T operator()(const unsigned long& i, const unsigned long& j) const
    {
      return static_cast<MATRIX_TYPE const*>(this)->get_entry(i, j);
    }
 
    inline T& operator()(const unsigned long& i, const unsigned long& j)
    {
      return static_cast<MATRIX_TYPE*>(this)->entry(i, j);
    }
 
    virtual void output(std::ostream& outfile) const
    {
      throw OomphLibError(
        "Output function is not implemented for this matrix class",
        OOMPH_CURRENT_FUNCTION,
        OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void output_bottom_right_zero_helper(
      std::ostream& outfile) const = 0;
 
    virtual void sparse_indexed_output_helper(std::ostream& outfile) const = 0;
 
 
    void sparse_indexed_output(
      std::ostream& outfile,
      const unsigned& precision = 0,
      const bool& output_bottom_right_zero = false) const
    {
      // Implemented as a wrapper around "sparse_indexed_output(std::ostream)"
      // so that only one output helper function is needed in derived classes.
 
      // We can't have separate functions for only "output_bottom_right_zero"
      // because people often write false as "0" and then C++ would pick the
      // wrong function.
 
      // If requested set the new precision and store the previous value.
      unsigned old_precision = 0;
      if (precision != 0)
      {
        old_precision = outfile.precision();
        outfile.precision(precision);
      }
 
      // Output as normal using the helper function defined in each matrix class
      sparse_indexed_output_helper(outfile);
 
      // If requested and there is no output for the last entry then output a
      // zero entry.
      if (output_bottom_right_zero && ncol() > 0 && nrow() > 0)
      {
        // Output as normal using the helper function defined
        // in each matrix class
        output_bottom_right_zero_helper(outfile);
      }
 
      // Restore the old value of the precision if we changed it
      if (precision != 0)
      {
        outfile.precision(old_precision);
      }
    }
 
    void sparse_indexed_output(
      std::string filename,
      const unsigned& precision = 0,
      const bool& output_bottom_right_zero = false) const
    {
      // Implemented as a wrapper around "sparse_indexed_output(std::ostream)"
      // so that only one output function needs to be written in matrix
      // subclasses.
 
      // Open file
      std::ofstream some_file(filename.c_str());
 
      // Output as normal
      sparse_indexed_output(some_file, precision, output_bottom_right_zero);
 
      // Close file
      some_file.close();
    }
  };
 
 
 
 
  // Forward definition of the linear solver class
  class LinearSolver;
 
  //=============================================================================
  //=============================================================================
  class DoubleMatrixBase
  {
  protected:
    // Pointer to a linear solver
    LinearSolver* Linear_solver_pt;
 
    // Pointer to a default linear solver
    LinearSolver* Default_linear_solver_pt;
 
  public:
    DoubleMatrixBase() : Linear_solver_pt(0), Default_linear_solver_pt(0) {}
 
    DoubleMatrixBase(const DoubleMatrixBase& matrix) = delete;
 
    void operator=(const DoubleMatrixBase&) = delete;
 
    virtual unsigned long nrow() const = 0;
 
    virtual unsigned long ncol() const = 0;
 
    virtual ~DoubleMatrixBase() {}
 
    virtual double operator()(const unsigned long& i,
                              const unsigned long& j) const = 0;
 
 
    LinearSolver*& linear_solver_pt()
    {
      return Linear_solver_pt;
    }
 
    LinearSolver* const& linear_solver_pt() const
    {
      return Linear_solver_pt;
    }
 
    void solve(DoubleVector& rhs);
 
    void solve(const DoubleVector& rhs, DoubleVector& soln);
 
    void solve(Vector<double>& rhs);
 
    void solve(const Vector<double>& rhs, Vector<double>& soln);
 
    virtual void residual(const DoubleVector& x,
                          const DoubleVector& b,
                          DoubleVector& residual_)
    {
      // compute residual = Ax
      this->multiply(x, residual_);
 
      // set residual to -residual (-Ax)
      unsigned nrow_local = residual_.nrow_local();
      double* residual_pt = residual_.values_pt();
      for (unsigned i = 0; i < nrow_local; i++)
      {
        residual_pt[i] = -residual_pt[i];
      }
 
      // set residual = b + residuals
      residual_ += b;
    }
 
    virtual double max_residual(const DoubleVector& x, const DoubleVector& rhs)
    {
      DoubleVector res;
      residual(x, rhs, res);
      return res.max();
    }
 
    virtual void multiply(const DoubleVector& x, DoubleVector& soln) const = 0;
 
    virtual void multiply_transpose(const DoubleVector& x,
                                    DoubleVector& soln) const = 0;
 
    // virtual void matrix_reduction(const double &alpha,
    //                              DoubleMatrixBase& reduced_matrix)=0;
  };
 
 
 
 
  //======================================================================
  //======================================================================
  template<class T>
  class DenseMatrix : public Matrix<T, DenseMatrix<T>>
  {
  protected:
    T* Matrixdata;
 
    unsigned long N;
 
    unsigned long M;
 
  public:
    DenseMatrix() : Matrixdata(0), N(0), M(0) {}
 
    DenseMatrix(const DenseMatrix& source_matrix)
    {
      // Set row and column lengths
      N = source_matrix.nrow();
      M = source_matrix.ncol();
      // Assign space for the data
      Matrixdata = new T[N * M];
      // Copy the data across from the other matrix
      for (unsigned long i = 0; i < N; i++)
      {
        for (unsigned long j = 0; j < M; j++)
        {
          Matrixdata[M * i + j] = source_matrix(i, j);
        }
      }
    }
 
    DenseMatrix& operator=(const DenseMatrix& source_matrix)
    {
      // Don't create a new matrix if the assignment is the identity
      if (this != &source_matrix)
      {
        // Check row and column length
        unsigned long n = source_matrix.nrow();
        unsigned long m = source_matrix.ncol();
        if ((N != n) || (M != m))
        {
          resize(n, m);
        }
        // Copy entries across from the other matrix
        for (unsigned long i = 0; i < N; i++)
        {
          for (unsigned long j = 0; j < M; j++)
          {
            (*this)(i, j) = source_matrix(i, j);
          }
        }
      }
      // Return reference to object itself (i.e. de-reference this pointer)
      return *this;
    }
 
    inline T& entry(const unsigned long& i, const unsigned long& j)
    {
#ifdef RANGE_CHECKING
      this->range_check(i, j);
#endif
      return Matrixdata[M * i + j];
    }
 
    inline T get_entry(const unsigned long& i, const unsigned long& j) const
    {
#ifdef RANGE_CHECKING
      this->range_check(i, j);
#endif
      return Matrixdata[M * i + j];
    }
 
    DenseMatrix(const unsigned long& n);
 
    DenseMatrix(const unsigned long& n, const unsigned long& m);
 
    DenseMatrix(const unsigned long& n,
                const unsigned long& m,
                const T& initial_val);
 
    virtual ~DenseMatrix()
    {
      delete[] Matrixdata;
      Matrixdata = 0;
    }
 
    inline unsigned long nrow() const
    {
      return N;
    }
 
    inline unsigned long ncol() const
    {
      return M;
    }
 
    void resize(const unsigned long& n)
    {
      resize(n, n);
    }
 
    void resize(const unsigned long& n, const unsigned long& m);
 
    void resize(const unsigned long& n,
                const unsigned long& m,
                const T& initial_value);
 
    void initialise(const T& val)
    {
      for (unsigned long i = 0; i < (N * M); ++i)
      {
        Matrixdata[i] = val;
      }
    }
 
    void output(std::ostream& outfile) const;
 
    void output(std::string filename) const;
 
    void indexed_output(std::ostream& outfile) const;
 
    void indexed_output(std::string filename) const;
 
    void output_bottom_right_zero_helper(std::ostream& outfile) const;
 
    void sparse_indexed_output_helper(std::ostream& outfile) const;
  };
 
 
 
 
  //================================================================
  //=================================================================
  template<class T, class MATRIX_TYPE>
  class SparseMatrix : public Matrix<T, MATRIX_TYPE>
  {
  protected:
    T* Value;
 
    unsigned long N;
 
    unsigned long M;
 
    unsigned long Nnz;
 
    static T Zero;
 
  public:
    SparseMatrix() : Value(0), N(0), M(0), Nnz(0) {}
 
    SparseMatrix(const SparseMatrix& source_matrix)
    {
      // Number of nonzero entries
      Nnz = source_matrix.nnz();
 
      // Number of rows
      N = source_matrix.nrow();
 
      // Number of columns
      M = source_matrix.ncol();
 
      // Values stored in C-style array
      Value = new T[Nnz];
 
      // Assign the values
      for (unsigned long i = 0; i < Nnz; i++)
      {
        Value[i] = source_matrix.value()[i];
      }
    }
 
    void operator=(const SparseMatrix&) = delete;
 
    virtual ~SparseMatrix()
    {
      delete[] Value;
      Value = 0;
    }
 
    T* value()
    {
      return Value;
    }
 
    const T* value() const
    {
      return Value;
    }
 
    inline unsigned long nrow() const
    {
      return N;
    }
 
    inline unsigned long ncol() const
    {
      return M;
    }
 
    inline unsigned long nnz() const
    {
      return Nnz;
    }
 
    virtual void output_bottom_right_zero_helper(std::ostream& outfile) const
    {
      std::string error_message = "SparseMatrix::output_bottom_right_zero_"
                                  "helper() is a virtual function.\n";
      error_message +=
        "It must be overloaded for specific sparse matrix storage formats\n";
 
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual void sparse_indexed_output_helper(std::ostream& outfile) const
    {
      std::string error_message =
        "SparseMatrix::sparse_indexed_output_helper() is a virtual function.\n";
      error_message +=
        "It must be overloaded for specific sparse matrix storage formats\n";
 
      throw OomphLibError(
        error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
  };
 
 
  //======================================================================
  //=====================================================================
  template<class T>
  class CRMatrix : public SparseMatrix<T, CRMatrix<T>>
  {
  public:
    CRMatrix() : SparseMatrix<T, CRMatrix<T>>()
    {
      Column_index = 0;
      Row_start = 0;
    }
 
 
    CRMatrix(const Vector<T>& value,
             const Vector<int>& column_index_,
             const Vector<int>& row_start_,
             const unsigned long& n,
             const unsigned long& m)
      : SparseMatrix<T, CRMatrix<T>>()
    {
      Column_index = 0;
      Row_start = 0;
      build(value, column_index_, row_start_, n, m);
    }
 
    CRMatrix(const CRMatrix& source_matrix)
      : SparseMatrix<T, CRMatrix<T>>(source_matrix)
    {
      // NNz, N and M are set the the copy constructor of the SparseMatrix
      // called above
      // Column indices stored in C-style array
      Column_index = new int[this->Nnz];
 
      // Assign:
      for (unsigned long i = 0; i < this->Nnz; i++)
      {
        Column_index[i] = source_matrix.column_index()[i];
      }
 
      // Row start:
      Row_start = new int[this->N + 1];
 
      // Assign:
      for (unsigned long i = 0; i <= this->N; i++)
      {
        Row_start[i] = source_matrix.row_start()[i];
      }
    }
 
    void operator=(const CRMatrix&) = delete;
 
    virtual ~CRMatrix()
    {
      delete[] Column_index;
      Column_index = 0;
      delete[] Row_start;
      Row_start = 0;
    }
 
    T get_entry(const unsigned long& i, const unsigned long& j) const
    {
#ifdef RANGE_CHECKING
      this->range_check(i, j);
#endif
      for (long k = Row_start[i]; k < Row_start[i + 1]; k++)
      {
        if (unsigned(Column_index[k]) == j)
        {
          return this->Value[k];
        }
      }
      return this->Zero;
    }
 
    T& entry(const unsigned long& i, const unsigned long& j)
    {
      std::string error_string =
        "Non-const access not provided for the CRMatrix<T> class\n";
      error_string +=
        "It is not possible to use round-bracket access: M(i,j)\n";
      error_string += "if M is not declared as const.\n";
      error_string += "The solution (albeit ugly) is to create const reference "
                      "to the matrix\n";
      error_string += " const CRMatrix<T>& read_M = M;\n";
      error_string += "Then read_M(i,j) is permitted\n";
 
      throw OomphLibError(
        error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
 
      // Dummy return
      T dummy;
      return dummy;
    }
 
    int* row_start()
    {
      return Row_start;
    }
 
    const int* row_start() const
    {
      return Row_start;
    }
 
    int* column_index()
    {
      return Column_index;
    }
 
    const int* column_index() const
    {
      return Column_index;
    }
 
    void output_bottom_right_zero_helper(std::ostream& outfile) const
    {
      int last_row_local = this->N - 1;
      int last_col = this->M - 1;
 
      // Use this strange thingy because of the CRTP discussed above.
      T last_value = this->operator()(last_row_local, last_col);
 
      if (last_value == T(0))
      {
        outfile << last_row_local << " " << last_col << " " << T(0)
                << std::endl;
      }
    }
 
    void sparse_indexed_output_helper(std::ostream& outfile) const
    {
      for (unsigned long i = 0; i < this->N; i++)
      {
        for (long j = Row_start[i]; j < Row_start[i + 1]; j++)
        {
          outfile << i << " " << Column_index[j] << " " << this->Value[j]
                  << std::endl;
        }
      }
    }
 
    void clean_up_memory();
 
    void build(const Vector<T>& value,
               const Vector<int>& column_index,
               const Vector<int>& row_start,
               const unsigned long& n,
               const unsigned long& m);
 
 
    void build_without_copy(T* value,
                            int* column_index,
                            int* row_start,
                            const unsigned long& nnz,
                            const unsigned long& n,
                            const unsigned long& m);
 
 
  protected:
    int* Column_index;
 
    int* Row_start;
  };
 
 
  // Forward definition for the superlu solver
  class SuperLUSolver;
 
 
  //=============================================================================
  //=============================================================================
  class CRDoubleMatrix : public Matrix<double, CRDoubleMatrix>,
                         public DoubleMatrixBase,
                         public DistributableLinearAlgebraObject
  {
  public:
    CRDoubleMatrix();
 
    CRDoubleMatrix(const LinearAlgebraDistribution* distribution_pt,
                   const unsigned& ncol,
                   const Vector<double>& value,
                   const Vector<int>& column_index,
                   const Vector<int>& row_start);
 
    CRDoubleMatrix(const LinearAlgebraDistribution* distribution_pt);
 
    CRDoubleMatrix(const CRDoubleMatrix& matrix);
 
    void operator=(const CRDoubleMatrix&) = delete;
 
    virtual ~CRDoubleMatrix();
 
    const Vector<int> get_index_of_diagonal_entries() const
    {
      // Check to see if the vector has been set up
      if (Index_of_diagonal_entries.size() == 0)
      {
        // Make the warning
        std::string err_strng =
          "The Index_of_diagonal_entries vector has not been ";
        err_strng += "set up yet. Run sort_entries() to set this vector up.";
 
        // Throw the warning
        OomphLibWarning(err_strng,
                        "CRDoubleMatrix::get_index_of_diagonal_entries()",
                        OOMPH_EXCEPTION_LOCATION);
      }
 
      // Return the vector
      return Index_of_diagonal_entries;
    } // End of index_of_diagonal_entries
 
    struct CRDoubleMatrixComparisonHelper
    {
      // Define the comparison operator
      bool operator()(const std::pair<int, double>& pair_1,
                      const std::pair<int, double>& pair_2)
      {
        // If the first argument of pair_1 is less than the first argument of
        // pair_2 then return TRUE otherwise return FALSE
        return (pair_1.first < pair_2.first);
      }
    } Comparison_struct;
 
    bool entries_are_sorted(const bool& doc_unordered_entries = false) const;
 
    void sort_entries();
 
    void build(const LinearAlgebraDistribution* distribution_pt,
               const unsigned& ncol,
               const Vector<double>& value,
               const Vector<int>& column_index,
               const Vector<int>& row_start);
 
    void build(const LinearAlgebraDistribution* distribution_pt);
 
    void build(const unsigned& ncol,
               const Vector<double>& value,
               const Vector<int>& column_index,
               const Vector<int>& row_start);
 
    void build_without_copy(const unsigned& ncol,
                            const unsigned& nnz,
                            double* value,
                            int* column_index,
                            int* row_start);
 
    void redistribute(const LinearAlgebraDistribution* const& dist_pt);
 
    void clear();
 
    inline unsigned long nrow() const
    {
      return DistributableLinearAlgebraObject::nrow();
    }
 
    inline unsigned long ncol() const
    {
      return CR_matrix.ncol();
    }
 
    void output_bottom_right_zero_helper(std::ostream& outfile) const
    {
      CR_matrix.output_bottom_right_zero_helper(outfile);
    }
 
    void sparse_indexed_output_helper(std::ostream& outfile) const
    {
      CR_matrix.sparse_indexed_output_helper(outfile);
    }
 
    void sparse_indexed_output_with_offset(std::string filename)
    {
      // Get offset
      unsigned first_row = distribution_pt()->first_row();
 
      // Open file
      std::ofstream some_file;
      some_file.open(filename.c_str());
      unsigned n = nrow_local();
      for (unsigned long i = 0; i < n; i++)
      {
        for (long j = row_start()[i]; j < row_start()[i + 1]; j++)
        {
          some_file << first_row + i << " " << column_index()[j] << " "
                    << value()[j] << std::endl;
        }
      }
      some_file.close();
    }
 
    inline double operator()(const unsigned long& i,
                             const unsigned long& j) const
    {
      return CR_matrix.get_entry(i, j);
    }
 
    int* row_start()
    {
      return CR_matrix.row_start();
    }
 
    const int* row_start() const
    {
      return CR_matrix.row_start();
    }
 
    int* column_index()
    {
      return CR_matrix.column_index();
    }
 
    const int* column_index() const
    {
      return CR_matrix.column_index();
    }
 
    double* value()
    {
      return CR_matrix.value();
    }
 
    const double* value() const
    {
      return CR_matrix.value();
    }
 
    inline unsigned long nnz() const
    {
      return CR_matrix.nnz();
    }
 
    virtual void ludecompose();
 
    virtual void lubksub(DoubleVector& rhs);
 
    void multiply(const DoubleVector& x, DoubleVector& soln) const;
 
    void multiply_transpose(const DoubleVector& x, DoubleVector& soln) const;
 
    void multiply(const CRDoubleMatrix& matrix_in,
                  CRDoubleMatrix& result) const;
 
    void matrix_reduction(const double& alpha, CRDoubleMatrix& reduced_matrix);
 
    unsigned& serial_matrix_matrix_multiply_method()
    {
      return Serial_matrix_matrix_multiply_method;
    }
 
    const unsigned& serial_matrix_matrix_multiply_method() const
    {
      return Serial_matrix_matrix_multiply_method;
    }
 
    unsigned& distributed_matrix_matrix_multiply_method()
    {
      return Distributed_matrix_matrix_multiply_method;
    }
 
    const unsigned& distributed_matrix_matrix_multiply_method() const
    {
      return Distributed_matrix_matrix_multiply_method;
    }
 
    bool built() const
    {
      return Built;
    }
 
    CRDoubleMatrix* global_matrix() const;
 
    void get_matrix_transpose(CRDoubleMatrix* result) const;
 
    double inf_norm() const;
 
    Vector<double> diagonal_entries() const;
 
    void add(const CRDoubleMatrix& matrix_in,
             CRDoubleMatrix& result_matrix) const;
 
  private:
    Vector<int> Index_of_diagonal_entries;
 
    unsigned Serial_matrix_matrix_multiply_method;
 
    unsigned Distributed_matrix_matrix_multiply_method;
 
    CRMatrix<double> CR_matrix;
 
    bool Built;
  };
 
 
 
 
  // Forward definition of the DenseLU class
  class DenseLU;
 
  //=================================================================
  //=================================================================
  class DenseDoubleMatrix : public DoubleMatrixBase, public DenseMatrix<double>
  {
  public:
    DenseDoubleMatrix();
 
    DenseDoubleMatrix(const unsigned long& n);
 
    DenseDoubleMatrix(const unsigned long& n, const unsigned long& m);
 
    DenseDoubleMatrix(const unsigned long& n,
                      const unsigned long& m,
                      const double& initial_val);
 
    DenseDoubleMatrix(const DenseDoubleMatrix& matrix) = delete;
 
    void operator=(const DenseDoubleMatrix&) = delete;
 
    inline unsigned long nrow() const
    {
      return DenseMatrix<double>::nrow();
    }
 
    inline unsigned long ncol() const
    {
      return DenseMatrix<double>::ncol();
    }
 
    inline double operator()(const unsigned long& i,
                             const unsigned long& j) const
    {
      return DenseMatrix<double>::get_entry(i, j);
    }
 
    inline double& operator()(const unsigned long& i, const unsigned long& j)
    {
      return DenseMatrix<double>::entry(i, j);
    }
 
    virtual ~DenseDoubleMatrix();
 
    virtual void ludecompose();
 
    virtual void lubksub(DoubleVector& rhs);
 
    virtual void lubksub(Vector<double>& rhs);
 
    void eigenvalues_by_jacobi(Vector<double>& eigen_val,
                               DenseMatrix<double>& eigen_vect) const;
 
    void multiply(const DoubleVector& x, DoubleVector& soln) const;
 
    void multiply_transpose(const DoubleVector& x, DoubleVector& soln) const;
 
    void matrix_reduction(const double& alpha,
                          DenseDoubleMatrix& reduced_matrix);
 
    void multiply(const DenseDoubleMatrix& matrix_in,
                  DenseDoubleMatrix& result);
  };
 
 
 
  //=================================================================
  //=================================================================
  template<class T>
  class RankThreeTensor
  {
  private:
    T* Tensordata;
 
    unsigned N;
 
    unsigned M;
 
    unsigned P;
 
    void range_check(const unsigned long& i,
                     const unsigned long& j,
                     const unsigned long& k) const
    {
      if (i >= N)
      {
        std::ostringstream error_message;
        error_message << "Range Error: i=" << i << " is not in the range (0,"
                      << N - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else if (j >= M)
      {
        std::ostringstream error_message;
        error_message << "Range Error: j=" << j << " is not in the range (0,"
                      << M - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else if (k >= P)
      {
        std::ostringstream error_message;
        error_message << "Range Error: k=" << k << " is not in the range (0,"
                      << P - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
 
 
  public:
    RankThreeTensor() : Tensordata(0), N(0), M(0), P(0) {}
 
    RankThreeTensor(const RankThreeTensor& source_tensor)
    {
      // Set row and column lengths
      N = source_tensor.nindex1();
      M = source_tensor.nindex2();
      P = source_tensor.nindex3();
      // Assign space for the data
      Tensordata = new T[N * M * P];
      // Copy the data across from the other matrix
      for (unsigned i = 0; i < N; i++)
      {
        for (unsigned j = 0; j < M; j++)
        {
          for (unsigned k = 0; k < P; k++)
          {
            Tensordata[P * (M * i + j) + k] = source_tensor(i, j, k);
          }
        }
      }
    }
 
    RankThreeTensor& operator=(const RankThreeTensor& source_tensor)
    {
      // Don't create a new matrix if the assignement is the identity
      if (this != &source_tensor)
      {
        // Check row and column length
        unsigned long n = source_tensor.nindex1();
        unsigned long m = source_tensor.nindex2();
        unsigned long p = source_tensor.nindex3();
        // Resie the tensor to be the same size as the old tensor
        if ((N != n) || (M != m) || (P != p))
        {
          resize(n, m, p);
        }
 
        // Copy entries across from the other matrix
        for (unsigned long i = 0; i < N; i++)
        {
          for (unsigned long j = 0; j < M; j++)
          {
            for (unsigned long k = 0; k < P; k++)
            {
              (*this)(i, j, k) = source_tensor(i, j, k);
            }
          }
        }
      }
      // Return reference to object itself (i.e. de-reference this pointer)
      return *this;
    }
 
 
    RankThreeTensor(const unsigned long& n)
    {
      // Set row and column lengths
      N = n;
      M = n;
      P = n;
      // Assign space for the n rows
      Tensordata = new T[N * M * P];
      // Initialise to zero if required
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
      initialise(T(0));
#endif
    }
 
    RankThreeTensor(const unsigned long& n_index1,
                    const unsigned long& n_index2,
                    const unsigned long& n_index3)
    {
      // Set row and column lengths
      N = n_index1;
      M = n_index2;
      P = n_index3;
      // Assign space for the n rows
      Tensordata = new T[N * M * P];
      // Initialise to zero if required
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
      initialise(T(0));
#endif
    }
 
 
    RankThreeTensor(const unsigned long& n_index1,
                    const unsigned long& n_index2,
                    const unsigned long& n_index3,
                    const T& initial_val)
    {
      // Set row and column lengths
      N = n_index1;
      M = n_index2;
      P = n_index3;
      // Assign space for the n rows
      Tensordata = new T[N * M * P];
      // Initialise to the initial value
      initialise(initial_val);
    }
 
    virtual ~RankThreeTensor()
    {
      delete[] Tensordata;
      Tensordata = 0;
    }
 
    void resize(const unsigned long& n)
    {
      resize(n, n, n);
    }
 
    void resize(const unsigned long& n_index1,
                const unsigned long& n_index2,
                const unsigned long& n_index3)
    {
      // If the sizes have not changed do nothing
      if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P))
      {
        return;
      }
      // Store old sizes
      unsigned long n_old = N, m_old = M, p_old = P;
      // Reassign the sizes
      N = n_index1;
      M = n_index2;
      P = n_index3;
      // Store triple pointer to old matrix data
      T* temp_tensor = Tensordata;
      // Re-create Tensordata in new size
      Tensordata = new T[N * M * P];
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
      initialise(T(0));
#endif
      // Transfer values
      unsigned long n_copy, m_copy, p_copy;
      n_copy = std::min(n_old, n_index1);
      m_copy = std::min(m_old, n_index2);
      p_copy = std::min(p_old, n_index3);
      // If matrix has values, transfer them to new matrix
      // Loop over rows
      for (unsigned long i = 0; i < n_copy; i++)
      {
        // Loop over columns
        for (unsigned long j = 0; j < m_copy; j++)
        {
          // Loop over columns
          for (unsigned long k = 0; k < p_copy; k++)
          {
            // Transfer values from temp_tensor
            Tensordata[M * P * i + P * j + k] =
              temp_tensor[m_old * p_old * i + p_old * j + k];
          }
        }
      }
      // Now kill storage for old tensor
      delete[] temp_tensor;
    }
 
    void resize(const unsigned long& n_index1,
                const unsigned long& n_index2,
                const unsigned long& n_index3,
                const T& initial_value)
    {
      // If the sizes have not changed do nothing
      if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P))
      {
        return;
      }
      // Store old sizes
      unsigned long n_old = N, m_old = M, p_old = P;
      // Reassign the sizes
      N = n_index1;
      M = n_index2;
      P = n_index3;
      // Store triple pointer to old matrix data
      T* temp_tensor = Tensordata;
      // Re-create Tensordata in new size
      Tensordata = new T[N * M * P];
      // Initialise the newly allocated storage
      initialise(initial_value);
 
      // Transfer values
      unsigned long n_copy, m_copy, p_copy;
      n_copy = std::min(n_old, n_index1);
      m_copy = std::min(m_old, n_index2);
      p_copy = std::min(p_old, n_index3);
      // If matrix has values, transfer them to new matrix
      // Loop over rows
      for (unsigned long i = 0; i < n_copy; i++)
      {
        // Loop over columns
        for (unsigned long j = 0; j < m_copy; j++)
        {
          // Loop over columns
          for (unsigned long k = 0; k < p_copy; k++)
          {
            // Transfer values from temp_tensor
            Tensordata[M * P * i + P * j + k] =
              temp_tensor[m_old * p_old * i + p_old * j + k];
          }
        }
      }
      // Now kill storage for old tensor
      delete[] temp_tensor;
    }
 
    void initialise(const T& val)
    {
      for (unsigned long i = 0; i < (N * M * P); ++i)
      {
        Tensordata[i] = val;
      }
    }
 
    unsigned long nindex1() const
    {
      return N;
    }
 
    unsigned long nindex2() const
    {
      return M;
    }
 
    unsigned long nindex3() const
    {
      return P;
    }
 
    inline T& operator()(const unsigned long& i,
                         const unsigned long& j,
                         const unsigned long& k)
    {
#ifdef RANGE_CHECKING
      this->range_check(i, j, k);
#endif
      return Tensordata[P * (M * i + j) + k];
    }
 
    inline T operator()(const unsigned long& i,
                        const unsigned long& j,
                        const unsigned long& k) const
    {
#ifdef RANGE_CHECKING
      this->range_check(i, j, k);
#endif
      return Tensordata[P * (M * i + j) + k];
    }
  };
 
 
 
  //=================================================================
  //=================================================================
  template<class T>
  class RankFourTensor
  {
  private:
    T* Tensordata;
 
    unsigned N;
 
    unsigned M;
 
    unsigned P;
 
    unsigned Q;
 
    void range_check(const unsigned long& i,
                     const unsigned long& j,
                     const unsigned long& k,
                     const unsigned long& l) const
    {
      if (i >= N)
      {
        std::ostringstream error_message;
        error_message << "Range Error: i=" << i << " is not in the range (0,"
                      << N - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else if (j >= M)
      {
        std::ostringstream error_message;
        error_message << "Range Error: j=" << j << " is not in the range (0,"
                      << M - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else if (k >= P)
      {
        std::ostringstream error_message;
        error_message << "Range Error: k=" << k << " is not in the range (0,"
                      << P - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else if (l >= Q)
      {
        std::ostringstream error_message;
        error_message << "Range Error: l=" << l << " is not in the range (0,"
                      << Q - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
 
  public:
    RankFourTensor() : Tensordata(0), N(0), M(0), P(0), Q(0) {}
 
    RankFourTensor(const RankFourTensor& source_tensor)
    {
      // Set row and column lengths
      N = source_tensor.nindex1();
      M = source_tensor.nindex2();
      P = source_tensor.nindex3();
      Q = source_tensor.nindex4();
 
      // Assign space for the data
      Tensordata = new T[N * M * P * Q];
 
      // Copy the data across from the other matrix
      for (unsigned i = 0; i < N; i++)
      {
        for (unsigned j = 0; j < M; j++)
        {
          for (unsigned k = 0; k < P; k++)
          {
            for (unsigned l = 0; l < Q; l++)
            {
              Tensordata[Q * (P * (M * i + j) + k) + l] =
                source_tensor(i, j, k, l);
            }
          }
        }
      }
    }
 
    RankFourTensor& operator=(const RankFourTensor& source_tensor)
    {
      // Don't create a new matrix if the assignement is the identity
      if (this != &source_tensor)
      {
        // Check row and column length
        unsigned long n = source_tensor.nindex1();
        unsigned long m = source_tensor.nindex2();
        unsigned long p = source_tensor.nindex3();
        unsigned long q = source_tensor.nindex4();
        // Resize the tensor to be the same size as the old tensor
        if ((N != n) || (M != m) || (P != p) || (Q != q))
        {
          resize(n, m, p, q);
        }
 
        // Copy entries across from the other matrix
        for (unsigned long i = 0; i < N; i++)
        {
          for (unsigned long j = 0; j < M; j++)
          {
            for (unsigned long k = 0; k < P; k++)
            {
              for (unsigned long l = 0; l < Q; l++)
              {
                (*this)(i, j, k, l) = source_tensor(i, j, k, l);
              }
            }
          }
        }
      }
      // Return reference to object itself (i.e. de-reference this pointer)
      return *this;
    }
 
 
    RankFourTensor(const unsigned long& n)
    {
      // Set row and column lengths
      N = n;
      M = n;
      P = n;
      Q = n;
      // Assign space for the n rows
      Tensordata = new T[N * M * P * Q];
      // Initialise to zero if required
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
      initialise(T(0));
#endif
    }
 
    RankFourTensor(const unsigned long& n_index1,
                   const unsigned long& n_index2,
                   const unsigned long& n_index3,
                   const unsigned long& n_index4)
    {
      // Set row and column lengths
      N = n_index1;
      M = n_index2;
      P = n_index3;
      Q = n_index4;
      // Assign space for the n rows
      Tensordata = new T[N * M * P * Q];
      // Initialise to zero if required
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
      initialise(T(0));
#endif
    }
 
 
    RankFourTensor(const unsigned long& n_index1,
                   const unsigned long& n_index2,
                   const unsigned long& n_index3,
                   const unsigned long& n_index4,
                   const T& initial_val)
    {
      // Set row and column lengths
      N = n_index1;
      M = n_index2;
      P = n_index3;
      Q = n_index4;
      // Assign space for the n rows
      Tensordata = new T[N * M * P * Q];
      // Initialise to the initial value
      initialise(initial_val);
    }
 
    virtual ~RankFourTensor()
    {
      delete[] Tensordata;
      Tensordata = 0;
    }
 
    void resize(const unsigned long& n)
    {
      resize(n, n, n, n);
    }
 
    void resize(const unsigned long& n_index1,
                const unsigned long& n_index2,
                const unsigned long& n_index3,
                const unsigned long& n_index4)
    {
      // If the sizes have not changed do nothing
      if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P) &&
          (n_index4 == Q))
      {
        return;
      }
      // Store old sizes
      unsigned long n_old = N, m_old = M, p_old = P, q_old = Q;
      // Reassign the sizes
      N = n_index1;
      M = n_index2;
      P = n_index3;
      Q = n_index4;
      // Store pointer to old matrix data
      T* temp_tensor = Tensordata;
      // Re-create Tensordata in new size
      Tensordata = new T[N * M * P * Q];
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
      initialise(T(0));
#endif
      // Transfer values
      unsigned long n_copy, m_copy, p_copy, q_copy;
      n_copy = std::min(n_old, n_index1);
      m_copy = std::min(m_old, n_index2);
      p_copy = std::min(p_old, n_index3);
      q_copy = std::min(q_old, n_index4);
      // If matrix has values, transfer them to new matrix
      // Loop over rows
      for (unsigned long i = 0; i < n_copy; i++)
      {
        // Loop over columns
        for (unsigned long j = 0; j < m_copy; j++)
        {
          // Loop over columns
          for (unsigned long k = 0; k < p_copy; k++)
          {
            // Loop over columns
            for (unsigned long l = 0; l < q_copy; l++)
            {
              // Transfer values from temp_tensor
              Tensordata[Q * (M * P * i + P * j + k) + l] =
                temp_tensor[q_old * (m_old * p_old * i + p_old * j + k) + l];
            }
          }
        }
      }
      // Now kill storage for old tensor
      delete[] temp_tensor;
    }
 
    void resize(const unsigned long& n_index1,
                const unsigned long& n_index2,
                const unsigned long& n_index3,
                const unsigned long& n_index4,
                const T& initial_value)
    {
      // If the sizes have not changed do nothing
      if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P) &&
          (n_index4 == Q))
      {
        return;
      }
      // Store old sizes
      unsigned long n_old = N, m_old = M, p_old = P, q_old = Q;
      // Reassign the sizes
      N = n_index1;
      M = n_index2;
      P = n_index3;
      Q = n_index4;
      // Store triple pointer to old matrix data
      T* temp_tensor = Tensordata;
      // Re-create Tensordata in new size
      Tensordata = new T[N * M * P * Q];
      // Initialise the newly allocated storage
      initialise(initial_value);
 
      // Transfer values
      unsigned long n_copy, m_copy, p_copy, q_copy;
      n_copy = std::min(n_old, n_index1);
      m_copy = std::min(m_old, n_index2);
      p_copy = std::min(p_old, n_index3);
      q_copy = std::min(q_old, n_index4);
      // If matrix has values, transfer them to new matrix
      // Loop over rows
      for (unsigned long i = 0; i < n_copy; i++)
      {
        // Loop over columns
        for (unsigned long j = 0; j < m_copy; j++)
        {
          // Loop over columns
          for (unsigned long k = 0; k < p_copy; k++)
          {
            // Loop over columns
            for (unsigned long l = 0; l < q_copy; l++)
            {
              // Transfer values from temp_tensor
              Tensordata[Q * (M * P * i + P * j + k) + l] =
                temp_tensor[q_old * (m_old * p_old * i + p_old * j + k) + l];
            }
          }
        }
      }
      // Now kill storage for old tensor
      delete[] temp_tensor;
    }
 
    void initialise(const T& val)
    {
      for (unsigned long i = 0; i < (N * M * P * Q); ++i)
      {
        Tensordata[i] = val;
      }
    }
 
    unsigned long nindex1() const
    {
      return N;
    }
 
    unsigned long nindex2() const
    {
      return M;
    }
 
    unsigned long nindex3() const
    {
      return P;
    }
 
    unsigned long nindex4() const
    {
      return Q;
    }
 
    inline T& operator()(const unsigned long& i,
                         const unsigned long& j,
                         const unsigned long& k,
                         const unsigned long& l)
    {
#ifdef RANGE_CHECKING
      this->range_check(i, j, k, l);
#endif
      return Tensordata[Q * (P * (M * i + j) + k) + l];
    }
 
    inline T operator()(const unsigned long& i,
                        const unsigned long& j,
                        const unsigned long& k,
                        const unsigned long& l) const
    {
#ifdef RANGE_CHECKING
      this->range_check(i, j, k, l);
#endif
      return Tensordata[Q * (P * (M * i + j) + k) + l];
    }
 
    inline T& raw_direct_access(const unsigned long& i)
    {
      return Tensordata[i];
    }
 
    inline const T& raw_direct_access(const unsigned long& i) const
    {
      return Tensordata[i];
    }
 
    unsigned offset(const unsigned long& i, const unsigned long& j) const
    {
      return (Q * (P * (M * i + j) + 0) + 0);
    }
  };
 
 
 
 
  //=================================================================
  //=================================================================
  template<class T>
  class RankFiveTensor
  {
  private:
    T* Tensordata;
 
    unsigned N;
 
    unsigned M;
 
    unsigned P;
 
    unsigned Q;
 
    unsigned R;
 
    void range_check(const unsigned long& i,
                     const unsigned long& j,
                     const unsigned long& k,
                     const unsigned long& l,
                     const unsigned long& m) const
    {
      if (i >= N)
      {
        std::ostringstream error_message;
        error_message << "Range Error: i=" << i << " is not in the range (0,"
                      << N - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else if (j >= M)
      {
        std::ostringstream error_message;
        error_message << "Range Error: j=" << j << " is not in the range (0,"
                      << M - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else if (k >= P)
      {
        std::ostringstream error_message;
        error_message << "Range Error: k=" << k << " is not in the range (0,"
                      << P - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else if (l >= Q)
      {
        std::ostringstream error_message;
        error_message << "Range Error: l=" << l << " is not in the range (0,"
                      << Q - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
      else if (m >= R)
      {
        std::ostringstream error_message;
        error_message << "Range Error: m=" << m << " is not in the range (0,"
                      << R - 1 << ")." << std::endl;
 
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
    }
 
  public:
    RankFiveTensor() : Tensordata(0), N(0), M(0), P(0), Q(0), R(0) {}
 
    RankFiveTensor(const RankFiveTensor& source_tensor)
    {
      // Set row and column lengths
      N = source_tensor.nindex1();
      M = source_tensor.nindex2();
      P = source_tensor.nindex3();
      Q = source_tensor.nindex4();
      R = source_tensor.nindex5();
 
      // Assign space for the data
      Tensordata = new T[N * M * P * Q * R];
 
      // Copy the data across from the other matrix
      for (unsigned i = 0; i < N; i++)
      {
        for (unsigned j = 0; j < M; j++)
        {
          for (unsigned k = 0; k < P; k++)
          {
            for (unsigned l = 0; l < Q; l++)
            {
              for (unsigned m = 0; m < R; m++)
              {
                Tensordata[R * (Q * (P * (M * i + j) + k) + l) + m] =
                  source_tensor(i, j, k, l, m);
              }
            }
          }
        }
      }
    }
 
    RankFiveTensor& operator=(const RankFiveTensor& source_tensor)
    {
      // Don't create a new matrix if the assignement is the identity
      if (this != &source_tensor)
      {
        // Check row and column length
        unsigned long n = source_tensor.nindex1();
        unsigned long m = source_tensor.nindex2();
        unsigned long p = source_tensor.nindex3();
        unsigned long q = source_tensor.nindex4();
        unsigned long r = source_tensor.nindex5();
        // Resize the tensor to be the same size as the old tensor
        if ((N != n) || (M != m) || (P != p) || (Q != q) || (R != r))
        {
          resize(n, m, p, q, r);
        }
 
        // Copy entries across from the other matrix
        for (unsigned long i = 0; i < N; i++)
        {
          for (unsigned long j = 0; j < M; j++)
          {
            for (unsigned long k = 0; k < P; k++)
            {
              for (unsigned long l = 0; l < Q; l++)
              {
                for (unsigned long m = 0; m < R; m++)
                {
                  (*this)(i, j, k, l, m) = source_tensor(i, j, k, l, m);
                }
              }
            }
          }
        }
      }
      // Return reference to object itself (i.e. de-reference this pointer)
      return *this;
    }
 
 
    RankFiveTensor(const unsigned long& n)
    {
      // Set row and column lengths
      N = n;
      M = n;
      P = n;
      Q = n;
      R = n;
      // Assign space for the n rows
      Tensordata = new T[N * M * P * Q * R];
      // Initialise to zero if required
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
      initialise(T(0));
#endif
    }
 
    RankFiveTensor(const unsigned long& n_index1,
                   const unsigned long& n_index2,
                   const unsigned long& n_index3,
                   const unsigned long& n_index4,
                   const unsigned long& n_index5)
    {
      // Set row and column lengths
      N = n_index1;
      M = n_index2;
      P = n_index3;
      Q = n_index4;
      R = n_index5;
      // Assign space for the n rows
      Tensordata = new T[N * M * P * Q * R];
      // Initialise to zero if required
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
      initialise(T(0));
#endif
    }
 
 
    RankFiveTensor(const unsigned long& n_index1,
                   const unsigned long& n_index2,
                   const unsigned long& n_index3,
                   const unsigned long& n_index4,
                   const unsigned long& n_index5,
                   const T& initial_val)
    {
      // Set row and column lengths
      N = n_index1;
      M = n_index2;
      P = n_index3;
      Q = n_index4;
      R = n_index5;
      // Assign space for the n rows
      Tensordata = new T[N * M * P * Q * R];
      // Initialise to the initial value
      initialise(initial_val);
    }
 
    virtual ~RankFiveTensor()
    {
      delete[] Tensordata;
      Tensordata = 0;
    }
 
    void resize(const unsigned long& n)
    {
      resize(n, n, n, n, n);
    }
 
    void resize(const unsigned long& n_index1,
                const unsigned long& n_index2,
                const unsigned long& n_index3,
                const unsigned long& n_index4,
                const unsigned long& n_index5)
    {
      // If the sizes have not changed do nothing
      if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P) &&
          (n_index4 == Q) && (n_index5 == R))
      {
        return;
      }
      // Store old sizes
      unsigned long n_old = N, m_old = M, p_old = P, q_old = Q, r_old = R;
      // Reassign the sizes
      N = n_index1;
      M = n_index2;
      P = n_index3;
      Q = n_index4;
      R = n_index5;
      // Store pointer to old matrix data
      T* temp_tensor = Tensordata;
      // Re-create Tensordata in new size
      Tensordata = new T[N * M * P * Q * R];
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
      initialise(T(0));
#endif
      // Transfer values
      unsigned long n_copy, m_copy, p_copy, q_copy, r_copy;
      n_copy = std::min(n_old, n_index1);
      m_copy = std::min(m_old, n_index2);
      p_copy = std::min(p_old, n_index3);
      q_copy = std::min(q_old, n_index4);
      r_copy = std::min(r_old, n_index5);
      // If matrix has values, transfer them to new matrix
      // Loop over rows
      for (unsigned long i = 0; i < n_copy; i++)
      {
        // Loop over columns
        for (unsigned long j = 0; j < m_copy; j++)
        {
          // Loop over columns
          for (unsigned long k = 0; k < p_copy; k++)
          {
            // Loop over columns
            for (unsigned long l = 0; l < q_copy; l++)
            {
              // Loop over columns
              for (unsigned long m = 0; m < r_copy; m++)
              {
                // Transfer values from temp_tensor
                Tensordata[R * (Q * (M * P * i + P * j + k) + l) + m] =
                  temp_tensor[r_old *
                                (q_old * (m_old * p_old * i + p_old * j + k) +
                                 l) +
                              m];
              }
            }
          }
        }
      }
      // Now kill storage for old tensor
      delete[] temp_tensor;
    }
 
    void resize(const unsigned long& n_index1,
                const unsigned long& n_index2,
                const unsigned long& n_index3,
                const unsigned long& n_index4,
                const unsigned long& n_index5,
                const T& initial_value)
    {
      // If the sizes have not changed do nothing
      if ((n_index1 == N) && (n_index2 == M) && (n_index3 == P) &&
          (n_index4 == Q) && (n_index5 == R))
      {
        return;
      }
      // Store old sizes
      unsigned long n_old = N, m_old = M, p_old = P, q_old = Q, r_old = R;
      // Reassign the sizes
      N = n_index1;
      M = n_index2;
      P = n_index3;
      Q = n_index4;
      R = n_index5;
      // Store triple pointer to old matrix data
      T* temp_tensor = Tensordata;
      // Re-create Tensordata in new size
      Tensordata = new T[N * M * P * Q * R];
      // Initialise the newly allocated storage
      initialise(initial_value);
 
      // Transfer values
      unsigned long n_copy, m_copy, p_copy, q_copy, r_copy;
      n_copy = std::min(n_old, n_index1);
      m_copy = std::min(m_old, n_index2);
      p_copy = std::min(p_old, n_index3);
      q_copy = std::min(q_old, n_index4);
      r_copy = std::min(r_old, n_index5);
      // If matrix has values, transfer them to new matrix
      // Loop over rows
      for (unsigned long i = 0; i < n_copy; i++)
      {
        // Loop over columns
        for (unsigned long j = 0; j < m_copy; j++)
        {
          // Loop over columns
          for (unsigned long k = 0; k < p_copy; k++)
          {
            // Loop over columns
            for (unsigned long l = 0; l < q_copy; l++)
            {
              // Loop over columns
              for (unsigned long m = 0; m < r_copy; m++)
              {
                // Transfer values from temp_tensor
                Tensordata[R * (Q * (M * P * i + P * j + k) + l) + m] =
                  temp_tensor[r_old *
                                (q_old * (m_old * p_old * i + p_old * j + k) +
                                 l) +
                              m];
              }
            }
          }
        }
      }
      // Now kill storage for old tensor
      delete[] temp_tensor;
    }
 
    void initialise(const T& val)
    {
      for (unsigned long i = 0; i < (N * M * P * Q * R); ++i)
      {
        Tensordata[i] = val;
      }
    }
 
    unsigned long nindex1() const
    {
      return N;
    }
 
    unsigned long nindex2() const
    {
      return M;
    }
 
    unsigned long nindex3() const
    {
      return P;
    }
 
    unsigned long nindex4() const
    {
      return Q;
    }
 
    unsigned long nindex5() const
    {
      return R;
    }
 
    inline T& operator()(const unsigned long& i,
                         const unsigned long& j,
                         const unsigned long& k,
                         const unsigned long& l,
                         const unsigned long& m)
    {
#ifdef RANGE_CHECKING
      this->range_check(i, j, k, l, m);
#endif
      return Tensordata[R * (Q * (P * (M * i + j) + k) + l) + m];
    }
 
    inline T operator()(const unsigned long& i,
                        const unsigned long& j,
                        const unsigned long& k,
                        const unsigned long& l,
                        const unsigned long& m) const
    {
#ifdef RANGE_CHECKING
      this->range_check(i, j, k, l, m);
#endif
      return Tensordata[R * (Q * (P * (M * i + j) + k) + l) + m];
    }
 
    inline T& raw_direct_access(const unsigned long& i)
    {
      return Tensordata[i];
    }
 
 
    inline const T& raw_direct_access(const unsigned long& i) const
    {
      return Tensordata[i];
    }
 
    unsigned offset(const unsigned long& i,
                    const unsigned long& j,
                    const unsigned long& k) const
    {
      return (R * (Q * (P * (M * i + j) + k) + 0) + 0);
    }
  };
 
 
 
  //=======================================================================
  //=======================================================================
  template<class T>
  class CCMatrix : public SparseMatrix<T, CCMatrix<T>>
  {
  public:
    CCMatrix() : SparseMatrix<T, CCMatrix<T>>()
    {
      Row_index = 0;
      Column_start = 0;
    }
 
 
    CCMatrix(const Vector<T>& value,
             const Vector<int>& row_index_,
             const Vector<int>& column_start_,
             const unsigned long& n,
             const unsigned long& m)
      : SparseMatrix<T, CCMatrix<T>>()
    {
      Row_index = 0;
      Column_start = 0;
      build(value, row_index_, column_start_, n, m);
    }
 
 
    CCMatrix(const CCMatrix& source_matrix)
      : SparseMatrix<T, CCMatrix<T>>(source_matrix)
    {
      // NNz, N and M are set the the copy constructor of the SparseMatrix
      // called above
 
      // Row indices stored in C-style array
      Row_index = new int[this->Nnz];
 
      // Assign:
      for (unsigned long i = 0; i < this->Nnz; i++)
      {
        Row_index[i] = source_matrix.row_index()[i];
      }
 
      // Column start:
      Column_start = new int[this->M + 1];
 
      // Assign:
      for (unsigned long i = 0; i <= this->M; i++)
      {
        Column_start[i] = source_matrix.column_start()[i];
      }
    }
 
    void operator=(const CCMatrix&) = delete;
 
 
    virtual ~CCMatrix()
    {
      delete[] Row_index;
      Row_index = 0;
      delete[] Column_start;
      Column_start = 0;
    }
 
    T get_entry(const unsigned long& i, const unsigned long& j) const
    {
#ifdef RANGE_CHECKING
      this->range_check(i, j);
#endif
      for (long k = Column_start[j]; k < Column_start[j + 1]; k++)
      {
        if (unsigned(Row_index[k]) == i)
        {
          return this->Value[k];
        }
      }
      return this->Zero;
    }
 
    T& entry(const unsigned long& i, const unsigned long& j)
    {
      std::string error_string =
        "Non-const access not provided for the CCMatrix<T> class\n";
      error_string +=
        "It is not possible to use round-bracket access: M(i,j)\n";
      error_string += "if M is not declared as const.\n";
      error_string += "The solution (albeit ugly) is to create const reference "
                      "to the matrix\n";
      error_string += " const CCMatrix<T>& read_M = M;\n";
      error_string += "Then read_M(i,j) is permitted\n";
 
      throw OomphLibError(
        error_string, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
 
      // Dummy return
      T dummy;
      return dummy;
    }
 
    int* column_start()
    {
      return Column_start;
    }
 
    const int* column_start() const
    {
      return Column_start;
    }
 
    int* row_index()
    {
      return Row_index;
    }
 
    const int* row_index() const
    {
      return Row_index;
    }
 
    void output_bottom_right_zero_helper(std::ostream& outfile) const
    {
      int last_row = this->N - 1;
      int last_col_local = this->M - 1;
 
      // Use this strange thingy because of the CRTP discussed above.
      T last_value = this->operator()(last_row, last_col_local);
 
      if (last_value == T(0))
      {
        outfile << last_row << " " << last_col_local << " " << T(0)
                << std::endl;
      }
    }
 
    void sparse_indexed_output_helper(std::ostream& outfile) const
    {
      for (unsigned long j = 0; j < this->N; j++)
      {
        for (long k = Column_start[j]; k < Column_start[j + 1]; k++)
        {
          outfile << Row_index[k] << " " << j << " " << this->Value[k]
                  << std::endl;
        }
      }
    }
 
    void clean_up_memory();
 
 
    void build(const Vector<T>& value,
               const Vector<int>& row_index,
               const Vector<int>& column_start,
               const unsigned long& n,
               const unsigned long& m);
 
    void build_without_copy(T* value,
                            int* row_index,
                            int* column_start,
                            const unsigned long& nnz,
                            const unsigned long& n,
                            const unsigned long& m);
 
 
  protected:
    int* Row_index;
 
    int* Column_start;
  };
 
 
 
  //=================================================================
  //=================================================================
  class CCDoubleMatrix : public DoubleMatrixBase, public CCMatrix<double>
  {
  public:
    CCDoubleMatrix();
 
    CCDoubleMatrix(const Vector<double>& value,
                   const Vector<int>& row_index_,
                   const Vector<int>& column_start_,
                   const unsigned long& n,
                   const unsigned long& m);
 
    CCDoubleMatrix(const CCDoubleMatrix& matrix) = delete;
 
    void operator=(const CCDoubleMatrix&) = delete;
 
    virtual ~CCDoubleMatrix();
 
    inline unsigned long nrow() const
    {
      return CCMatrix<double>::nrow();
    }
 
    inline unsigned long ncol() const
    {
      return CCMatrix<double>::ncol();
    }
 
    inline double operator()(const unsigned long& i,
                             const unsigned long& j) const
    {
      return CCMatrix<double>::get_entry(i, j);
    }
 
    virtual void ludecompose();
 
    virtual void lubksub(DoubleVector& rhs);
 
    void multiply(const DoubleVector& x, DoubleVector& soln) const;
 
    void multiply_transpose(const DoubleVector& x, DoubleVector& soln) const;
 
 
    void multiply(const CCDoubleMatrix& matrix_in, CCDoubleMatrix& result);
 
 
    void matrix_reduction(const double& alpha, CCDoubleMatrix& reduced_matrix);
 
    unsigned& matrix_matrix_multiply_method()
    {
      return Matrix_matrix_multiply_method;
    }
 
  private:
    unsigned Matrix_matrix_multiply_method;
  };
 
 
 
 
  //============================================================================
  //============================================================================
  template<class T>
  DenseMatrix<T>::DenseMatrix(const unsigned long& n)
  {
    // Set row and column lengths
    N = n;
    M = n;
    // Assign space for the n rows
    Matrixdata = new T[n * n];
    // Initialise to zero if required
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
    initialise(T(0));
#endif
  }
 
 
  //============================================================================
  //============================================================================
  template<class T>
  DenseMatrix<T>::DenseMatrix(const unsigned long& n, const unsigned long& m)
  {
    // Set row and column lengths
    N = n;
    M = m;
    // Assign space for the n rows
    Matrixdata = new T[n * m];
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
    initialise(T(0));
#endif
  }
 
  //============================================================================
  //============================================================================
  template<class T>
  DenseMatrix<T>::DenseMatrix(const unsigned long& n,
                              const unsigned long& m,
                              const T& initial_val)
  {
    // Set row and column lengths
    N = n;
    M = m;
    // Assign space for the n rows
    Matrixdata = new T[n * m];
    initialise(initial_val);
  }
 
 
  //============================================================================
  //============================================================================
  template<class T>
  void DenseMatrix<T>::resize(const unsigned long& n, const unsigned long& m)
  {
    // If the sizes are the same, do nothing
    if ((n == N) && (m == M))
    {
      return;
    }
    // Store old sizes
    unsigned long n_old = N, m_old = M;
    // Reassign the sizes
    N = n;
    M = m;
    // Store double pointer to old matrix data
    T* temp_matrix = Matrixdata;
 
    // Re-create Matrixdata in new size
    Matrixdata = new T[n * m];
    // Initialise to zero
#ifdef OOMPH_INITIALISE_DENSE_MATRICES
    initialise(T(0));
#endif
 
    // Transfer previously existing values
    unsigned long n_copy, m_copy;
    n_copy = std::min(n_old, n);
    m_copy = std::min(m_old, m);
 
    // If matrix has values, transfer them to new matrix
    // Loop over rows
    for (unsigned long i = 0; i < n_copy; i++)
    {
      // Loop over columns
      for (unsigned long j = 0; j < m_copy; j++)
      {
        // Transfer values from temp_matrix
        Matrixdata[m * i + j] = temp_matrix[m_old * i + j];
      }
    }
 
    // Now kill storage for old matrix
    delete[] temp_matrix;
  }
 
 
  //============================================================================
  //============================================================================
  template<class T>
  void DenseMatrix<T>::resize(const unsigned long& n,
                              const unsigned long& m,
                              const T& initial_value)
  {
    // If the size is not changed, just return
    if ((n == N) && (m == M))
    {
      return;
    }
    // Store old sizes
    unsigned long n_old = N, m_old = M;
    // Reassign the sizes
    N = n;
    M = m;
    // Store double pointer to old matrix data
    T* temp_matrix = Matrixdata;
    // Re-create Matrixdata in new size
    Matrixdata = new T[n * m];
    // Assign initial value (will use the newly allocated data)
    initialise(initial_value);
 
    // Transfering values
    unsigned long n_copy, m_copy;
    n_copy = std::min(n_old, n);
    m_copy = std::min(m_old, m);
    // If matrix has values, transfer them to temp_matrix
    // Loop over rows
    for (unsigned long i = 0; i < n_copy; i++)
    {
      // Loop over columns
      for (unsigned long j = 0; j < m_copy; j++)
      {
        // Transfer values to temp_matrix
        Matrixdata[m * i + j] = temp_matrix[m_old * i + j];
      }
    }
 
    // Now kill storage for old matrix
    delete[] temp_matrix;
  }
 
 
  //============================================================================
  //============================================================================
  template<class T>
  void DenseMatrix<T>::output(std::ostream& outfile) const
  {
    // Loop over the rows
    for (unsigned i = 0; i < N; i++)
    {
      // Loop over the columne
      for (unsigned j = 0; j < M; j++)
      {
        outfile << (*this)(i, j) << " ";
      }
      // Put in a newline
      outfile << std::endl;
    }
  }
 
 
  //============================================================================
  //============================================================================
  template<class T>
  void DenseMatrix<T>::output(std::string filename) const
  {
    // Open file
    std::ofstream some_file;
    some_file.open(filename.c_str());
 
    output(some_file);
    some_file.close();
  }
 
 
  //============================================================================
  //============================================================================
  template<class T>
  void DenseMatrix<T>::indexed_output(std::ostream& outfile) const
  {
    // Loop over the rows
    for (unsigned i = 0; i < N; i++)
    {
      // Loop over the columns
      for (unsigned j = 0; j < M; j++)
      {
        outfile << i << " " << j << " " << (*this)(i, j) << std::endl;
      }
    }
  }
 
 
  //============================================================================
  //============================================================================
  template<class T>
  void DenseMatrix<T>::indexed_output(std::string filename) const
  {
    // Open file
    std::ofstream some_file;
    some_file.open(filename.c_str());
    indexed_output(some_file);
    some_file.close();
  }
 
 
  //============================================================================
  //============================================================================
  template<class T>
  void DenseMatrix<T>::output_bottom_right_zero_helper(
    std::ostream& outfile) const
  {
    int last_row = this->N - 1;
    int last_col = this->M - 1;
 
    // Use this strange thingy because of the CRTP discussed above.
    T last_value = this->operator()(last_row, last_col);
 
    if (last_value == T(0))
    {
      outfile << last_row << " " << last_col << " " << T(0) << std::endl;
    }
  }
 
  //============================================================================
  //============================================================================
  template<class T>
  void DenseMatrix<T>::sparse_indexed_output_helper(std::ostream& outfile) const
  {
    // Loop over the rows
    for (unsigned i = 0; i < N; i++)
    {
      // Loop over the column
      for (unsigned j = 0; j < M; j++)
      {
        if ((*this)(i, j) != T(0))
        {
          outfile << i << " " << j << " " << (*this)(i, j) << std::endl;
        }
      }
    }
  }
 
 
 
 
  //=============================================================================
  //=============================================================================
  template<class T>
  void CCMatrix<T>::clean_up_memory()
  {
    // delete any previously allocated  storage
    if (this->Value != 0)
    {
      delete[] this->Value;
      this->Value = 0;
    }
    if (this->Row_index != 0)
    {
      delete[] this->Row_index;
      this->Row_index = 0;
    }
    if (this->Column_start != 0)
    {
      delete[] this->Column_start;
      this->Column_start = 0;
    }
    this->Nnz = 0;
    this->N = 0;
    this->M = 0;
  }
 
 
  //=============================================================================
  //=============================================================================
  template<class T>
  void CCMatrix<T>::build_without_copy(T* value,
                                       int* row_index,
                                       int* column_start,
                                       const unsigned long& nnz,
                                       const unsigned long& n,
                                       const unsigned long& m)
  {
    // Number of nonzero entries
    this->Nnz = nnz;
 
    // Number of rows
    this->N = n;
 
    // Number of columns
    this->M = m;
 
    // delete any previously allocated  storage
    if (this->Value != 0)
    {
      delete[] this->Value;
    }
    if (this->Row_index != 0)
    {
      delete[] this->Row_index;
    }
    if (this->Column_start != 0)
    {
      delete[] this->Column_start;
    }
 
    // set Value
    this->Value = value;
 
    // set Row_index
    this->Row_index = row_index;
 
    // set Column_start
    this->Column_start = column_start;
  }
 
 
  //===================================================================
  //===================================================================
  template<class T>
  void CCMatrix<T>::build(const Vector<T>& value,
                          const Vector<int>& row_index_,
                          const Vector<int>& column_start_,
                          const unsigned long& n,
                          const unsigned long& m)
  {
#ifdef PARANOID
    if (value.size() != row_index_.size())
    {
      std::ostringstream error_message;
      error_message << "length of value " << value.size()
                    << " and row_index vectors " << row_index_.size()
                    << " should match " << std::endl;
 
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
 
    // Number of nonzero entries
    this->Nnz = value.size();
 
    // Number of rows
    this->N = n;
 
    // Number of columns
    this->M = m;
 
    // We need to delete any previously allocated  storage
    if (this->Value != 0)
    {
      delete[] this->Value;
    }
    if (this->Row_index != 0)
    {
      delete[] this->Row_index;
    }
    if (this->Column_start != 0)
    {
      delete[] this->Column_start;
    }
 
    // Values stored in C-style array
    this->Value = new T[this->Nnz];
 
    // Row indices stored in C-style array
    this->Row_index = new int[this->Nnz];
 
    // Assign:
    for (unsigned long i = 0; i < this->Nnz; i++)
    {
      this->Value[i] = value[i];
      this->Row_index[i] = row_index_[i];
    }
 
    // Column start:
    // Find the size and aollcate
    unsigned long n_column_start = column_start_.size();
    this->Column_start = new int[n_column_start];
 
    // Assign:
    for (unsigned long i = 0; i < n_column_start; i++)
    {
      this->Column_start[i] = column_start_[i];
    }
  }
 
 
 
  //=============================================================================
  //=============================================================================
  template<class T>
  void CRMatrix<T>::clean_up_memory()
  {
    // delete any previously allocated  storage
    if (this->Value != 0)
    {
      delete[] this->Value;
      this->Value = 0;
    }
    if (this->Column_index != 0)
    {
      delete[] this->Column_index;
      this->Column_index = 0;
    }
    if (this->Row_start != 0)
    {
      delete[] this->Row_start;
      this->Row_start = 0;
    }
    this->Nnz = 0;
    this->N = 0;
    this->M = 0;
  }
 
 
  //=============================================================================
  //=============================================================================
  template<class T>
  void CRMatrix<T>::build_without_copy(T* value,
                                       int* column_index_,
                                       int* row_start_,
                                       const unsigned long& nnz,
                                       const unsigned long& n,
                                       const unsigned long& m)
  {
    // Number of nonzero entries
    this->Nnz = nnz;
 
    // Number of rows
    this->N = n;
 
    // Number of columns
    this->M = m;
 
    // delete any previously allocated  storage
    if (this->Value != 0)
    {
      delete[] this->Value;
    }
    if (this->Column_index != 0)
    {
      delete[] this->Column_index;
    }
    if (this->Row_start != 0)
    {
      delete[] this->Row_start;
    }
 
    // set Value
    this->Value = value;
 
    // set Column_index
    this->Column_index = column_index_;
 
    // set Row_start
    this->Row_start = row_start_;
  }
 
 
  //=================================================================
  //=================================================================
  template<class T>
  void CRMatrix<T>::build(const Vector<T>& value,
                          const Vector<int>& column_index_,
                          const Vector<int>& row_start_,
                          const unsigned long& n,
                          const unsigned long& m)
  {
#ifdef PARANOID
    if (value.size() != column_index_.size())
    {
      std::ostringstream error_message;
      error_message << "Must have the same number of values and column indices,"
                    << "we have " << value.size() << " values and "
                    << column_index_.size() << " column inidicies."
                    << std::endl;
      throw OomphLibError(
        error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
#endif
    // Number of nonzero entries
    this->Nnz = value.size();
 
    // Number of rows
    this->N = n;
 
    // Number of columns
    this->M = m;
 
    // We need to delete any previously allocated  storage
    if (this->Value != 0)
    {
      delete[] this->Value;
    }
    if (this->Column_index != 0)
    {
      delete[] this->Column_index;
    }
    if (this->Row_start != 0)
    {
      delete[] this->Row_start;
    }
 
    // Values stored in C-style array
    this->Value = new T[this->Nnz];
 
    // Column indices stored in C-style array
    this->Column_index = new int[this->Nnz];
 
    // Assign:
    for (unsigned long i = 0; i < this->Nnz; i++)
    {
      this->Value[i] = value[i];
      this->Column_index[i] = column_index_[i];
    }
 
    // Row start:
    // Find the size and allocate
    unsigned long n_row_start = row_start_.size();
    this->Row_start = new int[n_row_start];
 
    // Assign:
    for (unsigned long i = 0; i < n_row_start; i++)
    {
      this->Row_start[i] = row_start_[i];
    }
  }
 
 
  //=================================================================
  //=================================================================
  template<class T, class MATRIX_TYPE>
  T SparseMatrix<T, MATRIX_TYPE>::Zero = T(0);
 
 
  namespace RRR
  {
    extern std::string RayStr;
    extern bool RayBool;
  } // namespace RRR
 
  //=================================================================
  //=================================================================
  namespace CRDoubleMatrixHelpers
  {
    inline void deep_copy(const CRDoubleMatrix* const in_matrix_pt,
                          CRDoubleMatrix& out_matrix)
    {
#ifdef PARANOID
      // Is the out matrix built? We need an empty out matrix!
      if (out_matrix.built())
      {
        std::ostringstream err_msg;
        err_msg << "The result matrix has been built.\n"
                << "Please clear the matrix.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Check that the in matrix pointer is not null.
      if (in_matrix_pt == 0)
      {
        std::ostringstream err_msg;
        err_msg << "The in_matrix_pt is null.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
 
      // Check that the in matrix is built.
      if (!in_matrix_pt->built())
      {
        std::ostringstream err_msg;
        err_msg << "The in_matrix_pt is null.\n";
        throw OomphLibError(
          err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // First set the matrix matrix multiply methods (for both serial and
      // distributed)
      out_matrix.serial_matrix_matrix_multiply_method() =
        in_matrix_pt->serial_matrix_matrix_multiply_method();
 
      out_matrix.distributed_matrix_matrix_multiply_method() =
        in_matrix_pt->distributed_matrix_matrix_multiply_method();
 
 
      // The local nrow and nnz of the in matrix
      const unsigned in_nrow_local = in_matrix_pt->nrow_local();
      const unsigned long in_nnz = in_matrix_pt->nnz();
 
      // Storage for the values, column indices and row start
      double* out_values = new double[in_nnz];
      int* out_column_indices = new int[in_nnz];
      int* out_row_start = new int[in_nrow_local + 1];
 
      // The data to copy over
      const double* const in_values = in_matrix_pt->value();
      const int* const in_column_indices = in_matrix_pt->column_index();
      const int* const in_row_start = in_matrix_pt->row_start();
 
      // Copy the data
      std::copy(in_values, in_values + in_nnz, out_values);
 
      std::copy(
        in_column_indices, in_column_indices + in_nnz, out_column_indices);
 
      std::copy(
        in_row_start, in_row_start + (in_nrow_local + 1), out_row_start);
 
      // Build the matrix
      out_matrix.build(in_matrix_pt->distribution_pt());
 
      out_matrix.build_without_copy(in_matrix_pt->ncol(),
                                    in_nnz,
                                    out_values,
                                    out_column_indices,
                                    out_row_start);
 
      // The only thing we haven't copied over is the default linear solver
      // pointer, but I cannot figure out how to copy over a solver since
      // I do not know what it is.
    } // EoFunc deep_copy
 
    void create_uniformly_distributed_matrix(
      const unsigned& nrow,
      const unsigned& ncol,
      const OomphCommunicator* const comm_pt,
      const Vector<double>& values,
      const Vector<int>& column_indicies,
      const Vector<int>& row_start,
      CRDoubleMatrix& mat_out);
 
 
    double inf_norm(const DenseMatrix<CRDoubleMatrix*>& matrix_pt);
 
    double gershgorin_eigenvalue_estimate(
      const DenseMatrix<CRDoubleMatrix*>& matrix_pt);
 
    void concatenate(const DenseMatrix<CRDoubleMatrix*>& matrix_pt,
                     CRDoubleMatrix& result_matrix);
 
    void concatenate_without_communication(
      const Vector<LinearAlgebraDistribution*>& row_distribution_pt,
      const Vector<LinearAlgebraDistribution*>& col_distribution_pt,
      const DenseMatrix<CRDoubleMatrix*>& matrix_pt,
      CRDoubleMatrix& result_matrix);
 
    void concatenate_without_communication(
      const Vector<LinearAlgebraDistribution*>& block_distribution_pt,
      const DenseMatrix<CRDoubleMatrix*>& matrix_pt,
      CRDoubleMatrix& result_matrix);
 
  } // namespace CRDoubleMatrixHelpers
 
} // namespace oomph
#endif