LowLevelOperators.cpp

# include <cstdlib>
# include <iostream>
# include <fstream>
# include <iomanip>
# include <cmath>
#include "LowLevelOperators.hpp"
#ifdef USE_TBB_PARALLEL
#include <tbb/task_scheduler_init.h>
#include <tbb/blocked_range.h>
#include <tbb/parallel_reduce.h>
#include <tbb/tbb.h>
#include <tbb/tick_count.h>
#endif


//****************************************************************************80

void atx_cr ( const long int & n, const long int & nz_num, const long int * ia, const long int * ja, const double * a, const double *x,  double * w  )

//****************************************************************************80
//
//  Purpose:
//
//    ATX_CR computes A'*x for a matrix stored in sparse compressed row form.
//
//  Discussion:
//
//    The Sparse Compressed Row storage format is used.
//
//    The matrix A is assumed to be sparse.  To save on storage, only
//    the nonzero entries of A are stored.  The vector JA stores the
//    column index of the nonzero value.  The nonzero values are sorted
//    by row, and the compressed row vector IA then has the property that
//    the entries in A and JA that correspond to row I occur in indices
//    IA[I] through IA[I+1]-1.
//
//    For this version of MGMRES, the row and column indices are assumed
//    to use the C/C++ convention, in which indexing begins at 0.
//
//    If your index vectors IA and JA are set up so that indexing is based 
//    at 1, then each use of those vectors should be shifted down by 1.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    18 July 2007
//
//  Author:
//
//    Original C version by Lili Ju.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Richard Barrett, Michael Berry, Tony Chan, James Demmel,
//    June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
//    Charles Romine, Henk van der Vorst,
//    Templates for the Solution of Linear Systems:
//    Building Blocks for Iterative Methods,
//    SIAM, 1994,
//    ISBN: 0898714710,
//    LC: QA297.8.T45.
//
//    Tim Kelley,
//    Iterative Methods for Linear and Nonlinear Equations,
//    SIAM, 2004,
//    ISBN: 0898713528,
//    LC: QA297.8.K45.
//
//    Yousef Saad,
//    Iterative Methods for Sparse Linear Systems,
//    Second Edition,
//    SIAM, 2003,
//    ISBN: 0898715342,
//    LC: QA188.S17.
//
//  Parameters:
//
//    Input, long int N, the order of the system.
//
//    Input, long int NZ_NUM, the number of nonzeros.
//
//    Input, long int IA[N+1], JA[NZ_NUM], the row and column indices
//    of the matrix values.  The row vector has been compressed.
//
//    Input, double A[NZ_NUM], the matrix values.
//
//    Input, double X[N], the vector to be multiplied by A'.
//
//    Output, double W[N], the value of A'*X.
//
{

  for (long int i = 0; i < n; i++ )
  {
    w[i] = 0.0;
    long int k1 = ia[i];
    long int k2 = ia[i+1];
    for (long int  k = k1; k < k2; k++ )
    {
      w[ja[k]] = w[ja[k]] + a[k] * x[i];
    }
  }
}


//****************************************************************************80

void ax_cr ( const long int & n, const long int & nz_num, const long int * ia, const long int * ja, const double * a, const double *x,  double * w )

