frontal_solver.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//====================================================================
//Interface to HSL frontal solver (fortran)


#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif

//oomph-lib headers
#include "cfortran.h"
#include "frontal.h"
#include "frontal_solver.h"
#include "problem.h"
#include "mesh.h"
#include "matrices.h"
#include "assembly_handler.h"

namespace oomph
{


//====================================================================
/// Special solver for problems with one DOF -- HSL_MA42 can't handle 
/// that!
//====================================================================
void HSL_MA42::solve_for_one_dof(Problem* const &problem_pt, 
                                 DoubleVector &result)
{

 //Find the number of elements
 unsigned long n_el = problem_pt->mesh_pt()->nelement();

#ifdef PARANOID
 //Find the number of dofs (variables)
 int n_dof = problem_pt->ndof();
 // Check
 if (n_dof!=1)
  {
   throw OomphLibError(
    "You can only call this if the problem has just one dof!",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
#endif

 // Global jac and residuals are scalars!
 double global_jac=0.0;
 double global_res=0.0;

 //Locally cache pointer to assembly handler
 AssemblyHandler* const assembly_handler_pt = 
  problem_pt->assembly_handler_pt();

 
 // Do the assembly
 for(unsigned long e=0;e<n_el;e++)
  {
   //Get pointer to the element
   GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);

   //Get number of variables in element
   int nvar = assembly_handler_pt->ndof(elem_pt);
 
   // Only assemble if there's something to be assembled
   if (nvar>0)
    {
     //Set up matrices and Vector for call to get_residuals
     Vector<double> residuals(nvar);
     DenseMatrix<double> jacobian(nvar);

     //Get the residuals and jacobian
     assembly_handler_pt->get_jacobian(elem_pt,residuals,jacobian);
   
     // Add contribution
     global_jac+=jacobian(0,0);
     global_res+=residuals[0];
    }
  }
 
 // "Solve"
 result[0]=global_res/global_jac;
 
 //If we are storing the matrix, allocate the storage
 if(Enable_resolve)
  {
   //If storage has been allocated delete it
   if(W!=0) {delete[] W;}
   //Now allocate new storage
   W = new double[1];
   //Store the global jacobian
   W[0] = global_jac;
   //Set the number of degrees of freedom
   N_dof=1;
  }
 //Set the sign of the jacobian
 problem_pt->sign_of_jacobian() = 
  static_cast<int>(std::fabs(global_jac)/global_jac);
}

//====================================================================
/// Wrapper for HSL MA42 frontal solver
//====================================================================
void HSL_MA42::solve(Problem* const &problem_pt, DoubleVector &result)
{
#ifdef OOMPH_HAS_MPI
 if(problem_pt->communicator_pt()->nproc() > 1)
  {
   std::ostringstream error_message;
   error_message 
    << "HSL_MA42 solver cannot be used in parallel.\n"
    << "Please change to another linear solver.\n"
    << "If you want to use a frontal solver then try MumpsSolver\n";
   
   throw OomphLibError(error_message.str(),
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }
#endif
    

 //Find the number of elements
 unsigned long n_el = problem_pt->mesh_pt()->nelement();

 //Find the number of dofs (variables) and store for resolves
 N_dof = problem_pt->ndof();

 //Build the distribution .. this is a serial solver
 LinearAlgebraDistribution dist(problem_pt->communicator_pt(),
                                N_dof,false); 
 this->build_distribution(dist);

 //Cast the number of dofs into an integer for the HSL solver
 int n_dof = N_dof;

 // Special case: Just one variable! MA42 can't handle that
 if (n_dof==1) 
  {
   solve_for_one_dof(problem_pt,result);
   return;
  }

 // Reorder?
 if (this->Reorder_flag)
  {
   reorder_elements(problem_pt);
  }

 //Set the control flags
 int ifsize[5];
 double cntl[2], rinfo[2];
 
 //Set up memory for last and for ndf
 int ndf, nmaxe=2, nrhs=1, lrhs=1;

 //Memory for last and dx should be allocated dynamically from the heap
 //rather than the stack, otherwise one gets a nasty segmentation fault
 int *last = new int[n_dof];
 //Set up memory for dx
 double **dx = new double *;
 *dx = new double[n_dof];

 // Provide storage for exact values required for lenbuf array
 int exact_lenbuf[3]={0,0,0};
 bool exact_lenbuf_available=false;
 
 // Storage for recommendations for lenbuf so we can check how
 // close our factors are
 int lenbuf0_recommendation=0;
 int lenbuf1_recommendation=0;
 int lenbuf2_recommendation=0;

 // Provide storage
 int lenbuf[3];


 // Provide storage for exact values required for lenfle array
 int exact_lenfle[3]={0,0,0};
 bool exact_lenfle_available=false;

 // Storage for recommendation for lenfle so we can check how
 // close our factor is
 int lenfle0_recommendation=0;
 int lenfle1_recommendation=0;
 int lenfle2_recommendation=0;

 // Provide storage
 int lenfle[3];


 // Storage for recommendations for front size so we can check how
 // close our factors are
 double front0_recommendation=0;
 double front1_recommendation=0;


 // Flag to indicate if exact, required values for nfront are available
 // from previous unsucessful solve
 bool exact_nfront_available=false;

 // Storage for front sizes
 int nfront[2];
 
 //Locally cache pointer to assembly handler
 AssemblyHandler* const assembly_handler_pt = 
  problem_pt->assembly_handler_pt();


 // Loop over size allocation/solver until we've made all the arrays
 //=================================================================
 // big enough...
 //==============
 bool not_done=true;
 while (not_done)
  {

   //Initialise frontal solver
   //=========================
   MA42ID(Icntl,cntl,Isave);
   
   
   //Loop over the elements to setup variables associated with elements
   //==================================================================
   for(unsigned long e=0;e<n_el;e++)
    {
     //Pointer to the element
     GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);

     // MH: changed array to allocateable
     int* ivar;
     
     //Get number of variables in element
     int nvar = assembly_handler_pt->ndof(elem_pt);
     
     
     // Multiple equations
     //-------------------
     if (nvar!=1)
      {
       //Copy global equation numbers into local array
       ivar=new int[nvar];
       
       for(int i=0;i<nvar;i++) 
        {
         //Need to add one to equation numbers
         ivar[i] = 1 + assembly_handler_pt->eqn_number(elem_pt,i);
        }
      }

     // Just one equation
     //------------------
     else
      {

       //Copy global equation numbers into local array - minimum length: 2
       ivar=new int[2];

       // Here's the number of the only equation
       int only_eqn = assembly_handler_pt->eqn_number(elem_pt,0);

       //Need to add one to equation numbers
       ivar[0] = 1 + only_eqn;

       // Now add a dummy eqn either before or after the current one
       int dummy_eqn;
       if (only_eqn!=(n_dof-1))
        {
         dummy_eqn=only_eqn+1;
        }
       else
        {
         dummy_eqn=only_eqn-1;
        }
       ivar[1]=1+dummy_eqn;

       nvar=2;
      }
     
     
     //Now call the frontal routine
     //GB: added check to exclude case with no dofs
     if (nvar)
      {
       MA42AD(nvar,ivar,ndf,last,n_dof,Icntl,Isave,Info);
       //ALH : throw an exception if the frontal_solver fails
       if(Info[0] < 0) {throw NewtonSolverError(true);}
      }
     
     // Cleanup
     delete[] ivar;
     ivar=0;
     
    }
   

   //Loop over the elements again for symbolic factorisation
   //=======================================================
   for(unsigned long e=0;e<n_el;e++)
    {
     //Pointer to the element
     GeneralisedElement* elem_pt = problem_pt->mesh_pt()->element_pt(e);

     // MH: changed array to allocateable
     int* ivar;
     
     //Get number of variables in element
     int nvar = assembly_handler_pt->ndof(elem_pt);
     
     // Multiple equations
     //-------------------
     if (nvar!=1)
      {
       //Copy global equation numbers into local array
       ivar=new int[nvar];
       
       for(int i=0;i<nvar;i++) 
        {
         //Need to add one to equation numbers
         ivar[i] = 1 + assembly_handler_pt->eqn_number(elem_pt,i);
        }
      }
     
     // Just one equation
     //------------------
     else
      {
       
       //Copy global equation numbers into local array - minimum length: 2
       ivar=new int[2];

       // Here's the number of the only equation
       int only_eqn = assembly_handler_pt->eqn_number(elem_pt,0);

       //Need to add one to equation numbers
       ivar[0] = 1 + only_eqn;

       // Now add a dummy eqn either before or after the current one
       int dummy_eqn;
       if (only_eqn!=(n_dof-1))
        {
         dummy_eqn=only_eqn+1;
        }
       else
        {
         dummy_eqn=only_eqn-1;
        }
       ivar[1]=1+dummy_eqn;

       nvar=2;

      }

     //GB: added check to exclude case with no dofs
     if (nvar)
      {
       MA42JD(nvar,ivar,ndf,last,nmaxe,ifsize,Icntl,Isave,Info);
       
       //ALH : throw an exception if the frontal_solver fails
       if(Info[0] < 0) {throw NewtonSolverError(true);}
      }
     
     // Cleanup
     delete[] ivar;
     ivar=0;
     
    } // end of symbolic factorisation
   
   

   // Front sizes: "Somewhat larger than..."
   //---------------------------------------
   front0_recommendation=ifsize[0];
   front1_recommendation=ifsize[1];
   if (!exact_nfront_available)
    {
     nfront[0]=int(ceil(Front_factor*double(front0_recommendation)));
     nfront[1]=int(ceil(Front_factor*double(front1_recommendation)));
     
     if(n_dof < nfront[0]) nfront[0] = n_dof;   
     if(n_dof < nfront[1]) nfront[1] = n_dof;
    }
   //We have a problem if the front size is cocked
   

   // Sizes for lenbuf[]: "Somewhat larger than..."
   //----------------------------------------------
   lenbuf0_recommendation=ifsize[2] + ndf;
   //If we are storing the matrix, 
   //we need to allocate storage for lenbuf_1
   if(Enable_resolve) {lenbuf1_recommendation=ifsize[3];}
   //Otherwise set it to zero
   else {lenbuf1_recommendation=0;}
   lenbuf2_recommendation=ifsize[4];


   //Enable use of direct access files?
   //----------------------------------
   if (Use_direct_access_files)
    {
     if (Doc_stats)
      {
       oomph_info << "Using direct access files" << std::endl;
      }
     
     // Unit numbers for direct access files (middle one only used for
     // re-solve; set to zero as this is not used here).
     int istrm[3];
     istrm[0]=17;
     //If we are storing the matrix, set the stream
     if(Enable_resolve) {istrm[1]=19;}
     //otherwise, set it to zero
     else {istrm[1] = 0;}
     istrm[2]=18;
     
     //Set up the memory: Need memory "somewhat larger" than ...
     double factor=1.1;
     lenbuf[0] = int(ceil(factor*double(10*(nfront[0] + 1))));
     //If we are storing the matrix, set the value
     if(Enable_resolve)
      {lenbuf[1] = int(ceil(factor*double(10*(nfront[0]))));}
     //Otherwise, set to zero
     else {lenbuf[1] = 0;}
     lenbuf[2] = int(ceil(factor*double(10*(2*nfront[0]+5))));
     

     //Initial size allocation: Need memory "somewhat larger" than ...
     if (exact_lenfle_available)
      {
       lenfle[0] = exact_lenfle[0];
       lenfle[1] = exact_lenfle[1];
       lenfle[2] = exact_lenfle[2];
      }
     else
      {
       lenfle0_recommendation=ifsize[2]+ndf;
       lenfle1_recommendation=ifsize[3];
       lenfle2_recommendation=ifsize[4];
       lenfle[0]=int(ceil(Lenfle_factor*double(lenfle0_recommendation)));
       lenfle[1]=int(ceil(Lenfle_factor*double(lenfle1_recommendation)));
       lenfle[2]=int(ceil(Lenfle_factor*double(lenfle2_recommendation)));
      }
     
     //Setup direct access files
     MA42PD(istrm,lenbuf,lenfle,Icntl,Isave,Info);
    }
   else
    {
     if (Doc_stats)
      {
       oomph_info << "Not using direct access files" << std::endl;
      }
     
     //Initial size allocation: Need memory "somewhat larger" than ...
     if (exact_lenbuf_available)
      {
       lenbuf[0] = exact_lenbuf[0];
       lenbuf[1] = exact_lenbuf[1];
       lenbuf[2] = exact_lenbuf[2];
      }
     else
      {
       lenbuf[0] = int(ceil(Lenbuf_factor0*double(lenbuf0_recommendation)));
       //If we are storing the matrix, set the buffer size
       if(Enable_resolve)
        {lenbuf[1] = int(ceil(Lenbuf_factor1*double(lenbuf1_recommendation)));}
       //Otherwise its zero
       else {lenbuf[1] = 0;}
       lenbuf[2] = int(ceil(Lenbuf_factor2*double(lenbuf2_recommendation)));
      }     
    }



   if (Doc_stats)
    {
     oomph_info << "\n FRONT SIZE (min and actual): " 
               <<  ifsize[0] << " " << nfront[0] << std::endl;
    }
     
   
   // Initialise success
   bool success=true;
   
   // Workspace arrays: calculate amount in local variables initiall
   int lw = 1 + lenbuf[0] + lenbuf[1] + nfront[0]*nfront[1];
   if(lrhs*nfront[0] > nrhs*nfront[1])
    {
     lw += lrhs*nfront[0];
    }
   else
    {
     lw += nrhs*nfront[1];
    }
   
   int liw = lenbuf[2] + 2*nfront[0] + 4*nfront[1];

   //Set up memory for w and iw
   //Again allocate dynamically from the heap, rather than the stack   
   //If the required amount of storage has changed, reallocate
   //and store the values in the object member data
   if(lw != Lw)
    {
     //Set the new value of Lw (member data)
     Lw = lw;
     //Delete the previously allocated storage
     if(W) delete[] W;
     //Now allocate new storage
     W = new (std::nothrow) double [Lw];
     if (!W)
      {
       throw OomphLibError("Out of memory in allocation of w\n",
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
    }

   //If the required amount of storage has changed, reallocate
   if(liw != Liw)
    {
     //Set the new value of Liw (member data)
     Liw = liw;
     //Delete the previously allocated storgae
     if(IW!=0) {delete[] IW;}
     //Now allocate new storage
     IW = new (std::nothrow) int [Liw];
     if (!IW)
      {
       throw OomphLibError("Out of memory in allocation of iw\n",
                           OOMPH_CURRENT_FUNCTION,
                           OOMPH_EXCEPTION_LOCATION);
      }
    }

   //Doc
   if (Doc_stats)
    {
     double temp = (lenbuf[0]*sizeof(double) + lenbuf[2]*sizeof(int))/
      (1024.0*1024.0);
     oomph_info << "\n ESTIMATED MEMORY USAGE " << temp << "Mb" << std::endl;
    }
   

   //Now do the actual frontal assembly and solve loop
   //=================================================
   for(unsigned long e=0;e<n_el;e++)
    {
     //Get pointer to the element
     GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);
     
     //Get number of variables in element
     int nvar = assembly_handler_pt->ndof(elem_pt);
     
     // MH: changed array to allocatable
     int* ivar;
     
     // Multiple equations
     //-------------------
     if (nvar!=1)
      {
       //Copy global equation numbers into local array
       ivar=new int[nvar];
       
       for(int i=0;i<nvar;i++) 
        {
         //Need to add one to equation numbers
         ivar[i] = 1 + assembly_handler_pt->eqn_number(elem_pt,i);
        }
       nmaxe = nvar;
      }
     
     // Just one equation
     //------------------
     else
      {
       
       //Copy global equation numbers into local array - minimum length: 2
       ivar=new int[2];
       
       // Here's the number of the only equation
       int only_eqn = assembly_handler_pt->eqn_number(elem_pt,0);
       
       //Need to add one to equation numbers
       ivar[0] = 1 + only_eqn;
       
       // Now add a dummy eqn either before or after the current one
       int dummy_eqn;
       if (only_eqn!=(n_dof-1))
        {
         dummy_eqn=only_eqn+1;
        }
       else
        {
         dummy_eqn=only_eqn-1;
        }
       ivar[1]=1+dummy_eqn;
       
       nmaxe = 2;
      }
     
     //Set up matrices and Vector for call to get_residuals
     Vector<double> residuals(nvar);
     DenseMatrix<double> jacobian(nvar);
     
     //Get the residuals and jacobian
     assembly_handler_pt->get_jacobian(elem_pt,residuals,jacobian);
     
     // Add dummy rows and columns if required
     if (nvar==1)
      {
       double onlyjac=jacobian(0,0);
       jacobian.resize(2,2);
       jacobian(0,0)=onlyjac;
       jacobian(1,0)=0.0;
       jacobian(0,1)=0.0;
       jacobian(1,1)=0.0;
       
       double onlyres=residuals[0];
       residuals.resize(2);
       residuals[0]=onlyres;
       residuals[1]=0.0;
       
       nvar=2;
      }
     
     
     //GB: added check to exclude case with no dofs
     if (nvar)
      {
       //Now translate the residuals and jacobian into form to be
       //taken by stupid FORTRAN frontal solver
       //double avar[nvar][nmaxe], rhs[1][nmaxe];
       //Assign memory on the heap rather than the stack because the ndofs
       //could get too large
       double **avar = new double*[nvar];
       double *tmp = new double[nvar*nmaxe];
       for(int i=0;i<nvar;i++) {avar[i] = &tmp[i*nmaxe];}
       double **rhs = new double*[1];
       rhs[0] = new double[nmaxe];


       for(int i=0;i<nmaxe;i++)
        {
         rhs[0][i] = residuals[i];
         for(int j=0;j<nvar;j++)
          {
           //Note need to transpose here
           avar[j][i] = jacobian(i,j);
          }
        }
       
       //Call the frontal solver
       MA42BD(nvar,ivar,ndf,last,nmaxe,avar,nrhs,rhs,lrhs,n_dof,dx,nfront,
              lenbuf,Lw,W,Liw,IW,Icntl,cntl,Isave,Info,rinfo);
       

       // Error check:
       if (Info[0] < 0)
        {

         if (Doc_stats) oomph_info << "Error: Info[0]=" << Info[0] << std::endl;

         // Solve isn't going to be successful
         success=false;

         // Error: Entries in lenbuf too small -- this is covered
         if (Info[0]==-16)
          {
           if (Doc_stats) oomph_info << "lenbuf[] too small -- can recover..." 
                                    << std::endl;
          }
         else if (Info[0]==-12)
          {
           if (Doc_stats) oomph_info << "nfront[] too small -- can recover..." 
                                    << std::endl;
          }
         else if (Info[0]==-17)
          {
           if (Doc_stats) oomph_info << "lenfle[] too small -- can recover..." 
                                    << std::endl;
          }
         else
          {
           if (Doc_stats)
            {
             oomph_info<<"Can't recover from this error"<<std::endl;
            }
           throw NewtonSolverError(true);
          }
        }

       //Clean up the memory
       delete[] rhs[0]; rhs[0] = 0;
       delete[] rhs; rhs =0;
       delete[] avar[0]; avar[0] = 0;
       delete[] avar;
      }
     
     // Cleanup
     delete[] ivar;
     ivar=0;
     
    } // end of loop over elements for assemble/solve
   
   //Set the sign of the jacobian matrix
   problem_pt->sign_of_jacobian() = Info[1];
   
   //If we are not storing the matrix, then delete the assigned workspace
   if(!Enable_resolve)
    {
     //Free the workspace memory assigned and set stored values to zero
     delete[] IW; IW=0; Liw=0;
     delete[] W;  W=0; Lw=0;
     //Reset the number of dofs
     N_dof=0;
    }
   
   // We've done the elements -- did we encounter an error
   // that forces us to re-allocate workspace?
   if (success)
    {
     // We're done!
     not_done=false;
    }
   else
    {
     //Reset sizes for lenbuf 
     if (!Use_direct_access_files)
      {
       exact_lenbuf[0] = Info[4];
       exact_lenbuf[1] = Info[5];
       exact_lenbuf[2] = Info[6];
       exact_lenbuf_available=true;
      }

     // nfront[] has already been overwritten with required values -- don't
     // need to update anything. Just indicate that new values are 
     // available so they don't have to be re-assigned.
     exact_nfront_available=true;

    //Reset sizes for lenfle
     exact_lenfle[0] = lenbuf[0]*Info[9];
     exact_lenfle[1] = lenbuf[1]*Info[10];
     exact_lenfle[2] = lenbuf[2]*Info[11];
     exact_lenfle_available=true;
     
    }

  } // end of loop over shuffling of workspace array sizes until it all
    // fits...

 result.build(&dist);
 //Now load the values back into result
 for(int i=0;i<n_dof;i++) result[i] = dx[0][i];
 
 //Print the actual memory used
 if (Doc_stats)
  {
   if (!Use_direct_access_files)
    {
     double temp = (Info[4]*sizeof(double) + Info[6]*sizeof(int))/
      (1024.0*1024.0);
     oomph_info << " TOTAL MEMORY USED " << temp << "Mb" << std::endl;
    }
   oomph_info << std::endl;
   oomph_info << "lenbuf[]: " << lenbuf[0] << " " 
             << lenbuf[1] << " " << lenbuf[2] << " " << std::endl;
   oomph_info << "lenbuf[] factors required and initially specified:" 
             << std::endl;
   oomph_info << "lenbuf[0]: " << Info[4]/double(lenbuf0_recommendation) << " " 
             << Lenbuf_factor0 << std::endl;
   oomph_info << "lenbuf[1]: " << Info[5]/double(lenbuf1_recommendation) 
              << " " 
             << Lenbuf_factor1 << std::endl;

   oomph_info << "lenbuf[2]: " << Info[6]/double(lenbuf2_recommendation) << " " 
             << Lenbuf_factor2 << std::endl;
   oomph_info << "nfront[] factors required and initially specified:" 
             << std::endl;
   oomph_info << "nfront[0]: " << nfront[0]/double(front0_recommendation)
             << " " << Front_factor << std::endl;
   oomph_info << "nfront[1]: " << nfront[1]/double(front1_recommendation) 
             << " " << Front_factor << std::endl;
   if (Use_direct_access_files)
    {
     oomph_info << "lenfle[]: " << lenfle[0] << " " 
               << lenfle[1] << " " << lenfle[2] << " " << std::endl;
     oomph_info << "lenfle[] factors required and initially specified:" 
               << std::endl;
     oomph_info << "lenfle[0]: " 
               << Info[9]*lenbuf[0]/double(lenfle0_recommendation) 
               << " " << Lenfle_factor << std::endl;
     oomph_info << "lenfle[1]: " 
               << Info[10]*lenbuf[1]/double(lenfle1_recommendation)
               << " " << Lenfle_factor << std::endl;
     oomph_info << "lenfle[2]: " 
               << Info[11]*lenbuf[2]/double(lenfle2_recommendation)
               << " " << Lenfle_factor << std::endl;
    }
  }
 
 //Free the memory assigned
 delete[] *dx;
 delete dx;
 delete[] last;
}

//====================================================================
/// Wrapper for HSL MA42 frontal solver
//====================================================================
void HSL_MA42::resolve(const DoubleVector &rhs, DoubleVector &result)
{
 //If we haven't stored the factors complain
 if(W==0) 
  {
   throw OomphLibError("Resolve called without first solving",
                       OOMPH_CURRENT_FUNCTION,
                       OOMPH_EXCEPTION_LOCATION);
  }

 //Find the number of dofs (variables)
 int n_dof = N_dof;

 //Check that the number of DOFs is equal to the size of the RHS vector
#ifdef PARANOID
 if(n_dof != static_cast<int>(rhs.nrow()))
  {
   throw OomphLibError(
    "RHS does not have the same dimension as the linear system",
    OOMPH_CURRENT_FUNCTION,
    OOMPH_EXCEPTION_LOCATION);
  }
#endif

 //Special case: just one variable! MA42 can't handle it
 //Solve using the stored jacobian (single double)
 if(n_dof==1) {result[0]= rhs[0]/W[0]; return;}

 //Otherwise actually call the MA42 routine
 //We are solving the original matrix (not the transpose)
 int trans=0;
 //There is only one rhs
 int nrhs=1;
 //The size of the vectors is the number of degrees of freedom in the problem
 int lx = n_dof;

 //Allocate storage for the RHS
 double **b = new double*;
 *b = new double[n_dof];
 //Load in the RHS vector
 for(int i=0;i<n_dof;i++) {b[0][i] = rhs[i];}

 //Storage for the results
 double **x = new double*;
 *x = new double[n_dof];

 //Now call the resolver
 MA42CD(trans,nrhs,lx,b,x,Lw,W,Liw,IW,Icntl,Isave,Info);

 //If there has been an error throw it
 if(Info[0] < 0) {throw NewtonSolverError(true);}

 result.build(this->distribution_pt());
 //Put the answer into the result array
 for(int i=0;i<n_dof;++i) {result[i] = x[0][i];}

 //Delete the allocated storage
 delete[] *x;
 delete x;
 delete[] *b;
 delete b;
}


//=========================================================================
///Function to reorder the elements according to Sloan's algorithm
//=========================================================================
void HSL_MA42::reorder_elements(Problem* const &problem_pt)
{
 //Find the number of elements
 unsigned long n_el = problem_pt->mesh_pt()->nelement();

 //Find the number of dofs
 int n_dof = problem_pt->ndof();

 Vector<int> order(n_el);

 //Control parameters
 int icntl[10], info[15], direct, nsup, vars[n_dof], svar[n_dof];
 int liw, lw, *perm=0;
 double wt[3], rinfo[6];


 // Direct ordering: 1
 direct=1;

 // Initial guess for workspaces (deliberately too small so the
 // routine fails and returns the required sizes
 liw=1;
 lw=1;

 //Initialise things here
 MC63ID(icntl);


 // Suppress printing of error and warning messages?
 if (!Doc_stats)
  {
   icntl[0]=-1;
   icntl[1]=-1;
  }

 //Set up iw and w
 int *iw = new int [liw];
 double *w = new double [lw];

 //Set up the vectors that hold the element connectivity information
 //Can initialise the first entry of eltptr to 1
 Vector<int> eltvar, eltptr(1,1);

 //Locally cache pointer to assembly handler
 AssemblyHandler* const assembly_handler_pt = 
  problem_pt->assembly_handler_pt();
 
 //Now loop over all the elements and assemble eltvar and eltptr
 for(unsigned long e=0;e<n_el;e++)
  {
   //Set up the default order
   order[e] = e;

   //Get pointer to the element
   GeneralisedElement *elem_pt = problem_pt->mesh_pt()->element_pt(e);

   //Get the number of variables in the element via the assembly handler
   int nvar = assembly_handler_pt->ndof(elem_pt);

   // Multiple equations
   //-------------------
   if (nvar!=1)
    {
     //Copy the equation numbers into the global array
     for(int i=0;i<nvar;i++) 
      {
       //Need to add one to equation numbers
       eltvar.push_back(1 + assembly_handler_pt->eqn_number(elem_pt,i));
      }
     eltptr.push_back(eltptr.back()+nvar);
    }
   // Just one equation
   //------------------
   else
    {
     // Here's the number of the only equation
     int only_eqn=assembly_handler_pt->eqn_number(elem_pt,0);

     //Need to add one to equation numbers
     eltvar.push_back(1 + only_eqn);

     // Now add a dummy eqn either before or after the current one
     int dummy_eqn;
     if (only_eqn!=(n_dof-1)) {dummy_eqn=only_eqn+1;}
     else {dummy_eqn=only_eqn-1;}
     
     eltvar.push_back(1+dummy_eqn);
     
     eltptr.push_back(eltptr.back()+2);
    }

  } //End of loop over the elements

 //Set the value of n_e (number of entries in the element variable list
 int n_e = eltvar.size();
 
 do
   {
    //Call the reordering routine
    MC63AD(direct,n_dof,n_el,n_e,&eltvar[0],&eltptr[0],&order[0],
           perm,nsup,vars,svar,wt,liw,iw,lw,w,icntl,info,rinfo);
 
    //If not enough space in iw or w, allocate enough and try again
    if(info[0] == -4)
      {
       if (Doc_stats) oomph_info << " Reallocating liw to " 
                                  << info[4] << std::endl;
       delete[] iw;
       liw = info[4];
       iw = new int[liw];
      }
    
    //If not enough space in w, allocate enough and try again
    if(info[0] == -2)
     {
      if (Doc_stats) oomph_info << " Reallocating lw to " 
                                 << info[5] << std::endl;
      delete[] w;
      lw = info[5];
      w = new double [lw];
     }
    
   } while(info[0] < 0);


 if (Doc_stats)
  {
   oomph_info <<"\n Initial wavefront details  :\n      max "
             << rinfo[0] << "\tmean "
             << rinfo[1] << "\tprofile " << rinfo[2];
   oomph_info <<"\n Reordered wavefront details:\n      max "
             << rinfo[3] << "\tmean "
             << rinfo[4] << "\tprofile " << rinfo[5];
   oomph_info << std::endl;
  }

 //Free the memory allocated
 delete[] w;
 delete[] iw;



 // Store re-ordered elements in temporary vector
 Vector<GeneralisedElement*> tmp_element_pt(n_el);
 for (unsigned e=0;e<n_el;e++)
  {
   tmp_element_pt[e]=problem_pt->mesh_pt()->element_pt(order[e]-1);
  }
 for (unsigned e=0;e<n_el;e++)
  {
   problem_pt->mesh_pt()->element_pt(e)=tmp_element_pt[e];
  }

}

}