//****************************************************************************80
//
//  Purpose:
//
//    AX_CR computes A*x for a matrix stored in sparse compressed row form.
//
//  Discussion:
//
//    The Sparse Compressed Row storage format is used.
//
//    The matrix A is assumed to be sparse.  To save on storage, only
//    the nonzero entries of A are stored.  The vector JA stores the
//    column index of the nonzero value.  The nonzero values are sorted
//    by row, and the compressed row vector IA then has the property that
//    the entries in A and JA that correspond to row I occur in indices
//    IA[I] through IA[I+1]-1.
//
//    For this version of MGMRES, the row and column indices are assumed
//    to use the C/C++ convention, in which indexing begins at 0.
//
//    If your index vectors IA and JA are set up so that indexing is based 
//    at 1, then each use of those vectors should be shifted down by 1.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    18 July 2007
//
//  Author:
//
//    Original C version by Lili Ju.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Richard Barrett, Michael Berry, Tony Chan, James Demmel,
//    June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
//    Charles Romine, Henk van der Vorst,
//    Templates for the Solution of Linear Systems:
//    Building Blocks for Iterative Methods,
//    SIAM, 1994,
//    ISBN: 0898714710,
//    LC: QA297.8.T45.
//
//    Tim Kelley,
//    Iterative Methods for Linear and Nonlinear Equations,
//    SIAM, 2004,
//    ISBN: 0898713528,
//    LC: QA297.8.K45.
//
//    Yousef Saad,
//    Iterative Methods for Sparse Linear Systems,
//    Second Edition,
//    SIAM, 2003,
//    ISBN: 0898715342,
//    LC: QA188.S17.
//
//  Parameters:
//
//    Input, long int N, the order of the system.
//
//    Input, long int NZ_NUM, the number of nonzeros.
//
//    Input, long int IA[N+1], JA[NZ_NUM], the row and column indices
//    of the matrix values.  The row vector has been compressed.
//
//    Input, double A[NZ_NUM], the matrix values.
//
//    Input, double X[N], the vector to be multiplied by A.
//
//    Output, double W[N], the value of A*X.
//
{
#ifdef USE_TBB_PARALLEL
  ax_cr_parallel axc_par;
  axc_par.aij=a;
  axc_par.I=ia;
  axc_par.J=ja;
  axc_par.xA=x;
  axc_par.resArray=w;
  size_t length=static_cast<size_t>(n);
  tbb::parallel_for(tbb::blocked_range<size_t>(0,static_cast<size_t>(length)), axc_par);
#else
  for (long int i = 0; i < n; i++ )
  {
    w[i] = 0.0;
    long int k1 = ia[i];
    long int k2 = ia[i+1];
    for (long int k = k1; k < k2; k++ )
    {
      w[i] = w[i] + a[k] * x[ja[k]];
    }
  }
#endif  
}


//****************************************************************************80

void diagonal_pointer_cr ( const long int & n, const long int & nz_num, const long int * ia, const long int * ja, long int * ua )

//****************************************************************************80
//
//  Purpose:
//
//    DIAGONAL_POINTER_CR finds diagonal entries in a sparse compressed row matrix.
//
//  Discussion:
//
//    The matrix A is assumed to be stored in compressed row format.  Only
//    the nonzero entries of A are stored.  The vector JA stores the
//    column index of the nonzero value.  The nonzero values are sorted
//    by row, and the compressed row vector IA then has the property that
//    the entries in A and JA that correspond to row I occur in indices
//    IA[I] through IA[I+1]-1.
//
//    The array UA can be used to locate the diagonal elements of the matrix.
//
//    It is assumed that every row of the matrix includes a diagonal element,
//    and that the elements of each row have been ascending sorted.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    18 July 2007
//
//  Author:
//
//    Original C version by Lili Ju.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, long int N, the order of the system.
//
//    Input, long int NZ_NUM, the number of nonzeros.
//
//    Input, long int IA[N+1], JA[NZ_NUM], the row and column indices
//    of the matrix values.  The row vector has been compressed.  On output,
//    the order of the entries of JA may have changed because of the sorting.
//
//    Output, long int UA[N], the index of the diagonal element of each row.
//
{
#ifdef USE_TBB_PARALLEL
  size_t length=static_cast<size_t>(n);
  diagonal_pointer_cr_parallel dppar;
  dppar.inputArrayRow=ia;
  dppar.inputArrayCol=ja;
  dppar.outputArray=ua;
  tbb::parallel_for(tbb::blocked_range<size_t>(0,length), dppar);
#else
  for (long int i = 0; i < n; i++ )
  {
    ua[i] = -1;
    long int j1 = ia[i];
    long int j2 = ia[i+1];

    for (long int j = j1; j < j2; j++ )
    {
      if( ja[j] == i ) 
      {
        ua[i] = j;
        break;
      }
    }
  }
#endif  
}


//****************************************************************************80

void ilu_cr ( const long int & n, const long int & nz_num, const long int * ia, const long int * ja, const double * a, long int * ua,  double * l )

//****************************************************************************80
//
//  Purpose:
//
//    ILU_CR computes the incomplete LU factorization of a matrix.
//
//  Discussion:
//
//    The matrix A is assumed to be stored in compressed row format.  Only
//    the nonzero entries of A are stored.  The vector JA stores the
//    column index of the nonzero value.  The nonzero values are sorted
//    by row, and the compressed row vector IA then has the property that
//    the entries in A and JA that correspond to row I occur in indices
//    IA[I] through IA[I+1]-1.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    25 July 2007
//
//  Author:
//
//    Original C version by Lili Ju.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, long int N, the order of the system.
//
//    Input, long int NZ_NUM, the number of nonzeros.
//
//    Input, long int IA[N+1], JA[NZ_NUM], the row and column indices
//    of the matrix values.  The row vector has been compressed.
//
//    Input, double A[NZ_NUM], the matrix values.
//
//    Input, long int UA[N], the index of the diagonal element of each row.
//
//    Output, double L[NZ_NUM], the ILU factorization of A.
//
{
    long int *iw= new long int[n];
    //
    //  Copy A.
    //
    vectorCopy<double>(a, l,  nz_num);

    for (long int i = 0; i < n; i++ ) 
    {
        //
        //  IW points to the nonzero entries in row I.
        //
        initializeVector<long int>(iw, -1, n);
        for (long int k = ia[i]; k <= ia[i+1] - 1; k++ ) 
        {
            iw[ja[k]] = k;
        }

        long int j = ia[i];
        long int jrow=0;
        do 
        {
            jrow = ja[j];
            if ( i <= jrow )
            {
                break;
            }
            double tl = l[j] * l[ua[jrow]];
            l[j] = tl;
            for (long int  jj = ua[jrow] + 1; jj <= ia[jrow+1] - 1; jj++ ) 
            {
                long int jw = iw[ja[jj]];
                if ( jw != -1 ) 
                {
                    l[jw] = l[jw] - tl * l[jj];
                }
            }
            j = j + 1;
        } while ( j <= ia[i+1] - 1 );
        ua[i] = j;
        if ( jrow != i ) 
        {
            std::cout << "\n";
            std::cout << "ILU_CR - Fatal error!\n";
            std::cout << "  JROW != I\n";
            std::cout << "  JROW = " << jrow << "\n";
            std::cout << "  I    = " << i << "\n";
            exit ( 1 );
        }

        if ( l[j] == 0.0 ) 
        {
            std::cout << "\n";
            std::cout << "ILU_CR - Fatal error!\n";
            std::cout << "  Zero pivot on step I = " << i << "\n";
            std::cout << "  L[" << j << "] = 0.0\n";
            exit ( 1 );
        }
        l[j] = 1.0 / l[j];
    }

    for (long int k = 0; k < n; k++ ) 
    {
        l[ua[k]] = 1.0 / l[ua[k]];
    }
    delete [] iw;
    iw=NULL;
}


//****************************************************************************80

void lus_cr ( const long int & n, const long int & nz_num, const long int *ia, const long int *ja, const double *l, const long int *ua,   const double *r, double * z )

//****************************************************************************80
//
//  Purpose:
//
//    LUS_CR applies the incomplete LU preconditioner.
//
//  Discussion:
//
//    The linear system M * Z = R is solved for Z.  M is the incomplete
//    LU preconditioner matrix, and R is a vector supplied by the user.
//    So essentially, we're solving L * U * Z = R.
//
//    The matrix A is assumed to be stored in compressed row format.  Only
//    the nonzero entries of A are stored.  The vector JA stores the
//    column index of the nonzero value.  The nonzero values are sorted
//    by row, and the compressed row vector IA then has the property that
//    the entries in A and JA that correspond to row I occur in indices
//    IA[I] through IA[I+1]-1.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    18 July 2007
//
//  Author:
//
//    Original C version by Lili Ju.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, long int N, the order of the system.
//
//    Input, long int NZ_NUM, the number of nonzeros.
//
//    Input, long int IA[N+1], JA[NZ_NUM], the row and column indices
//    of the matrix values.  The row vector has been compressed.
//
//    Input, double L[NZ_NUM], the matrix values.
//
//    Input, long int UA[N], the index of the diagonal element of each row.
//
//    Input, double R[N], the right hand side.
//
//    Output, double Z[N], the solution of the system M * Z = R.
//
{
  
    double *w= new double[n];
    vectorCopy<double>(r, w, n);

//
//  Solve L * w = w where L is unit lower triangular.
//
  for ( long int i = 1; i < n; i++ )
  {
    for ( long int j = ia[i]; j < ua[i]; j++ )
    {
      w[i] = w[i] - l[j] * w[ja[j]];
    }
  }
//
//  Solve U * w = w, where U is upper triangular.
//
  for (long int i = n - 1; 0 <= i; i-- ) 
  {
    for (long int j = ua[i] + 1; j < ia[i+1]; j++ ) 
    {
      w[i] = w[i] - l[j] * w[ja[j]];
    }
    w[i] = w[i] / l[ua[i]];
  }
//
//  Copy Z out.
//
    vectorCopy<double>(w,z, n);
    delete [] w;
    w=NULL;
}


//****************************************************************************80

void mult_givens ( const double & c, const double & s, const long int & k, double * g )

//****************************************************************************80
//
//  Purpose:
//
//    MULT_GIVENS applies a Givens rotation to two successive entries of a vector.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    08 August 2006
//
//  Author:
//
//    Original C version by Lili Ju.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Richard Barrett, Michael Berry, Tony Chan, James Demmel,
//    June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
//    Charles Romine, Henk van der Vorst,
//    Templates for the Solution of Linear Systems:
//    Building Blocks for Iterative Methods,
//    SIAM, 1994,
//    ISBN: 0898714710,
//    LC: QA297.8.T45.
//
//    Tim Kelley,
//    Iterative Methods for Linear and Nonlinear Equations,
//    SIAM, 2004,
//    ISBN: 0898713528,
//    LC: QA297.8.K45.
//
//    Yousef Saad,
//    Iterative Methods for Sparse Linear Systems,
//    Second Edition,
//    SIAM, 2003,
//    ISBN: 0898715342,
//    LC: QA188.S17.
//
//  Parameters:
//
//    Input, double C, S, the cosine and sine of a Givens
//    rotation.
//
//    Input, long int K, indicates the location of the first vector entry.
//
//    Input/output, double G[K+2], the vector to be modified.  On output,
//    the Givens rotation has been applied to entries G(K) and G(K+1).
//
{

  double g1 = c * g[k] - s * g[k+1];
  double g2 = s * g[k] + c * g[k+1];
  g[k]   = g1;
  g[k+1] = g2;
}







void pmgmres_ilu_cr ( const long int & n, const long int & nz_num, 
                      const long int * ia,  long int * ja,  double *a,   double * x, const double *rhs, 
                      const long int & itr_max, const long int & mr, const double & tol_abs, 
                      const double & tol_rel, short int verbose )


//****************************************************************************80
//
//  Purpose:
//
//    PMGMRES_ILU_CR applies the preconditioned restarted GMRES algorithm.
//
//  Discussion:
//
//    The matrix A is assumed to be stored in compressed row format.  Only
//    the nonzero entries of A are stored.  The vector JA stores the
//    column index of the nonzero value.  The nonzero values are sorted
//    by row, and the compressed row vector IA then has the property that
//    the entries in A and JA that correspond to row I occur in indices
//    IA[I] through IA[I+1]-1.
//
//    This routine uses the incomplete LU decomposition for the
//    preconditioning.  This preconditioner requires that the sparse
//    matrix data structure supplies a storage position for each diagonal
//    element of the matrix A, and that each diagonal element of the
//    matrix A is not zero.
//
//    Thanks to Jesus Pueblas Sanchez-Guerra for supplying two
//    corrections to the code on 31 May 2007.
//
//
//    This implementation of the code stores the doubly-dimensioned arrays
//    H and V as vectors.  However, it follows the C convention of storing
//    them by rows, rather than my own preference for storing them by
//    columns.   I may come back and change this some time.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    26 July 2007
//
//  Author:
//
//    Original C version by Lili Ju.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Richard Barrett, Michael Berry, Tony Chan, James Demmel,
//    June Donato, Jack Dongarra, Victor Eijkhout, Roidan Pozo,
//    Charles Romine, Henk van der Vorst,
//    Templates for the Solution of Linear Systems:
//    Building Blocks for Iterative Methods,
//    SIAM, 1994.
//    ISBN: 0898714710,
//    LC: QA297.8.T45.
//
//    Tim Kelley,
//    Iterative Methods for Linear and Nonlinear Equations,
//    SIAM, 2004,
//    ISBN: 0898713528,
//    LC: QA297.8.K45.
//
//    Yousef Saad,
//    Iterative Methods for Sparse Linear Systems,
//    Second Edition,
//    SIAM, 2003,
//    ISBN: 0898715342,
//    LC: QA188.S17.
//
//  Parameters:
//
//    Input, long int N, the order of the linear system.
//
//    Input, long int NZ_NUM, the number of nonzero matrix values.
//
//    Input, long int IA[N+1], JA[NZ_NUM], the row and column indices
//    of the matrix values.  The row vector has been compressed.
//
//    Input, double A[NZ_NUM], the matrix values.
//
//    Input/output, double X[N]; on input, an approximation to
//    the solution.  On output, an improved approximation.
//
//    Input, double RHS[N], the right hand side of the linear system.
//
//    Input, long int ITR_MAX, the maximum number of (outer) iterations to take.
//
//    Input, long int MR, the maximum number of (inner) iterations to take.
//    MR must be less than N.
//
//    Input, double TOL_ABS, an absolute tolerance applied to the
//    current residual.
//
//    Input, double TOL_REL, a relative tolerance comparing the
//    current residual to the initial residual.
//
{
    double delta = 1.0e-03;
    double rho_tol=0.0;  
    double *c = new double[mr+1];
    double *g = new double[mr+1];
    double *h = new double[(mr+1)*mr];
    double *l= new double[ia[n]+1];
    double *r= new double[n];
    double *s= new double[mr+1];
    long int *ua= new long int[n];
    double *v= new double[(mr+1)*n];
    double *y= new double[mr+1];

    if ( verbose>1 )
    {
        std::cout<<"Rearranging CR..."<<std::endl;
    }
    rearrange_cr ( n, nz_num, ia, ja, a );

    if ( verbose>1 )
    {
        std::cout<<"...done, diagonal pointer to cr..."<<std::endl;
    }
    diagonal_pointer_cr ( n, nz_num, ia, ja, ua );

    if ( verbose>1 )
    {
        std::cout<<"...done, ilu_cr..."<<std::endl;
    }
    ilu_cr ( n, nz_num, ia, ja, a, ua, l );

    if ( verbose>0 )
    {
        if ( verbose>1 )
        {
            std::cout<<"...done"<<std::endl;
        }
        std::cout << "\n";
        std::cout << "PMGMRES_ILU_CR\n";
        std::cout << "  Number of unknowns = " << n << "\n";
    }
    double rho = 1.e32;
    long int  itr_used = 0;
    for (long int itr = 0; itr < itr_max; itr++ ) 
    {
        ax_cr ( n, nz_num, ia, ja, a, x, r );
        vxmvVectorUpdating<double>(r, rhs,n);
        lus_cr ( n, nz_num, ia, ja, l, ua, r, r );
        rho = sqrt ( r8vec_dot ( n, r, r ) );
        if ( verbose>1 )
        {
            std::cout << "  ITR = " << itr << "  Residual = " << rho << "\n";
        }

        if ( itr == 0 )
        {
            rho_tol = rho * tol_rel;
        }
        double invrho=1.0/rho;
        yaxsFunc<double>(r, v, invrho,n);
        g[0] = rho;
        
        for (long int i = 1; i < mr + 1; i++ ) 
        {
            g[i] = 0.0;
        }

        for (long int i = 0; i < mr + 1; i++ ) 
        {
            for (long int j = 0; j < mr; j++ ) 
            {
                h[i*(mr)+j] = 0.0;
            }
        }
        long int k_copy=-1;
        for (long int  k = 0; k < mr; k++ )
        {
            k_copy = k;
            ax_cr ( n, nz_num, ia, ja, a, v+k*n, v+(k+1)*n ); 
            lus_cr ( n, nz_num, ia, ja, l, ua, v+(k+1)*n, v+(k+1)*n );
            double av = sqrt ( r8vec_dot ( n, v+(k+1)*n, v+(k+1)*n ) );
            for ( long int j = 0; j <= k; j++ ) 
            {
                h[j*mr+k] = r8vec_dot ( n, v+(k+1)*n, v+j*n );
                for (long int  i = 0; i < n; i++ ) 
                {
                    v[(k+1)*n+i] = v[(k+1)*n+i] - h[j*mr+k] * v[j*n+i];
                }
            }
            h[(k+1)*mr+k] = sqrt ( r8vec_dot ( n, v+(k+1)*n, v+(k+1)*n ) );
            if ( ( av + delta * h[(k+1)*mr+k]) == av ) 
            {
                for (long int j = 0; j < k + 1; j++ )
                {
                    double htmp = r8vec_dot ( n, v+(k+1)*n, v+j*n );
                    h[j*mr+k] = h[j*mr+k] + htmp;
                    for (long int  i = 0; i < n; i++ ) 
                    {
                        v[(k+1)*n+i] = v[(k+1)*n+i] - htmp * v[j*n+i];
                    }
                }
                h[(k+1)*mr+k] = sqrt ( r8vec_dot ( n, v+(k+1)*n, v+(k+1)*n ) );
            }
            if ( h[(k+1)*mr+k] != 0.0 )
            {
                for (long int i = 0; i < n; i++ )
                {
                    v[(k+1)*n+i] = v[(k+1)*n+i] / h[(k+1)*mr+k];
                } 
            }

            if ( 0 < k )  
            {
                for (long int  i = 0; i < k + 2; i++ ) 
                {
                    y[i] = h[i*mr+k];
                }
                for (long  j = 0; j < k; j++ ) 
                {
                    mult_givens ( c[j], s[j], j, y );
                }
                for (long int  i = 0; i < k + 2; i++ )
                {
                    h[i*mr+k] = y[i];
                }
            }
            double mu = sqrt ( h[k*mr+k] * h[k*mr+k] + h[(k+1)*mr+k] * h[(k+1)*mr+k] );
            c[k] = h[k*mr+k] / mu;
            s[k] = -h[(k+1)*mr+k] / mu;
            h[k*mr+k] = c[k] * h[k*mr+k] - s[k] * h[(k+1)*mr+k];
            h[(k+1)*mr+k] = 0.0;
            mult_givens ( c[k], s[k], k, g );
            rho = fabs ( g[k+1] );
            itr_used = itr_used + 1;
            if ( verbose>2 )
            {
                std::cout << "  K   = " << k << "  Residual = " << rho << "\n";
            }

            if ( rho <= rho_tol && rho <= tol_abs )
            {
                break;
            }
        }
        y[k_copy] = g[k_copy] / h[k_copy*mr+k_copy];
        for ( long int i = k_copy - 1; 0 <= i; i-- )
        {
            y[i] = g[i];
            for ( long int j = i + 1; j < k_copy + 1; j++ ) 
            {
                y[i] = y[i] - h[i*mr+j] * y[j];
            }
            y[i] = y[i] / h[i*mr+i];
        }
        xxvyFunc<double>(x, y, v, k_copy, n);
        if ( rho <= rho_tol && rho <= tol_abs )
        {
            break;
        }
    }//end on solver iterations

    if ( verbose>0 )
    {
        std::cout << "\n";;
        std::cout << "PMGMRES_ILU_CR:\n";
        std::cout << "  Iterations = " << itr_used << "\n";
        std::cout << "  Final residual = " << rho << "\n";
    }

    delete [] c;
    delete [] g;
    delete [] h;
    delete [] l;
    delete [] r;
    delete [] s;
    delete [] ua;
    delete [] v;
    delete [] y;

    c=NULL;
    g=NULL;
    h=NULL;
    l=NULL;
    r=NULL;
    s=NULL;
    ua=NULL;
    v=NULL;
    y=NULL;
}







//****************************************************************************80

double r8vec_dot ( const long int & n, const double * a1, const double * a2 )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_DOT computes the dot product of a pair of R8VEC's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    03 July 2005
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, long int N, the number of entries in the vectors.
//
//    Input, double A1[N], A2[N], the two vectors to be considered.
//
//    Output, double R8VEC_DOT, the dot product of the vectors.
//
{
  double value=0.0;

#ifdef USE_TBB_PARALLEL
  size_t length=static_cast<size_t>(n);
  r8vec_dot_parallel r8dot(a1,a2);
  tbb::parallel_reduce(tbb::blocked_range<size_t>(0,length), r8dot );
  value=r8dot.getSum();
#else
  for (long int  i = 0; i < n; i++ )
  {
    value = value + a1[i] * a2[i];
  }
#endif
  return(value);
}

//****************************************************************************80

void rearrange_cr( const long int & n, const long int & nz_num, const long int * ia, long int * ja, double * a )

//****************************************************************************80
//
//  Purpose:
//
//    REARRANGE_CR sorts a sparse compressed row matrix.
//
//  Discussion:
//
//    This routine guarantees that the entries in the CR matrix
//    are properly sorted.
//
//    After the sorting, the entries of the matrix are rearranged in such
//    a way that the entries of each column are listed in ascending order
//    of their column values.
//
//    The matrix A is assumed to be stored in compressed row format.  Only
//    the nonzero entries of A are stored.  The vector JA stores the
//    column index of the nonzero value.  The nonzero values are sorted
//    by row, and the compressed row vector IA then has the property that
//    the entries in A and JA that correspond to row I occur in indices
//    IA[I] through IA[I+1]-1.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    18 July 2007
//
//  Author:
//
//    Original C version by Lili Ju.
//    C++ version by John Burkardt.
//
//  Parameters:
//
//    Input, long int N, the order of the system.
//
//    Input, long int NZ_NUM, the number of nonzeros.
//
//    Input, long int IA[N+1], the compressed row index.
//
//    Input/output, long int JA[NZ_NUM], the column indices.  On output,
//    the order of the entries of JA may have changed because of the sorting.
//
//    Input/output, double A[NZ_NUM], the matrix values.  On output, the
//    order of the entries may have changed because of the sorting.
//
{
  for (long int i = 0; i < n; i++ )
  {
    long int j1 = ia[i];
    long int j2 = ia[i+1];
    rearrange_cr_row(j1,j2,ja,a);
  }
}


void rearrange_cr_row ( const long int & j1,const long int & j2, long int * ja, double * a )
{
    long int is = j2 - j1;
    for (long int  k = 1; k < is; k++ ) 
    {
      for (long int  j = j1; j < j2 - k; j++ ) 
      {
        if ( ja[j+1] < ja[j] ) 
        {
          long int itemp = ja[j+1];
          ja[j+1] =  ja[j];
          ja[j] =  itemp;

          double dtemp = a[j+1];
          a[j+1] =  a[j];
          a[j] = dtemp;
        }
      }
    }
}