# include <cmath>
# include <cstdlib>
# include <cstring>
# include <ctime>
# include <fstream>
# include <iomanip>
# include <iostream>

using namespace std;

# include "hb_io.hpp"

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

bool ch_eqi ( char c1, char c2 )

//****************************************************************************80
//
//  Purpose:
//
//    CH_EQI is true if two characters are equal, disregarding case.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 June 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char C1, C2, the characters to compare.
//
//    Output, bool CH_EQI, is true if the two characters are equal,
//    disregarding case.
//
{
  if ( 97 <= c1 && c1 <= 122 ) 
  {
    c1 = c1 - 32;
  } 
  if ( 97 <= c2 && c2 <= 122 ) 
  {
    c2 = c2 - 32;
  }     

  return ( c1 == c2 );
}
//****************************************************************************80

bool ch_is_digit ( char c )

//****************************************************************************80
//
//  Purpose:
//
//    CH_IS_DIGIT returns TRUE if a character is a decimal digit.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 December 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char C, the character to be analyzed.
//
//    Output, bool CH_IS_DIGIT, is TRUE if C is a digit.
//
{
  if ( '0' <= c && c <= '9' )
  {
    return true;
  }
  else
  {
    return false;
  }
}
//****************************************************************************80

bool ch_is_format_code ( char c )

//****************************************************************************80
//
//  Purpose:
//
//    CH_IS_FORMAT_CODE returns TRUE if a character is a FORTRAN format code.
//
//  Discussion:
//
//    The format codes accepted here are not the only legal format
//    codes in FORTRAN90.  However, they are more than sufficient
//    for my needs!
//
//  Table:
//
//    A  Character
//    B  Binary digits
//    D  Real number, exponential representation
//    E  Real number, exponential representation
//    F  Real number, fixed point
//    G  General format
//    I  Integer
//    L  Logical variable
//    O  Octal digits
//    Z  Hexadecimal digits
//    *  Free format
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 November 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char C, the character to be analyzed.
//
//    Output, bool CH_IS_FORMAT_CODE, is TRUE if C is a FORTRAN format code.
//
{
       if ( ch_eqi ( c, 'A' ) ) 
  {
    return true;
  }
  else if ( ch_eqi ( c, 'B' ) )
  {
    return true;
  }
  else if ( ch_eqi ( c, 'D' ) )
  {
    return true;
  }
  else if ( ch_eqi ( c, 'E' ) )
  {
    return true;
  }
  else if ( ch_eqi ( c, 'F' ) )
  {
    return true;
  }
  else if ( ch_eqi ( c, 'G' ) )
  {
    return true;
  }
  else if ( ch_eqi ( c, 'I' ) )
  {
    return true;
  }
  else if ( ch_eqi ( c, 'L' ) )
  {
    return true;
  }
  else if ( ch_eqi ( c, 'O' ) )
  {
    return true;
  }
  else if ( ch_eqi ( c, 'Z' ) )
  {
    return true;
  }
  else if ( c == '*' )
  {
    return true;
  }
  else
  {
    return false;
  }
}
//****************************************************************************80

int ch_to_digit ( char c )

//****************************************************************************80
//
//  Purpose:
//
//    CH_TO_DIGIT returns the integer value of a base 10 digit.
//
//  Example:
//
//     C   DIGIT
//    ---  -----
//    '0'    0
//    '1'    1
//    ...  ...
//    '9'    9
//    ' '    0
//    'X'   -1
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 June 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char C, the decimal digit, '0' through '9' or blank are legal.
//
//    Output, int CH_TO_DIGIT, the corresponding integer value.  If C was
//    'illegal', then DIGIT is -1.
//
{
  int digit;

  if ( '0' <= c && c <= '9' )
  {
    digit = c - '0';
  }
  else if ( c == ' ' )
  {
    digit = 0;
  }
  else
  {
    digit = -1;
  }

  return digit;
}
//****************************************************************************80

void hb_exact_read ( ifstream &input, int nrow, int nrhs, int rhscrd, 
  char *rhsfmt, char *rhstyp, double exact[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_EXACT_READ reads the exact solution vectors in an HB file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ifstream &INPUT, the unit from which data is read.
//
//    Input, int NROW, the number of rows or variables.
//
//    Input, int NRHS, the number of right hand sides.
//
//    Input, int RHSCRD, the number of lines in the file for
//    right hand sides.
//
//    Input, char *RHSFMT, the 20 character format for reading values
//    of the right hand side.
//
//    Input, char *RHSTYP, the 3 character right hand side type.
//    First character is F for full storage or M for same as matrix.
//    Second character is G if starting "guess" vectors are supplied.
//    Third character is X if exact solution vectors are supplied.
//    Ignored if NRHS = 0.
//
//    Output, double EXACT[NROW*NRHS], the exact solution vectors.
//
{
  char code;
  int i;
  int j;
  int jhi;
  int jlo;
  int khi;
  int klo;
  char line[255];
  int line_num;
  int m;
  int r;
  char *s;
  int w;

  if ( 0 < rhscrd )
  {
    if ( rhstyp[2] == 'X' )
    {
      s_to_format ( rhsfmt, &r, &code, &w, &m );

      line_num = 1 + ( nrow * nrhs - 1 ) / r;

      jhi = 0;
      for ( i = 1; i <= line_num; i++ )
      {
        input.getline ( line, sizeof ( line ) );
        jlo = jhi + 1;
        jhi = i4_min ( jlo + r - 1, nrow * nrhs );

        khi = 0;
        for ( j = jlo; j <= jhi; j++ )
        {
          klo = khi + 1;
          khi = i4_min ( klo + w - 1, strlen ( line ) );
          s = s_substring ( line, klo, khi );
          exact[j-1] = atof ( s );
        }
      }
    }
  }

  return;
}
//****************************************************************************80

void hb_exact_write ( ofstream &output, int nrow, int nrhs, int rhscrd, 
  char *rhsfmt, char *rhstyp, double exact[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_EXACT_WRITE writes the exact solution vectors to an HB file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ofstream &OUTPUT, the unit to which data is written.
//
//    Input, int NROW, the number of rows or variables.
//
//    Input, int NRHS, the number of right hand sides.
//
//    Input, int RHSCRD, the number of lines in the file for
//    right hand sides.
//
//    Input, char *RHSFMT, the 20 character format for reading values
//    of the right hand side.
//
//    Input, char *RHSTYP, the 3 character right hand side type.
//    First character is F for full storage or M for same as matrix.
//    Second character is G if starting "guess" vectors are supplied.
//    Third character is X if exact solution vectors are supplied.
//    Ignored if NRHS = 0.
//
//    Input, double EXACT[NROW*NRHS], the exact solution vectors.
//
{
  char code;
  int i;
  int j;
  int jhi;
  int jlo;
  int line_num;
  int m;
  int r;
  int w;

  if ( 0 < rhscrd )
  {
    if ( rhstyp[2] == 'X' )
    {
      s_to_format ( rhsfmt, &r, &code, &w, &m );
      line_num = 1 + ( nrow * nrhs - 1 ) / r;

      jhi = 0;
      for ( i = 1; i <= line_num; i++ )
      {
        jlo = jhi + 1;
        jhi = i4_min ( jlo + r - 1, nrow * nrhs );

        for ( j = jlo; j <= jhi; j++ )
        {
          output << setw(w) << exact[j-1];
        }
        output << "\n";
      }
    }
  }
  return;
}
//****************************************************************************80

void hb_file_read ( ifstream &input, char **title, char **key, int *totcrd, 
  int *ptrcrd, int *indcrd, int *valcrd, int *rhscrd, char **mxtype, int *nrow, 
  int *ncol, int *nnzero, int *neltvl, char **ptrfmt, char **indfmt, char **valfmt, 
  char **rhsfmt, char **rhstyp, int *nrhs, int *nrhsix, int **colptr, 
  int **rowind, double **values, double **rhsval, int **rhsptr, int **rhsind,  
  double **rhsvec, double **guess, double **exact )

//****************************************************************************80
//
//  Purpose:
//
//    HB_FILE_READ reads an HB file.
//
//  Discussion:
//
//    This routine reads all the information from an HB file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    15 August 2016
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ifstream &INPUT, the unit from which the data is read.
//
//    Output, char *TITLE, a 72 character title for the matrix.
//
//    Output, char *KEY, an 8 character identifier for the matrix.
//
//    Output, int *TOTCRD, the total number of lines of data.
//
//    Output, int *PTRCRD, the number of input lines for pointers.
//
//    Output, int *INDCRD, the number of input lines for row indices.
//
//    Output, int *VALCRD, the number of input lines for numerical values.
//
//    Output, int *RHSCRD, the number of input lines for right hand sides.
//
//    Output, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Output, int *NROW, the number of rows or variables.
//
//    Output, int *NCOL, the number of columns or elements.
//
//    Output, int *NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Output, int *NELTVL, the number of finite element matrix entries,
//    set to 0 in the case of assembled matrices.
//
//    Output, char *PTRFMT, the 16 character format for reading pointers.
//
//    Output, char *INDFMT, the 16 character format for reading indices.
//
//    Output, char *VALFMT, the 20 character format for reading values.
//
//    Output, char *RHSFMT, the 20 character format for reading values
//    of the right hand side.
//
//    Output, char *RHSTYP, the 3 character right hand side type.
//    First character is F for full storage or M for same as matrix.
//    Second character is G if starting "guess" vectors are supplied.
//    Third character is X if exact solution vectors are supplied.
//
//    Output, int *NRHS, the number of right hand sides.
//
//    Output, int *NRHSIX, the number of entries of storage for right
//    hand side values, in the case where RHSTYP[0] = 'M' and
//    MXTYPE[2] = 'A'.
//
//    Output, int COLPTR[NCOL+1], COLPTR[I-1] points to the location of
//    the first entry of column I in the sparse matrix structure.
//
//    If MXTYPE[2] == 'A':
//
//      Output, int ROWIND[NNZERO], the row index of each item.
//
//    If MXTYPE[2] == 'F':
//
//      Output, int ROWIND[NELTVL], the row index of each item.
//
//    If RHSTYP[0] == 'F':
//
//      Output, double RHSVAL[NROW*NRHS], contains NRHS dense right hand
//      side vectors.
//
//      Output, int RHSPTR[], is not used.
//
//      Output, int RHSIND[], is not used.
//
//      Output, int RHSVEC[], is not used.
//
//    If RHSTYP[0] = 'M' and MXTYPE[2] = 'A':
//
//      Output, double RHSVAL[], is not used.
//
//      Output, int RHSPTR[NRHS+1], RHSPTR[I-1] points to the location of
//      the first entry of right hand side I in the sparse right hand
//      side vector.
//
//      Output, int RHSIND[NRHSIX], indicates, for each entry of
//      RHSVEC, the corresponding row index.
//
//      Output, double RHSVEC[NRHSIX], contains the value of the right hand
//      side entries.
//
//    If RHSTYP[0] = 'M' and MXTYPE[2] = 'E':
//
//      Output, double RHSVAL[NNZERO*NRHS], contains NRHS unassembled
//      finite element vector right hand sides.
//
//      Output, int RHSPTR[], is not used.
//
//      Output, int RHSIND[], is not used.
//
//      Output, double RHSVEC[], is not used.
//
//    Output, double GUESS[NROW*NRHS], the starting guess vectors.
//
//    Output, double EXACT[NROW*NRHS], the exact solution vectors.
//
{
//
//  Read the header block.
//
  hb_header_read ( input, title, key, totcrd, ptrcrd, indcrd, 
    valcrd, rhscrd, mxtype, nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, 
    valfmt, rhsfmt, rhstyp, nrhs, nrhsix );
//
//  Read the matrix structure.
//
  if ( 0 < ptrcrd )
  {
    if ( (*colptr) )
    {
      delete [] (*colptr);
    }
    (*colptr) = new int[*(ncol)+1];
  }

  if ( 0 < indcrd )
  {
    if ( (*rowind) )
    {
      delete [] (*rowind);
    }

    if ( (*mxtype)[2] == 'A' )
    {
      (*rowind) = new int[*nnzero];
    }
    else if ( (*mxtype)[2] == 'E' )
    {
      (*rowind) = new int[*neltvl];
    }
    else
    {
      cerr << "\n";
      cerr << "HB_FILE_READ - Fatal error!\n";
      cerr << "  Illegal value of MXTYPE character 3!\n";
      exit ( 1 );
    }

  }

  hb_structure_read ( input, *ncol, *mxtype, *nnzero, *neltvl, 
    *ptrcrd, *ptrfmt, *indcrd, *indfmt, *colptr, *rowind );
//
//  Read the matrix values.
//
  if ( 0 < valcrd )
  {
    if ( *values )
    {
      delete [] (*values);
    }

    if ( (*mxtype)[2] == 'A' )
    {
      *values = new double[*nnzero];
    }
    else if ( (*mxtype)[2] == 'E' )
    {
      *values = new double[*neltvl];
    }
    else
    {
      cerr << "\n";
      cerr << "HB_FILE_READ - Fatal error!\n";
      cerr << "  Illegal value of MXTYPE character 3!\n";
      exit ( 1 );
    }

    hb_values_read ( input, *valcrd, *mxtype, *nnzero, *neltvl, 
      *valfmt, *values );
  }
//
//  Read the right hand sides.
//
  if ( 0 < *rhscrd )
  {
    if ( (*rhstyp)[0] == 'F' )
    {
      if ( *rhsval )
      {
        delete [] *rhsval;
      }

      *rhsval = new double[(*nrow)*(*nrhs)];
    }
    else if ( (*rhstyp)[0] == 'M' && (*mxtype)[2] == 'A' )
    {
      if ( *rhsptr )
      {
        delete [] *rhsptr;
      }

      *rhsptr = new int[*nrhs+1];

      if ( *rhsind )
      {
        delete [] *rhsind;
      }

      *rhsind = new int[*nrhsix];

      if ( *rhsvec )
      {
        delete [] *rhsvec;
      }

      *rhsvec = new double[(*nrhsix)];
    }
    else if ( (*rhstyp)[0] == 'M' && (*mxtype)[2] == 'E' )
    {
      if ( *rhsval )
      {
        delete [] *rhsval;
      }

      *rhsval = new double[(*nnzero)*(*nrhs)];
    }
    else
    {
      cerr << "\n";
      cerr << "HB_FILE_READ - Fatal error!\n";
      cerr << "  Illegal combination of RHSTYP character 1\n";
      cerr << "  and MXTYPE character 3!\n";
      exit ( 1 );
    }

    hb_rhs_read ( input, *nrow, *nnzero, *nrhs, *nrhsix, 
      *rhscrd, *ptrfmt, *indfmt, *rhsfmt, *mxtype, *rhstyp, *rhsval, 
      *rhsind, *rhsptr, *rhsvec );
//
//  Read the starting guesses.
//
    if ( (*rhstyp)[1] == 'G' )
    {
      if ( *guess )
      {
        delete [] *guess;
      }

      *guess = new double[(*nrow)*(*nrhs)];

      hb_guess_read ( input, *nrow, *nrhs, *rhscrd, *rhsfmt, *rhstyp, *guess );
    }
//
//  Read the exact solutions.
//
    if ( (*rhstyp)[2] == 'X' )
    {
      if ( *exact )
      {
        delete [] *exact;
      }

      *exact = new double[(*nrow)*(*nrhs)];

      hb_exact_read ( input, *nrow, *nrhs, *rhscrd, *rhsfmt, *rhstyp, *exact );
    }
  }

  return;
}
//****************************************************************************80

void hb_file_write ( ofstream &output, char *title, char *key, int totcrd, 
  int ptrcrd, int indcrd, int valcrd, int rhscrd, char *mxtype, int nrow, 
  int ncol, int nnzero, int neltvl, char *ptrfmt, char *indfmt, char *valfmt, 
  char *rhsfmt, char *rhstyp, int nrhs, int nrhsix, int colptr[],
  int rowind[], double values[], double rhsval[], int rhsptr[], int rhsind[], 
  double rhsvec[], double guess[], double exact[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_FILE_WRITE writes an HB file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ofsteam &OUTPUT, the unit to which data is written.
//
//    Input, char *TITLE, a 72 character title for the matrix.
//
//    Input, char *KEY, an 8 character identifier for the matrix.
//
//    Input, int TOTCRD, the total number of lines of data.
//
//    Input, int PTRCRD, the number of input lines for pointers.
//
//    Input, int INDCRD, the number of input lines for row indices.
//
//    Input, int VALCRD, the number of input lines for numerical values.
//
//    Input, int RHSCRD, the number of input lines for right hand sides.
//
//    Input, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Input, int NROW, the number of rows or variables.
//
//    Input, int NCOL, the number of columns or elements.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Input, int NELTVL, the number of finite element matrix entries,
//    set to 0 in the case of assembled matrices.
//
//    Input, char *PTRFMT, the 16 character format for reading pointers.
//
//    Input, char *INDFMT, the 16 character format for reading indices.
//
//    Input, char *VALFMT, the 20 character format for reading values.
//
//    Input, char *RHSFMT, the 20 character format for reading values
//    of the right hand side.
//
//    Input, char *RHSTYP, the 3 character right hand side type.
//    First character is F for full storage or M for same as matrix.
//    Second character is G if starting "guess" vectors are supplied.
//    Third character is X if exact solution vectors are supplied.
//    Ignored if NRHS = 0.
//
//    Input, int NRHS, the number of right hand sides.
//
//    Input, int NRHSIX, the number of row indices (set to 0
//    in the case of unassembled matrices.)  Ignored if NRHS = 0.
//
//    Input, int COLPTR[NCOL+1], COLPTR(I) points to the location of
//    the first entry of column I in the sparse matrix structure.
//
//    If MXTYPE[2] == 'A':
//
//      Input, int ROWIND[NNZERO], the row index of each item.
//
//      Input, double VALUES[NNZERO], the nonzero values of the matrix.
//
//    If MXTYPE[2] == 'E':
//
//      Input, int ROWIND[NELTVL], the row index of each item.
//
//      Input, double VALUES[NELTVL], the nonzero values of the matrix.
//
//    If RHSTYP[0] == 'F':
//
//      Input, double RHSVAL[NROW*NRHS], contains NRHS dense right hand
//      side vectors.
//
//      Input, int RHSPTR[], is not used.
//
//      Input, int RHSIND[], is not used.
//
//      Input, double RHSVEC[], is not used.
//
//    If RHSTYP[0] = 'M' and MXTYPE[2] = 'A':
//
//      Input, double RHSVAL[], is not used.
//
//      Input, int RHSPTR[NRHS+1], RHSPTR(I) points to the location of
//      the first entry of right hand side I in the sparse right hand
//      side vector.
//
//      Input, int RHSIND[NRHSIX], indicates, for each entry of
//      RHSVEC, the corresponding row index.
//
//      Input, double RHSVEC[NRHSIX], contains the value of the right hand
//      side entries.
//
//    If RHSTYP[0] = 'M' and MXTYPE[2] = 'E':
//
//      Input, double RHSVAL[NNZERO*NRHS], contains NRHS unassembled
//      finite element vector right hand sides.
//
//      Input, int RHSPTR[], is not used.
//
//      Input, int RHSIND[], is not used.
//
//      Input, double RHSVEC[], is not used.
//
//    Input, double GUESS[NROW*NRHS], the starting guess vectors.
//
//    Input, double EXACT[NROW*NRHS], the exact solution vectors.
//
{
//
//  Write the header block.
//
  hb_header_write ( output, title, key, totcrd, ptrcrd, indcrd, 
    valcrd, rhscrd, mxtype, nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, 
    valfmt, rhsfmt, rhstyp, nrhs, nrhsix );
//
//  Write the matrix structure.
//
  hb_structure_write ( output, ncol, mxtype, nnzero, neltvl, 
    ptrfmt, indfmt, colptr, rowind );
//
//  Write the matrix values.
//
  hb_values_write ( output, valcrd, mxtype, nnzero, neltvl, 
    valfmt, values );
//
//  Write the right hand sides.
//
  hb_rhs_write ( output, nrow, nnzero, nrhs, nrhsix, 
    rhscrd, ptrfmt, indfmt, rhsfmt, mxtype, rhstyp, rhsval, 
    rhsind, rhsptr, rhsvec );
//
//  Write the starting guesses.
//
  hb_guess_write ( output, nrow, nrhs, rhscrd, rhsfmt, rhstyp, guess );
//
//  Write the exact solutions.
//
  hb_exact_write ( output, nrow, nrhs, rhscrd, rhsfmt, rhstyp, exact );

  return;
}
//****************************************************************************80

void hb_guess_read ( ifstream &input, int nrow, int nrhs, int rhscrd, 
  char *rhsfmt, char *rhstyp, double guess[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_GUESS_READ reads the starting guess vectors in an HB file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ifstream &INPUT, the unit from which data is read.
//
//    Input, int NROW, the number of rows or variables.
//
//    Input, int NRHS, the number of right hand sides.
//
//    Input, int RHSCRD, the number of lines in the file for
//    right hand sides.
//
//    Input, char *RHSFMT, the 20 character format for reading values
//    of the right hand side.
//
//    Input, char *RHSTYP, the 3 character right hand side type.
//    First character is F for full storage or M for same as matrix.
//    Second character is G if starting "guess" vectors are supplied.
//    Third character is X if exact solution vectors are supplied.
//    Ignored if NRHS = 0.
//
//    Output, double GUESS[NROW*NRHS], the starting guess vectors.
//
{
  char code;
  int i;
  int j;
  int jhi;
  int jlo;
  int khi;
  int klo;
  char line[255];
  int line_num;
  int m;
  int r;
  char *s;
  int w;

  if ( 0 < rhscrd )
  {
    if ( rhstyp[1] == 'G' )
    {
      s_to_format ( rhsfmt, &r, &code, &w, &m );

      line_num = 1 + ( nrow * nrhs - 1 ) / r;

      jhi = 0;
      for ( i = 1; i <= line_num; i++ )
      {
        input.getline ( line, sizeof ( line ) );
        jlo = jhi + 1;
        jhi = i4_min ( jlo + r - 1, nrow * nrhs );

        khi = 0;
        for ( j = jlo; j <= jhi; j++ )
        {
          klo = khi + 1;
          khi = i4_min ( klo + w - 1, strlen ( line ) );
          s = s_substring ( line, klo, khi );
          guess[j-1] = atof ( s );
        }
      }
    }
  }

  return;
}
//****************************************************************************80

void hb_guess_write ( ofstream &output, int nrow, int nrhs, int rhscrd, 
  char *rhsfmt, char *rhstyp, double guess[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_GUESS_WRITE writes the starting guess vectors to an HB file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ofstream &OUTPUT, the unit to which data is written.
//
//    Input, int NROW, the number of rows or variables.
//
//    Input, int NRHS, the number of right hand sides.
//
//    Input, int RHSCRD, the number of lines in the file for
//    right hand sides.
//
//    Input, char *RHSFMT, the 20 character format for reading values
//    of the right hand side.
//
//    Input, char *RHSTYP, the 3 character right hand side type.
//    First character is F for full storage or M for same as matrix.
//    Second character is G if starting "guess" vectors are supplied.
//    Third character is X if exact solution vectors are supplied.
//    Ignored if NRHS = 0.
//
//    Input, double GUESS[NROW*NRHS], the starting guess vectors.
//
{
  char code;
  int i;
  int j;
  int jhi;
  int jlo;
  int line_num;
  int m;
  int r;
  int w;

  if ( 0 < rhscrd )
  {
    if ( rhstyp[1] == 'G' )
    {
      s_to_format ( rhsfmt, &r, &code, &w, &m );
      line_num = 1 + ( nrow * nrhs - 1 ) / r;

      jhi = 0;
      for ( i = 1; i <= line_num; i++ )
      {
        jlo = jhi + 1;
        jhi = i4_min ( jlo + r - 1, nrow * nrhs );

        for ( j = jlo; j <= jhi; j++ )
        {
          output << setw(w) << guess[j-1];
        }
        output << "\n";
      }
    }
  }
  return;
}
//****************************************************************************80

void hb_header_print ( char *title, char *key, int totcrd, int ptrcrd, 
  int indcrd, int valcrd, int rhscrd, char *mxtype, int nrow, int ncol, 
  int nnzero, int neltvl, char *ptrfmt, char *indfmt, char *valfmt, 
  char *rhsfmt, char *rhstyp, int nrhs, int nrhsix )

//****************************************************************************80
//
//  Purpose:
//
//    HB_HEADER_PRINT prints the header of an HB file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 April 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, char *TITLE, a 72 character title for the matrix.
//
//    Input, char *KEY, an 8 character identifier for the matrix.
//
//    Input, int TOTCRD, the total number of lines of data.
//
//    Input, int PTRCRD, the number of input lines for pointers.
//
//    Input, int INDCRD, the number of input lines for row indices.
//
//    Input, int VALCRD, the number of input lines for numerical values.
//
//    Input, int RHSCRD, the number of input lines for right hand sides.
//
//    Input, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Input, int NROW, the number of rows or variables.
//
//    Input, int NCOL, the number of columns or elements.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Input, int NELTVL, the number of finite element matrix entries,
//    set to 0 in the case of assembled matrices.
//
//    Input, char *PTRFMT, the 16 character format for reading pointers.
//
//    Input, char *INDFMT, the 16 character format for reading indices.
//
//    Input, char *VALFMT, the 20 character format for reading values.
//
//    Input, char *RHSFMT, the 20 character format for reading values
//    of the right hand side.
//
//    Input, char *RHSTYP, the 3 character right hand side type.
//    First character is F for full storage or M for same as matrix.
//    Second character is G if starting "guess" vectors are supplied.
//    Third character is X if exact solution vectors are supplied.
//
//    Input, int NRHS, the number of right hand sides.
//
//    Input, int NRHSIX, the number of entries of storage for right
//    hand side values, in the case where RHSTYP[0] = 'M' and
//    MXTYPE[2] = 'A'.
//
{
  cout << "\n";
  cout << title << "\n";
  cout << "\n";
  cout << "  TOTCRD = " << totcrd << "\n";
  cout << "  PTRCRD = " << ptrcrd << "\n";
  cout << "  INDCRD = " << indcrd << "\n";
  cout << "  VALCRD = " << valcrd << "\n";
  cout << "  RHSCRD = " << rhscrd << "\n";
  cout << "\n";
  cout << "  KEY =    '" << key    << "'.\n";
  cout << "  MXTYPE = '" << mxtype << "'.\n";
  cout << "  RHSTYP = '" << rhstyp << "'.\n";
  cout << "\n";
  cout << "  NROW =   " << nrow   << "\n";
  cout << "  NCOL =   " << ncol   << "\n";
  cout << "  NNZERO = " << nnzero << "\n";
  cout << "  NELTVL = " << neltvl << "\n";
  cout << "  NRHS =   " << nrhs   << "\n";
  cout << "  NRHSIX = " << nrhsix << "\n";
  cout << "\n";
  cout << "  PTRFMT = '" << ptrfmt << "'.\n";
  cout << "  INDFMT = '" << indfmt << "'.\n";
  cout << "  VALFMT = '" << valfmt << "'.\n";
  cout << "  RHSFMT = '" << rhsfmt << "'.\n";

  return;
}
//****************************************************************************80

void hb_header_read ( ifstream &input, char **title, char **key, int *totcrd, 
  int *ptrcrd, int *indcrd, int *valcrd, int *rhscrd, char **mxtype, int *nrow, 
  int *ncol, int *nnzero, int *neltvl, char **ptrfmt, char **indfmt, char **valfmt, 
  char **rhsfmt, char **rhstyp, int *nrhs, int *nrhsix )

//****************************************************************************80
//
//  Purpose:
//
//    HB_HEADER_READ reads the header of an HB file.
//
//  Discussion:
//
//    The user should already have opened the file, and positioned it
//    to the first record.
//
//    Thanks to Alexander Zilyakov for pointing out that the line 
//      if ( 0 < rhscrd )
//    should be 
//      if ( 0 < *rhscrd )
//    15 August 2016.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    15 August 2016
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ifstream &INPUT, the unit from which data is read.
//
//    Output, char *TITLE, a 72 character title for the matrix.
//
//    Output, char *KEY, an 8 character identifier for the matrix.
//
//    Output, int *TOTCRD, the total number of lines of data.
//
//    Output, int *PTRCRD, the number of input lines for pointers.
//
//    Output, int *INDCRD, the number of input lines for row indices.
//
//    Output, int *VALCRD, the number of input lines for numerical values.
//
//    Output, int *RHSCRD, the number of input lines for right hand sides.
//
//    Output, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Output, int *NROW, the number of rows or variables.
//
//    Output, int *NCOL, the number of columns or elements.
//
//    Output, int *NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Output, int *NELTVL, the number of finite element matrix entries,
//    set to 0 in the case of assembled matrices.
//
//    Output, char *PTRFMT, the 16 character format for reading pointers.
//
//    Output, char *INDFMT, the 16 character format for reading indices.
//
//    Output, char *VALFMT, the 20 character format for reading values.
//
//    Output, char *RHSFMT, the 20 character format for reading values
//    of the right hand side.
//
//    Output, char *RHSTYP, the 3 character right hand side type.
//    First character is F for full storage or M for same as matrix.
//    Second character is G if starting "guess" vectors are supplied.
//    Third character is X if exact solution vectors are supplied.
//
//    Output, int *NRHS, the number of right hand sides.
//
//    Output, int *NRHSIX, the number of entries of storage for right
//    hand side values, in the case where RHSTYP[0] = 'M' and
//    MXTYPE[2] = 'A'.
//
{
  char *field;
  char line[255];
//
//  Read line 1.
//
  input.getline ( line, sizeof ( line ) );

  if ( input.eof() )
  {
    cerr << "\n";
    cerr << "HB_HEADER_READ - Fatal error!\n";
    cerr << "  I/O error reading header line 1.\n";
    exit ( 1 );
  }

  *title = s_substring ( line, 1, 72 );
  s_trim ( *title );

  *key = s_substring ( line, 73, 80 );
  s_trim ( *key );
//
//  Read line 2.
//
  input.getline ( line, sizeof ( line ) );

  if ( input.eof() )
  {
    cerr << "\n";
    cerr << "HB_HEADER_READ - Fatal error!\n";
    cerr << "  I/O error reading header line 2.\n";
    exit ( 1 );
  }

  field = s_substring ( line,  1, 14 );
  *totcrd = atoi ( field );

  field = s_substring ( line, 15, 28 );
  *ptrcrd = atoi ( field );

  field = s_substring ( line, 29, 42 );
  *indcrd = atoi ( field );

  field = s_substring ( line, 43, 56 );
  *valcrd = atoi ( field );

  field = s_substring ( line, 57, 70 );
  *rhscrd = atoi ( field );
//
//  Read line 3.
//
  input.getline ( line, sizeof ( line ) );
  
  if ( input.eof() )
  {
    cerr << "\n";
    cerr << "HB_HEADER_READ - Fatal error!\n";
    cerr << "  I/O error reading header line 3.\n";
    exit ( 1 );
  }

  *mxtype = s_substring ( line, 1, 3 );
  s_trim ( *mxtype );

  field = s_substring ( line, 15, 28 );
  *nrow = atoi ( field );

  field = s_substring ( line, 29, 42 );
  *ncol = atoi ( field );

  field = s_substring ( line, 43, 56 );
  *nnzero = atoi ( field );

  field = s_substring ( line, 57, 70 );
  *neltvl = atoi ( field );
//
//  Read line 4.
//
  input.getline ( line, sizeof ( line ) );
  
  if ( input.eof() )
  {
    cerr << "\n";
    cerr << "HB_HEADER_READ - Fatal error!\n";
    cerr << "  I/O error reading header line 4.\n";
    exit ( 1 );
  }

  *ptrfmt = s_substring ( line,  1, 16 );
  s_trim ( *ptrfmt );

  *indfmt = s_substring ( line, 17, 32 );
  s_trim ( *indfmt );

  *valfmt = s_substring ( line, 33, 52 );
  s_trim ( *valfmt );

  *rhsfmt = s_substring ( line, 53, 72 );
  s_trim ( *rhsfmt );
//
//  Read line 5.
//
  if ( 0 < *rhscrd )
  {
    input.getline ( line, sizeof ( line ) );
  
    if ( input.eof ( ) )
    {
      cerr << "\n";
      cerr << "HB_HEADER_READ - Fatal error!\n";
      cerr << "  I/O error reading header line 5.\n";
      exit ( 1 );
    }

    *rhstyp = s_substring ( line, 1, 3 );
    s_trim ( *rhstyp );

    field = s_substring ( line, 15, 28 );
    *nrhs = atoi ( field );

    field = s_substring ( line, 29, 42 );
    *nrhsix = atoi ( field );
  }
  else
  {
    *rhstyp = NULL;
    *nrhs = 0;
    *nrhsix = 0;
  }

  return;
}
//****************************************************************************80

void hb_header_write ( ofstream &output, char *title, char *key, int totcrd, 
  int ptrcrd, int indcrd, int valcrd, int rhscrd, char *mxtype, int nrow, 
  int ncol, int nnzero, int neltvl, char *ptrfmt, char *indfmt, char *valfmt, 
  char *rhsfmt, char *rhstyp, int nrhs, int nrhsix )

//****************************************************************************80
//
//  Purpose:
//
//    HB_HEADER_WRITE writes the header of an HB file.
//
//  Discussion:
//
//    The user should already have opened the file, and positioned it
//    to the first record.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 April 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ofstream &OUTPUT, the unit to which data is written.
//
//    Input, char *TITLE, a 72 character title for the matrix.
//
//    Input, char *KEY, an 8 character identifier for the matrix.
//
//    Input, int TOTCRD, the total number of lines of data.
//
//    Input, int PTRCRD, the number of input lines for pointers.
//
//    Input, int INDCRD, the number of input lines for row indices.
//
//    Input, int VALCRD, the number of input lines for numerical values.
//
//    Input, int RHSCRD, the number of input lines for right hand sides.
//
//    Input, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Input, int NROW, the number of rows or variables.
//
//    Input, int NCOL, the number of columns or elements.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Input, int NELTVL, the number of finite element matrix entries,
//    set to 0 in the case of assembled matrices.
//
//    Input, char *PTRFMT, the 16 character format for reading pointers.
//
//    Input, char *INDFMT, the 16 character format for reading indices.
//
//    Input, char *VALFMT, the 20 character format for reading values.
//
//    Input, char *RHSFMT, the 20 character format for reading values
//    of the right hand side.
//
//    Input, char *RHSTYP, the 3 character right hand side type.
//    First character is F for full storage or M for same as matrix.
//    Second character is G if starting "guess" vectors are supplied.
//    Third character is X if exact solution vectors are supplied.
//    Ignored if NRHS = 0.
//
//    Input, int NRHS, the number of right hand sides.
//
//    Input, int NRHSIX, the number of row indices (set to 0
//    in the case of unassembled matrices.)  Ignored if NRHS = 0.
//
{
  output << setiosflags(ios::left)
         << setw(72) << title
         << setw(8)  << key 
         << resetiosflags(ios::left) << "\n";

  output << setw(14) << totcrd
         << setw(14) << ptrcrd
         << setw(14) << indcrd
         << setw(14) << valcrd
         << setw(14) << rhscrd << "\n";

  output << setiosflags(ios::left)
         << setw(3)  << mxtype
         << resetiosflags(ios::left)
         << "           "
         << setw(14) << nrow
         << setw(14) << ncol
         << setw(14) << nnzero
         << setw(14) << neltvl << "\n";

  output << setiosflags(ios::left) 
         << setw(16) << ptrfmt
         << setw(16) << indfmt
         << setw(20) << valfmt
         << setw(20) << rhsfmt 
         << resetiosflags(ios::left) << "\n";

  if ( 0 < rhscrd )
  {
    output << setiosflags(ios::left) 
           << setw(3)  << rhstyp
           << resetiosflags(ios::left) 
           << "           "
           << setw(14) << nrhs
           << setw(14) << nrhsix << "\n";
  }

  return;
}
//****************************************************************************80

double *hb_matvec_a_mem ( int nrow, int ncol, int nnzero, int nrhs, 
  int colptr[], int rowind[], double values[], double exact[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_MATVEC_A_MEM multiplies an assembled Harwell Boeing matrix times a vector.
//
//  Discussion:
//
//    In this "A_MEM" version of the routine, the matrix is assumed to be in
//    "assembled" form, and all the data is assumed to be small enough
//    to reside completely in memory; the matrix and multiplicand vectors
//    are assumed to have been read into memory before this routine is called.
//
//    It is assumed that MXTYPE(3:3) = 'A', that is, that the matrix is
//    stored in the "assembled" format.
//
//    Also, the storage used for the vectors X and the products A*X
//    corresponds to RHSTYP(1:1) = 'F', that is, the "full" storage mode
//    for vectors.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, int NROW, the number of rows or variables.
//
//    Input, int NCOL, the number of columns or elements.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  
//
//    Input, int NRHS, the number of right hand sides.
//
//    Input, int COLPTR[NCOL+1], COLPTR(I) points to the location of
//    the first entry of column I in the sparse matrix structure.
//
//    Input, int ROWIND[NNZERO], the row index of each item.
//
//    Input, double VALUES[NNZERO], the nonzero values of the matrix.
//
//    Input, double EXACT[NCOL*NRHS], contains NRHS dense vectors.
//
//    Output, double HB_MATVEC_A_MEM[NROW*NRHS], the product vectors A*X.
//
{
  int column;
  int k;
  double *rhsval;
  int rhs;
  int row;

  rhsval = new double[nrow*nrhs];
//
//  Zero out the result vectors.
//
  for ( rhs = 1; rhs <= nrhs; rhs++ )
  {
    for ( row = 1; row <= nrow; row++ )
    {
      rhsval[row-1+(rhs-1)*nrow] = 0.0E+00;
    }
  }
//
//  For each column J of the matrix,
//
  for ( column = 1; column <= ncol; column++ )
  {
//
//  For nonzero entry K
//
    for ( k = colptr[column-1]; k <= colptr[column]-1; k++ )
    {
      row = rowind[k-1];
//
//  For each right hand side vector:
//
//    B(I,1:NRHS) = B(I,1:NRHS) + A(I,J) * X(J,1:NRHS)
//
      for ( rhs = 1; rhs <= nrhs; rhs++ )
      {
        rhsval[row-1+(rhs-1)*nrow] = rhsval[row-1+(rhs-1)*nrow] 
          + values[k-1] * exact[column-1+(rhs-1)*ncol];
      }
    }
  }

  return rhsval;
}
//****************************************************************************80

void hb_rhs_read ( ifstream &input, int nrow, int nnzero, int nrhs, int nrhsix, 
  int rhscrd, char *ptrfmt, char *indfmt, char *rhsfmt, char *mxtype, 
  char *rhstyp, double rhsval[], int rhsind[], int rhsptr[], double rhsvec[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_RHS_READ reads the right hand side information in an HB file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ifstream &INPUT, the unit from which data is read.
//
//    Input, int NROW, the number of rows or variables.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Input, int NRHS, the number of right hand sides.
//
//    Input, int NRHSIX, the number of entries of storage for right
//    hand side values, in the case where RHSTYP[0] = 'M' and
//    MXTYPE[2] = 'A'.
//
//    Input, int RHSCRD, the number of lines in the file for
//    right hand sides.
//
//    Input, char *PTRFMT, the 16 character format for reading pointers.
//
//    Input, char *INDFMT, the 16 character format for reading indices.
//
//    Input, char *RHSFMT, the 20 character format for reading values
//    of the right hand side.
//
//    Input, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Input, char *RHSTYP, the 3 character right hand side type.
//    First character is F for full storage or M for same as matrix.
//    Second character is G if starting "guess" vectors are supplied.
//    Third character is X if exact solution vectors are supplied.
//    Ignored if NRHS = 0.
//
//    If RHSTYP[0] == 'F':
//
//      Output, double RHSVAL[NROW*NRHS], contains NRHS dense right hand
//      side vectors.
//
//      Output, int RHSPTR[], is not used.
//
//      Output, int RHSIND[], is not used.
//
//      Output, int RHSVEC[], is not used.
//
//    If RHSTYP[0] = 'M' and MXTYPE[2] = 'A':
//
//      Output, double RHSVAL[], is not used.
//
//      Output, int RHSPTR[NRHS+1], RHSPTR[I-1] points to the location of
//      the first entry of right hand side I in the sparse right hand
//      side vector.
//
//      Output, int RHSIND[NRHSIX], indicates, for each entry of
//      RHSVEC, the corresponding row index.
//
//      Output, double RHSVEC[NRHSIX], contains the value of the right hand
//      side entries.
//
//    If RHSTYP[0] = 'M' and MXTYPE[2] = 'E':
//
//      Output, double RHSVAL[NNZERO*NRHS], contains NRHS unassembled
//      finite element vector right hand sides.
//
//      Output, int RHSPTR[], is not used.
//
//      Output, int RHSIND[], is not used.
//
//      Output, double RHSVEC[], is not used.
//
{
  char code;
  int i;
  int j;
  int jhi;
  int jlo;
  int khi;
  int klo;
  char line[255];
  int line_num;
  int m;
  int r;
  char *s;
  int w;
//
//  Read the right hand sides.
//    case F                             = "full" or "dense";
//    case not F + matrix storage is "A" = sparse pointer RHS
//    case not F + matrix storage is "E" = finite element RHS
//
  if ( 0 < rhscrd )
  {
//
//  Dense right hand sides:
//
    if ( rhstyp[0] == 'F' )
    {
      s_to_format ( rhsfmt, &r, &code, &w, &m );

      line_num = 1 + ( nrow * nrhs - 1 ) / r;

      jhi = 0;
      for ( i = 1; i <= line_num; i++ )
      {
        input.getline ( line, sizeof ( line ) );
        jlo = jhi + 1;
        jhi = i4_min ( jlo + r - 1, nrow * nrhs );

        khi = 0;
        for ( j = jlo; j <= jhi; j++ )
        {
          klo = khi + 1;
          khi = i4_min ( klo + w - 1, strlen ( line ) );
          s = s_substring ( line, klo, khi );
          rhsval[j-1] = atof ( s );
        }
      }
    }
//
//  Sparse right-hand sides stored like the matrix.
//  Read pointer array, indices, and values.
//
    else if ( rhstyp[0] == 'M' )
    {
      if ( mxtype[2] == 'A' )
      {
        s_to_format ( ptrfmt, &r, &code, &w, &m );

        line_num = 1 + ( nrhs + 1 - 1 ) / r;

        jhi = 0;
        for ( i = 1; i <= line_num; i++ )
        {
          input.getline ( line, sizeof ( line ) );
          jlo = jhi + 1;
          jhi = i4_min ( jlo + r - 1, nrhs + 1 );
 
          khi = 0;
          for ( j = jlo; j <= jhi; j++ )
          { 
            klo = khi + 1;
            khi = i4_min ( klo + w - 1, strlen ( line ) );
            s = s_substring ( line, klo, khi );
            rhsptr[j-1] = atoi ( s );
          }
        }

        s_to_format ( indfmt, &r, &code, &w, &m );

        line_num = 1 + ( nrhsix - 1 ) / r;

        jhi = 0;
        for ( i = 1; i <= line_num; i++ )
        {
          input.getline ( line, sizeof ( line ) );
          jlo = jhi + 1;
          jhi = i4_min ( jlo + r - 1, nnzero );
 
          khi = 0;
          for ( j = jlo; j <= jhi; j++ )
          { 
            klo = khi + 1;
            khi = i4_min ( klo + w - 1, strlen ( line ) );
            s = s_substring ( line, klo, khi );
            rhsind[j-1] = atoi ( s );
          }
        }

        s_to_format ( rhsfmt, &r, &code, &w, &m );

        line_num = 1 + ( nrhsix - 1 ) / r;

        jhi = 0;
        for ( i = 1; i <= line_num; i++ )
        {
          input.getline ( line, sizeof ( line ) );
          jlo = jhi + 1;
          jhi = i4_min ( jlo + r - 1, nrhsix );

          khi = 0;
          for ( j = jlo; j <= jhi; j++ )
          {
            klo = khi + 1;
            khi = i4_min ( klo + w - 1, strlen ( line ) );
            s = s_substring ( line, klo, khi );
            rhsvec[j-1] = atof ( s );
          }
        }
      }
//
//  Sparse right hand sides in finite element format.
//
      else if ( mxtype[2] == 'E' )
      {
        s_to_format ( rhsfmt, &r, &code, &w, &m );

        line_num = 1 + ( nnzero * nrhs - 1 ) / r;

        jhi = 0;
        for ( i = 1; i <= line_num; i++ )
        {
          input.getline ( line, sizeof ( line ) );
          jlo = jhi + 1;
          jhi = i4_min ( jlo + r - 1, nnzero * nrhs );

          khi = 0;
          for ( j = jlo; j <= jhi; j++ )
          {
            klo = khi + 1;
            khi = i4_min ( klo + w - 1, strlen ( line ) );
            s = s_substring ( line, klo, khi );
            rhsval[j-1] = atof ( s );
          }
        }
      }
      else
      {
        cerr << "\n";
        cerr << "HB_RHS_READ - Fatal error!\n";
        cerr << "  Illegal value of MXTYPE character 3!\n";
        exit ( 1 );
      }

    }
//
//  0 < RHSCRD, but RHSTYP not recognized.
//
    else
    {
      cerr << "\n";
      cerr << "HB_RHS_READ - Fatal error!\n";
      cerr << "  Illegal value of RHSTYP character 1!\n";
      exit ( 1 );
    }

  }

  return;
}
//****************************************************************************80

void hb_rhs_write ( ofstream &output, int nrow, int nnzero, int nrhs, int nrhsix, 
  int rhscrd, char *ptrfmt, char *indfmt, char *rhsfmt, char *mxtype, 
  char *rhstyp, double rhsval[], int rhsind[], int rhsptr[], double rhsvec[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_RHS_WRITE writes the right hand side information to an HB file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ofstream &OUTPUT, the unit to which data is written.
//
//    Input, int NROW, the number of rows or variables.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Input, int NRHS, the number of right hand sides.
//
//    Input, int NRHSIX, the number of entries of storage for right
//    hand side values, in the case where RHSTYP[0] = 'M' and
//    MXTYPE[2] = 'A'.
//
//    Input, int RHSCRD, the number of lines in the file for
//    right hand sides.
//
//    Input, char *PTRFMT, the 16 character format for reading pointers.
//
//    Input, char *INDFMT, the 16 character format for reading indices.
//
//    Input, char *RHSFMT, the 20 character format for reading values
//    of the right hand side.
//
//    Input, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Input, char *RHSTYP, the 3 character right hand side type.
//    First character is F for full storage or M for same as matrix.
//    Second character is G if starting "guess" vectors are supplied.
//    Third character is X if exact solution vectors are supplied.
//    Ignored if NRHS = 0.
//
//    If RHSTYP[0] == 'F':
//
//      Input, double RHSVAL[NROW*NRHS], contains NRHS dense right hand
//      side vectors.
//
//      Input, int RHSPTR[], is not used.
//
//      Input, int RHSIND[], is not used.
//
//      Input, double RHSVEC[], is not used.
//
//    If RHSTYP[0] = 'M' and MXTYPE[2] = 'A':
//
//      Input, double RHSVAL[], is not used.
//
//      Input, int RHSPTR[NRHS+1], RHSPTR(I) points to the location of
//      the first entry of right hand side I in the sparse right hand
//      side vector.
//
//      Input, int RHSIND[NRHSIX], indicates, for each entry of
//      RHSVEC, the corresponding row index.
//
//      Input, double RHSVEC[NRHSIX], contains the value of the right hand
//      side entries.
//
//    If RHSTYP[0] = 'M' and MXTYPE[2] = 'E':
//
//      Input, double RHSVAL[NNZERO*NRHS], contains NRHS unassembled
//      finite element vector right hand sides.
//
//      Input, int RHSPTR[], is not used.
//
//      Input, int RHSIND[], is not used.
//
//      Input, double RHSVEC[], is not used.
//
{
  char code;
  int i;
  int j;
  int jhi;
  int jlo;
  int line_num;
  int m;
  int r;
  int w;
//
//  Read the right hand sides.
//    case F                             = "full" or "dense";
//    case not F + matrix storage is "A" = sparse pointer RHS
//    case not F + matrix storage is "E" = finite element RHS
//
  if ( 0 < rhscrd )
  {
//
//  Dense right hand sides:
//
    if ( rhstyp[0] == 'F' )
    {
      s_to_format ( rhsfmt, &r, &code, &w, &m );
      line_num = 1 + ( nrow * nrhs - 1 ) / r;

      jhi = 0;
      for ( i = 1; i <= line_num; i++ )
      {
        jlo = jhi + 1;
        jhi = i4_min ( jlo + r - 1, nrow * nrhs );

        for ( j = jlo; j <= jhi; j++ )
        {
          output << setw(w) << rhsval[j-1];
        }
        output << "\n";
      }
    }
//
//  Sparse right-hand sides stored like the matrix.
//  Read pointer array, indices, and values.
//
    else if ( rhstyp[0] == 'M' )
    {
      if ( mxtype[2] == 'A' )
      {
        s_to_format ( ptrfmt, &r, &code, &w, &m );
        line_num = 1 + ( nrhs + 1 - 1 ) / r;

        jhi = 0;
        for ( i = 1; i <= line_num; i++ )
        {
          jlo = jhi + 1;
          jhi = i4_min ( jlo + r - 1, nrhs + 1 );

          for ( j = jlo; j <= jhi; j++ )
          {
            output << setw(w) << rhsptr[j-1];
          }
          output << "\n";
        }

        s_to_format ( indfmt, &r, &code, &w, &m );
        line_num = 1 + ( nrhsix - 1 ) / r;

        jhi = 0;
        for ( i = 1; i <= line_num; i++ )
        {
          jlo = jhi + 1;
          jhi = i4_min ( jlo + r - 1, nrhsix );

          for ( j = jlo; j <= jhi; j++ )
          {
            output << setw(w) << rhsind[j-1];
          }
          output << "\n";
        }

        s_to_format ( rhsfmt, &r, &code, &w, &m );
        line_num = 1 + ( nrhsix - 1 ) / r;

        jhi = 0;
        for ( i = 1; i <= line_num; i++ )
        {
          jlo = jhi + 1;
          jhi = i4_min ( jlo + r - 1, nrhsix );

          for ( j = jlo; j <= jhi; j++ )
          {
            output << setw(w) << rhsvec[j-1];
          }
          output << "\n";
        }
      }
//
//  Sparse right hand sides in finite element format.
//
      else if ( mxtype[2] == 'E' )
      {
        s_to_format ( rhsfmt, &r, &code, &w, &m );
        line_num = 1 + ( nnzero * nrhs - 1 ) / r;

        jhi = 0;
        for ( i = 1; i <= line_num; i++ )
        {
          jlo = jhi + 1;
          jhi = i4_min ( jlo + r - 1, nnzero * nrhs );

          for ( j = jlo; j <= jhi; j++ )
          {
            output << setw(w) << rhsval[j-1];
          }
          output << "\n";
        }
      }
      else
      {
        cerr << "\n";
        cerr << "HB_RHS_WRITE - Fatal error!\n";
        cerr << "  Illegal value of MXTYPE character 3!\n";
        exit ( 1 );
      }
    }
  }
  else
  {
    cerr << "\n";
    cerr << "HB_RHS_WRITE - Fatal error!\n";
    cerr << "  Illegal value of RHSTYP character 1!\n";
    exit ( 1 );
  }

  return;
}
//****************************************************************************80

void hb_structure_print ( int ncol, char *mxtype, int nnzero, int neltvl, 
  int colptr[], int rowind[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_STRUCTURE_PRINT prints the structure of an HB matrix.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 April 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, int NCOL, the number of columns.
//
//    Input, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Input, int NELTVL, the number of finite element matrix entries,
//    set to 0 in the case of assembled matrices.
//
//    Input, int COLPTR[NCOL+1], COLPTR[I-1] points to the location of
//    the first entry of column I in the sparse matrix structure.
//
//    If MXTYPE[2] == 'A':
//
//      Input, int ROWIND[NNZERO], the row index of each item.
//
//    If MXTYPE[2] == 'F':
//
//      Input, int ROWIND[NELTVL], the row index of each item.
//
{
  int j;
  int k;
  int khi;
  int klo;

  if ( mxtype[2] == 'A' )
  {
    cout << "\n";
    cout << "Column Begin   End   ----------------------------------------\n";
    cout << "\n";
    for ( j = 1; j <= i4_min ( ncol, 10 ); j++ )
    {
      if ( colptr[j]-1 < colptr[j-1] )
      {
        cout << setw(6) << j << "   EMPTY\n";
      }
      else
      {
        for ( klo = colptr[j-1]; klo <= colptr[j]-1; klo = klo + 10 )
        {
          khi = i4_min ( klo + 9, colptr[j]-1 );
          if ( klo == colptr[j-1] )
          {
            cout << setw(6) << j
                 << setw(6) << colptr[j-1]
                 << setw(6) << colptr[j]-1 << "   ";
          }
          for ( k = klo; k <= khi; k++ )
          {
            cout << setw(4) << rowind[k-1];
          }
          cout << "\n";
        }
      }
    }
    cout << "                     ----------------------------------------\n";
  }
  else if ( mxtype[2] == 'E' )
  {
    cout << "\n";
    cout << "Column Begin   End   ----------------------------------------\n";
    cout << "                     ----------------------------------------\n";

    cout << "\n";
    cout << "  I haven't thought about how to print an\n";
    cout << "  unassembled matrix yet!\n";
  }
  else
  {
    cerr << "\n";
    cerr << "HB_STRUCTURE_PRINT - Fatal error!\n";
    cerr << "  Illegal value of MXTYPE character #3 = " << mxtype[2] << "\n";;
    exit ( 1 );
  }

  return;
}
//****************************************************************************80

void hb_structure_read ( ifstream &input, int ncol, char *mxtype, int nnzero, 
  int neltvl, int ptrcrd, char *ptrfmt, int indcrd, char *indfmt, 
  int colptr[], int rowind[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_STRUCTURE_READ reads the structure of an HB matrix.
//
//  Discussion:
//
//    The user should already have opened the file, and positioned it
//    to just after the header records.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 April 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ifstream &INPUT, the unit from which data is read.
//
//    Input, int NCOL, the number of columns.
//
//    Input, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Input, int NELTVL, the number of finite element matrix entries,
//    set to 0 in the case of assembled matrices.
//
//    Input, int PTRCRD, the number of input lines for pointers.
//
//    Input, char *PTRFMT, the 16 character format for reading pointers.
//
//    Input, int INDCRD, the number of input lines for indices.
//
//    Input, char *INDFMT, the 16 character format for reading indices.
//
//    Output, int COLPTR[NCOL+1], COLPTR[I-1] points to the location of
//    the first entry of column I in the sparse matrix structure.
//
//    If MXTYPE[2] == 'A':
//
//      Output, int ROWIND[NNZERO], the row index of each item.
//
//    If MXTYPE[2] == 'F':
//
//      Output, int ROWIND[NELTVL], the row index of each item.
//
{
  char code;
  int i;
  int j;
  int jhi;
  int jlo;
  int khi;
  int klo;
  char line[255];
  int line_num;
  int m;
  int number;
  int r;
  char *s;
  int w;

  s_to_format ( ptrfmt, &r, &code, &w, &m );

  if ( mxtype[2] == 'A' )
  {
    line_num = 1 + ( ( ncol + 1 ) - 1 ) / r;
  }
  else
  {
    line_num = 1 + ( ( ncol     ) - 1 ) / r;
  }

  jhi = 0;
  for ( i = 1; i <= line_num; i++ )
  {
    input.getline ( line, sizeof ( line ) );
    jlo = jhi + 1;
    if ( mxtype[2] == 'A' )
    {
      jhi = i4_min ( jlo + r - 1, ncol + 1 );
    }
    else
    {
      jhi = i4_min ( jlo + r - 1, ncol     );
    }
    khi = 0;
    for ( j = jlo; j <= jhi; j++ )
    {
      klo = khi + 1;
      khi = i4_min ( klo + w - 1, strlen ( line ) );
      s = s_substring ( line, klo, khi );
      colptr[j-1] = atoi ( s );
    }
  }

  if ( mxtype[2] == 'A' )
  {
    s_to_format ( indfmt, &r, &code, &w, &m );

    line_num = 1 + ( nnzero - 1 ) / r;

    jhi = 0;
    for ( i = 1; i <= line_num; i++ )
    {
      input.getline ( line, sizeof ( line ) );
      jlo = jhi + 1;
      jhi = i4_min ( jlo + r - 1, nnzero );
 
      khi = 0;
      for ( j = jlo; j <= jhi; j++ )
      { 
        klo = khi + 1;
        khi = i4_min ( klo + w - 1, strlen ( line ) );
        s = s_substring ( line, klo, khi );
        rowind[j-1] = atoi ( s );
      }
    }
  }
  else if ( mxtype[2] == 'E' )
  {
    s_to_format ( indfmt, &r, &code, &w, &m );

    number = colptr[ncol-1] - colptr[0];
    line_num = 1 + ( number - 1 ) / r;

    jhi = 0;
    for ( i = 1; i <= line_num; i++ )
    {
      input.getline ( line, sizeof ( line ) );
      jlo = jhi + 1;
      jhi = i4_min ( jlo + r - 1, number );
 
      khi = 0;
      for ( j = jlo; j <= jhi; j++ )
      { 
        klo = khi + 1;
        khi = i4_min ( klo + w - 1, strlen ( line ) );
        s = s_substring ( line, klo, khi );
        rowind[j-1] = atoi ( s );
      }
    }
  }
  else
  {
    cerr << "\n";
    cerr << "HB_STRUCTURE_READ - Fatal error!\n";
    cerr << "  Illegal value of MXTYPE character 3.\n";
    exit ( 1 );
  }

  return;
}
//****************************************************************************80

void hb_structure_write ( ofstream &output, int ncol, char *mxtype, 
  int nnzero, int neltvl, char *ptrfmt, char *indfmt, int colptr[], 
  int rowind[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_STRUCTURE_WRITE writes the structure of an HB matrix.
//
//  Discussion:
//
//    If the user is creating an HB file, then the user should
//    already have opened the file, and written the header records.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 April 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ofstream &OUTPUT, the unit to which data is written.
//
//    Input, int NCOL, the number of columns.
//
//    Input, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Input, int NELTVL, the number of finite element matrix entries,
//    set to 0 in the case of assembled matrices.
//
//    Input, char *PTRFMT, the 16 character format for reading pointers.
//
//    Input, char *INDFMT, the 16 character format for reading indices.
//
//    Input, int COLPTR[NCOL+1], COLPTR[I-1] points to the location of
//    the first entry of column I in the sparse matrix structure.
//
//    If MXTYPE[2] == 'A':
//
//      Input, int ROWIND[NNZERO], the row index of each item.
//
//    If MXTYPE[2] == 'F':
//
//      Input, int ROWIND[NELTVL], the row index of each item.
//
{
  char code;
  int j;
  int m;
  int number;
  int r;
  int w;

  s_to_format ( ptrfmt, &r, &code, &w, &m );

  for ( j = 1; j <= ncol + 1; j++ )
  {
    output << setw(w) << colptr[j-1];
    if ( j % r == 0 )
    {
      output << "\n";
    }
  }
  if ( ( ncol + 1 ) % r != 0 )
  {
    output << "\n";
  }

  s_to_format ( indfmt, &r, &code, &w, &m );

  if ( mxtype[2] == 'A' )
  {
    for ( j = 1; j <= nnzero; j++ )
    {
      output << setw(w) << rowind[j-1];
      if ( j % r == 0 )
      {
        output << "\n";
      }
    }
    if ( nnzero % r != 0 )
    {
      output << "\n";
    }
  }
  else if ( mxtype[2] == 'E' )
  {
    number = colptr[ncol-1] - colptr[0];

    for ( j = 1; j <= number; j++ )
    {
      output << setw(w) << rowind[j-1];
      if ( j % r == 0 )
      {
        output << "\n";
      }
    }
    if ( number % r != 0 )
    {
      output << "\n";
    }
  }
  else
  {
    cerr << "\n";
    cerr << "HB_STRUCTURE_WRITE - Fatal error!\n";
    cerr << "  Illegal value of MXTYPE character 3.\n";
    exit ( 1 );
  }

  return;
}
//****************************************************************************80

int *hb_ua_colind ( int ncol, int colptr[], int nnzero )

//****************************************************************************80
//
//  Purpose:
//
//    HB_UA_COLUMN_INDEX creates a column index for an unsymmetric assembled matrix.
//
//  Discussion:
//
//    It is assumed that the input data corresponds to a Harwell-Boeing
//    matrix which is unsymmetric, and assembled.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    22 September 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, integer NCOL, the number of columns.
//
//    Input, integer COLPTR[NCOL+1], COLPTR(I) points to the location of
//    the first entry of column I in the sparse matrix structure.
//
//    Input, int NNZERO, the number of nonzeros.
//
//    Output, int HB_UA_COLIND[NNZERO], the column index of each matrix value.
//
{
  int *colind;
  int i;
  int j;

  colind = new int[nnzero];

  for ( i = 1; i <= ncol; i++ )
  {
    for ( j = colptr[i-1]; j <= colptr[i] - 1; j++ )
    {
      colind[j-1] = i;
    }
  }

  return colind;
}
//****************************************************************************80

void hb_values_print ( int ncol, int colptr[], char *mxtype, int nnzero, 
  int neltvl, double values[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_VALUES_PRINT prints the values of an HB matrix.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, int NCOL, the number of columns.
//
//    Input, int COLPTR[NCOL+1], COLPTR[I-1] points to the location of
//    the first entry of column I in the sparse matrix structure.
//
//    Input, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Input, int NELTVL, the number of finite element matrix entries,
//    set to 0 in the case of assembled matrices.
//
//    If MXTYPE[2] == 'A':
//
//      Input, double VALUES[NNZERO], the nonzero values of the matrix.
//
//    If MXTYPE[2] == 'E':
//
//      Input, double VALUES[NELTVL], the nonzero values of the matrix.
//
{
  int j;
  int k;
  int khi;
  int klo;

  if ( mxtype[2] == 'A' )
  {
    cout << "\n";
    cout << "Column Begin   End   ----------------------------------------\n";
    cout << "\n";

    for ( j = 1; j <= ncol; j++ )
    {
      if ( 5 < j && j < ncol )
      {
        continue;
      }

      if ( j == ncol && 6 < ncol )
      {
        cout << "Skipping intermediate columns...)\n";
      }

      if ( colptr[j]-1 < colptr[j-1] )
      {
        cout << setw(6) << j << "   EMPTY\n";
      }
      else
      {
        for ( klo = colptr[j-1]; klo <= colptr[j]-1; klo = klo + 5 )
        {
          khi = i4_min ( klo + 4, colptr[j]-1 );
          if ( klo == colptr[j-1] )
          {
            cout << setw(5) << j
                 << setw(5) << colptr[j-1]
                 << setw(5) << colptr[j]-1 << "   ";
          }
          for ( k = klo; k <= khi; k++ )
          {
            cout << setw(12) << values[k-1];
          }
          cout << "\n";
        }
      }
    }

    cout << "                     ----------------------------------------\n";
  }
  else if ( mxtype[2] == 'E' )
  {
    cout << "\n";
    cout << "Column Begin   End   ----------------------------------------\n";
    cout << "                     ----------------------------------------\n";

    cout << "\n";
    cout << "I haven't thought about how to print an\n";
    cout << "unassembled matrix yet!\n";
  }
  else
  {
    cerr << "\n";
    cerr << "HB_VALUES_PRINT - Fatal error!\n";
    cerr << "  Illegal value of MXTYPE character 3 = " << mxtype[2] << "\n";
    exit ( 1 );
  }

  return;
}
//****************************************************************************80

void hb_values_read ( ifstream &input, int valcrd, char *mxtype, int nnzero,
  int neltvl, char *valfmt, double values[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_VALUES_READ reads the values of an HB matrix.
//
//  Discussion:
//
//    The user should already have opened the file, and positioned it
//    to just after the header and structure records.
//
//    Values are contained in an HB file if the VALCRD parameter
//    is nonzero.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ifstream &INPUT, the unit from which data is read.
//
//    Input, int VALCRD, the number of input lines for numerical values.
//
//    Input, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Input, int NELTVL, the number of finite element matrix entries,
//    set to 0 in the case of assembled matrices.
//
//    Input, char *VALFMT, the 20 character format for reading values.
//
//    If MXTYPE[2] == 'A':
//
//      Output, double VALUES[NNZERO], the nonzero values of the matrix.
//
//    If MXTYPE[2] == 'E':
//
//      Output, double VALUES[NELTVL], the nonzero values of the matrix.
//
{
  char code;
  int i;
  int j;
  int jhi;
  int jlo;
  int khi;
  int klo;
  char line[255];
  int line_num;
  int m;
  int r;
  char *s;
  int w;

  s_to_format ( valfmt, &r, &code, &w, &m );

//
//  Read the matrix values.
//    case "A" = assembled;
//    case "E" = unassembled finite element matrices.
//
  if ( 0 < valcrd )
  {
    if ( mxtype[2] == 'A' )
    {
      line_num = 1 + ( nnzero - 1 ) / r;
    }
    else if ( mxtype[2] == 'E' )
    {
      line_num = 1 + ( neltvl - 1 ) / r;
    }
    else
    {
      cerr << "\n";
      cerr << "HB_VALUES_READ - Fatal error!\n";
      cerr << "  Illegal value of MXTYPE character 3.\n";
      exit ( 1 );
    }

    jhi = 0;
    for ( i = 1; i <= line_num; i++ )
    {
      input.getline ( line, sizeof ( line ) );
      jlo = jhi + 1;
      if ( mxtype[2] == 'A' )
      {
        jhi = i4_min ( jlo + r - 1, nnzero );
      }
      else
      {
        jhi = i4_min ( jlo + r - 1, neltvl );
      }
      khi = 0;
      for ( j = jlo; j <= jhi; j++ )
      {
        klo = khi + 1;
        khi = i4_min ( klo + w - 1, strlen ( line ) );
        s = s_substring ( line, klo, khi );
        values[j-1] = atof ( s );
      }
    }
  }

  return;
}
//****************************************************************************80

void hb_values_write ( ofstream &output, int valcrd, char *mxtype, 
  int nnzero, int neltvl, char *valfmt, double values[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_VALUES_WRITE writes the values of an HB matrix.
//
//  Discussion:
//
//    If the user is creating an HB file, then the user should already
//    have opened the file, and written the header and structure records.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, ofstream &OUTPUT, the unit to which data is written.
//
//    Input, int VALCRD, the number of input lines for numerical values.
//
//    Input, char *MXTYPE, the 3 character matrix type.
//    First character is R for Real, C for complex, P for pattern only.
//    Second character is S for symmetric, U for unsymmetric, H for
//      Hermitian, Z for skew symmetric, R for rectangular.
//    Third character is A for assembled and E for unassembled
//      finite element matrices.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  In the case of unassembled finite
//    element matrices, in which the right hand side vectors are also
//    stored as unassembled finite element vectors, this is the total
//    number of entries in a single unassembled right hand side vector.
//
//    Input, int NELTVL, the number of finite element matrix entries,
//    set to 0 in the case of assembled matrices.
//
//    Input, char *VALFMT, the 20 character format for reading values.
//
//    If MXTYPE[2] == 'A':
//
//      Input, double VALUES[NNZERO], the nonzero values of the matrix.
//
//    If MXTYPE[2] == 'E':
//
//      Input, double VALUES[NELTVL], the nonzero values of the matrix.
//
{
  char code;
  int i;
  int j;
  int jhi;
  int jlo;
  int line_num;
  int m;
  int r;
  int w;

  if ( 0 < valcrd )
  {
    s_to_format ( valfmt, &r, &code, &w, &m );

    if ( mxtype[2] == 'A' )
    {
      line_num = 1 + ( nnzero - 1 ) / r;
    }
    else if ( mxtype[2] == 'E' )
    {
      line_num = 1 + ( neltvl - 1 ) / r;
    }
    else
    {
      cerr << "\n";
      cerr << "HB_VALUES_WRITE - Fatal error!\n";
      cerr << "  Illegal value of MXTYPE character 3.\n";
      exit ( 1 );
    }

    jhi = 0;
    for ( i = 1; i <= line_num; i++ )
    {
      jlo = jhi + 1;
      if ( mxtype[2] == 'A' )
      {
        jhi = i4_min ( jlo + r - 1, nnzero );
      }
      else
      {
        jhi = i4_min ( jlo + r - 1, neltvl );
      }

      for ( j = jlo; j <= jhi; j++ )
      {
        output << setw(w) << values[j-1];
      }
      output << "\n";
    }
  }

  return;
}
//****************************************************************************80

double *hb_vecmat_a_mem ( int nrow, int ncol, int nnzero, int nrhs, 
  int colptr[], int rowind[], double values[], double exact[] )

//****************************************************************************80
//
//  Purpose:
//
//    HB_VECMAT_A_MEM multiplies a vector times an assembled Harwell Boeing matrix.
//
//  Discussion:
//
//    In this "A_MEM" version of the routine, the matrix is assumed to be in
//    "assembled" form, and all the data is assumed to be small enough
//    to reside completely in memory; the matrix and multiplicand vectors
//    are assumed to have been read into memory before this routine is called.
//
//    It is assumed that MXTYPE(3:3) = 'A', that is, that the matrix is
//    stored in the "assembled" format.
//
//    Also, the storage used for the vectors X and the products A*X
//    corresponds to RHSTYP(1:1) = 'F', that is, the "full" storage mode
//    for vectors.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    21 January 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Iain Duff, Roger Grimes, John Lewis,
//    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
//    October 1992.
//
//  Parameters:
//
//    Input, int NROW, the number of rows or variables.
//
//    Input, int NCOL, the number of columns or elements.
//
//    Input, int NNZERO.  In the case of assembled sparse matrices,
//    this is the number of nonzeroes.  
//
//    Input, int NRHS, the number of right hand sides.
//
//    Input, int COLPTR[NCOL+1], COLPTR(I) points to the location of
//    the first entry of column I in the sparse matrix structure.
//
//    Input, int ROWIND[NNZERO], the row index of each item.
//
//    Input, double VALUES[NNZERO], the nonzero values of the matrix.
//
//    Input, double EXACT[NROW*NRHS], contains NRHS dense vectors.
//
//    Output, double HB_VECMAT_A_MEM[NCOL*NRHS], the product vectors A'*X.
//
{
  int column;
  int k;
  double *rhsval;
  int rhs;
  int row;

  rhsval = new double[ncol*nrhs];
//
//  Zero out the result vectors.
//
  for ( rhs = 1; rhs <= nrhs; rhs++ )
  {
    for ( column = 1; column <= ncol; column++ )
    {
      rhsval[column-1+(rhs-1)*ncol] = 0.0E+00;
    }
  }
//
//  For each column J of the matrix,
//
  for ( column = 1; column <= ncol; column++ )
  {
//
//  For nonzero entry K
//
    for ( k = colptr[column-1]; k <= colptr[column]-1; k++ )
    {
      row = rowind[k-1];
//
//  For each right hand side vector:
//
//    B(J,1:NRHS) = B(J,1:NRHS) + X(I,1:NRHS) * A(I,J)
//
      for ( rhs = 1; rhs <= nrhs; rhs++ )
      {
        rhsval[column-1+(rhs-1)*ncol] = rhsval[column-1+(rhs-1)*ncol] 
          + values[k-1] * exact[row-1+(rhs-1)*nrow];
      }
    }
  }

  return rhsval;
}
//****************************************************************************80

int i4_max ( int i1, int i2 )

//****************************************************************************80
//
//  Purpose:
//
//    I4_MAX returns the maximum of two I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 October 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I1, I2, are two integers to be compared.
//
//    Output, int I4_MAX, the larger of I1 and I2.
//
//
{
  if ( i2 < i1 )
  {
    return i1;
  }
  else
  {
    return i2;
  }

}
//****************************************************************************80

int i4_min ( int i1, int i2 )

//****************************************************************************80
//
//  Purpose:
//
//    I4_MIN returns the smaller of two I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 October 1998
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int I1, I2, two integers to be compared.
//
//    Output, int I4_MIN, the smaller of I1 and I2.
//
//
{
  if ( i1 < i2 )
  {
    return i1;
  }
  else
  {
    return i2;
  }

}
//****************************************************************************80

void i4vec_print ( int n, int a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    I4VEC_PRINT prints an I4VEC.
//
//  Discussion:
//
//    An I4VEC is a vector of I4's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    14 November 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of components of the vector.
//
//    Input, int A[N], the vector to be printed.
//
//    Input, string TITLE, a title.
//
{
  int i;

  cout << "\n";
  cout << title << "\n";
  cout << "\n";
  for ( i = 0; i < n; i++ ) 
  {
    cout << "  " << setw(8) << i 
         << ": " << setw(8) << a[i]  << "\n";
  }
  return;
}
//****************************************************************************80

void i4vec_print_part ( int n, int a[], int max_print, string title )

//****************************************************************************80
//
//  Purpose:
//
//    I4VEC_PRINT_PART prints "part" of an I4VEC.
//
//  Discussion:
//
//    The user specifies MAX_PRINT, the maximum number of lines to print.
//
//    If N, the size of the vector, is no more than MAX_PRINT, then
//    the entire vector is printed, one entry per line.
//
//    Otherwise, if possible, the first MAX_PRINT-2 entries are printed,
//    followed by a line of periods suggesting an omission,
//    and the last entry.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    27 February 2010
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries of the vector.
//
//    Input, int A[N], the vector to be printed.
//
//    Input, int MAX_PRINT, the maximum number of lines
//    to print.
//
//    Input, string TITLE, a title.
//
{
  int i;

  if ( max_print <= 0 )
  {
    return;
  }

  if ( n <= 0 )
  {
    return;
  }

  cout << "\n";
  cout << title << "\n";
  cout << "\n";

  if ( n <= max_print )
  {
    for ( i = 0; i < n; i++ )
    {
      cout << "  " << setw(8) << i
           << "  " << setw(8) << a[i] << "\n";
    }
  }
  else if ( 3 <= max_print )
  {
    for ( i = 0; i < max_print - 2; i++ )
    {
      cout << "  " << setw(8) << i
           << ": " << setw(8) << a[i] << "\n";
    }
    cout << "  ........  ........\n";
    i = n - 1;
    cout << "  " << setw(8) << i
         << ": " << setw(8) << a[i] << "\n";
  }
  else
  {
    for ( i= 0; i < max_print - 1; i++ )
    {
      cout << "  " << setw(8) << i
           << ": " << setw(8) << a[i] << "\n";
    }
    i = max_print - 1;
    cout << "  " << setw(8) << i
         << ": " << setw(8) << a[i]
         << "  " << "...more entries...\n";
  }

  return;
}
//****************************************************************************80

void r8mat_print ( int m, int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_PRINT prints an R8MAT.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//    Entry A(I,J) is stored as A[I+J*M]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    10 September 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows in A.
//
//    Input, int N, the number of columns in A.
//
//    Input, double A[M*N], the M by N matrix.
//
//    Input, string TITLE, a title.
//
{
  r8mat_print_some ( m, n, a, 1, 1, m, n, title );

  return;
}
//****************************************************************************80

void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, string title )

//****************************************************************************80
//
//  Purpose:
//
//    R8MAT_PRINT_SOME prints some of an R8MAT.
//
//  Discussion:
//
//    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
//    in column-major order.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    26 June 2013
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int M, the number of rows of the matrix.
//    M must be positive.
//
//    Input, int N, the number of columns of the matrix.
//    N must be positive.
//
//    Input, double A[M*N], the matrix.
//
//    Input, int ILO, JLO, IHI, JHI, designate the first row and
//    column, and the last row and column to be printed.
//
//    Input, string TITLE, a title.
//
{
# define INCX 5

  int i;
  int i2hi;
  int i2lo;
  int j;
  int j2hi;
  int j2lo;

  cout << "\n";
  cout << title << "\n";

  if ( m <= 0 || n <= 0 )
  {
    cout << "\n";
    cout << "  (None)\n";
    return;
  }
//
//  Print the columns of the matrix, in strips of 5.
//
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    if ( n < j2hi )
    {
      j2hi = n;
    }
    if ( jhi < j2hi )
    {
      j2hi = jhi;
    }
    cout << "\n";
//
//  For each column J in the current range...
//
//  Write the header.
//
    cout << "  Col:    ";
    for ( j = j2lo; j <= j2hi; j++ )
    {
      cout << setw(7) << j - 1 << "       ";
    }
    cout << "\n";
    cout << "  Row\n";
    cout << "\n";
//
//  Determine the range of the rows in this strip.
//
    if ( 1 < ilo )
    {
      i2lo = ilo;
    }
    else
    {
      i2lo = 1;
    }
    if ( ihi < m )
    {
      i2hi = ihi;
    }
    else
    {
      i2hi = m;
    }

    for ( i = i2lo; i <= i2hi; i++ )
    {
//
//  Print out (up to) 5 entries in row I, that lie in the current strip.
//
      cout << setw(5) << i - 1 << ": ";
      for ( j = j2lo; j <= j2hi; j++ )
      {
        cout << setw(12) << a[i-1+(j-1)*m] << "  ";
      }
      cout << "\n";
    }
  }

  return;
# undef INCX
}
//****************************************************************************80

void r8vec_print ( int n, double a[], string title )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_PRINT prints an R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    16 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of components of the vector.
//
//    Input, double A[N], the vector to be printed.
//
//    Input, string TITLE, a title.
//
{
  int i;

  cout << "\n";
  cout << title << "\n";
  cout << "\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  " << setw(8)  << i
         << ": " << setw(14) << a[i]  << "\n";
  }

  return;
}
//****************************************************************************80

void r8vec_print_part ( int n, double a[], int max_print, string title )

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_PRINT_PART prints "part" of an R8VEC.
//
//  Discussion:
//
//    The user specifies MAX_PRINT, the maximum number of lines to print.
//
//    If N, the size of the vector, is no more than MAX_PRINT, then
//    the entire vector is printed, one entry per line.
//
//    Otherwise, if possible, the first MAX_PRINT-2 entries are printed,
//    followed by a line of periods suggesting an omission,
//    and the last entry.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    27 February 2010
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries of the vector.
//
//    Input, double A[N], the vector to be printed.
//
//    Input, int MAX_PRINT, the maximum number of lines
//    to print.
//
//    Input, string TITLE, a title.
//
{
  int i;

  if ( max_print <= 0 )
  {
    return;
  }

  if ( n <= 0 )
  {
    return;
  }

  cout << "\n";
  cout << title << "\n";
  cout << "\n";

  if ( n <= max_print )
  {
    for ( i = 0; i < n; i++ )
    {
      cout << "  " << setw(8) << i
           << "  " << setw(14) << a[i] << "\n";
    }
  }
  else if ( 3 <= max_print )
  {
    for ( i = 0; i < max_print - 2; i++ )
    {
      cout << "  " << setw(8) << i
           << ": " << setw(14) << a[i] << "\n";
    }
    cout << "  ........  ..............\n";
    i = n - 1;
    cout << "  " << setw(8) << i
         << ": " << setw(14) << a[i] << "\n";
  }
  else
  {
    for ( i = 0; i < max_print - 1; i++ )
    {
      cout << "  " << setw(8) << i
           << ": " << setw(14) << a[i] << "\n";
    }
    i = max_print - 1;
    cout << "  " << setw(8) << i
         << ": " << setw(14) << a[i]
         << "  " << "...more entries...\n";
  }

  return;
}
//****************************************************************************80

int s_len_trim ( char *s )

//****************************************************************************80
//
//  Purpose:
//
//    S_LEN_TRIM returns the length of a string to the last nonblank.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    26 April 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *S, a pointer to a string.
//
//    Output, int S_LEN_TRIM, the length of the string to the last nonblank.
//    If S_LEN_TRIM is 0, then the string is entirely blank.
//
{
  int n;
  char* t;

  n = strlen ( s );
  t = s + strlen ( s ) - 1;

  while ( 0 < n ) 
  {
    if ( *t != ' ' )
    {
      return n;
    }
    t--;
    n--;
  }

  return n;
}
//****************************************************************************80

char *s_substring ( char *s, int a, int b )

//****************************************************************************80
//
//  Purpose:
//
//    S_SUBSTRING returns a substring of a given string.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    10 April 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *S, a pointer to a string.
//
//    Input, int A, B, the indices of the first and last character of S to copy.
//    These are 1-based indices!  B should be 
//
//    Output, char *S_SUBSTRING, a pointer to the substring.
//
{
  int i;
  int j;
  char *t;

  t = new char[b+2-a];

  j = 0;
  for ( i = a; i <= b; i++ )
  {
    t[j] = s[i-1];
    j = j + 1;
  }
  t[j] = '\0';

  return t;
}
//****************************************************************************80

void s_to_format ( char *s, int *r, char *code, int *w, int *m )

//****************************************************************************80
//
//  Purpose:
//
//    S_TO_FORMAT reads a FORTRAN format from a string.
//
//  Discussion:
//
//    This routine will read as many characters as possible until it reaches
//    the end of the string, or encounters a character which cannot be
//    part of the format.  This routine is limited in its ability to
//    recognize FORTRAN formats.  In particular, we are only expecting
//    a single format specification, and cannot handle extra features
//    such as 'ES' and 'EN' codes, '5X' spacing, and so on.
//
//    Legal input is:
//
//       0 nothing
//       1 blanks
//       2 optional '('
//       3 blanks
//       4 optional repeat factor R
//       5 blanks
//       6 CODE ( 'A', 'B', 'E', 'F', 'G', 'I', 'L', 'O', 'Z', '*' )
//       7 blanks
//       8 width W
//       9 optional decimal point
//      10 optional mantissa M
//      11 blanks
//      12 optional ')'
//      13 blanks
//
//  Example:
//
//    S                 R   CODE   W    M
//
//    'I12              1   I      12   0
//    'E8.0'            1   E       8   0
//    'F10.5'           1   F      10   5
//    '2G14.6'          2   G      14   6
//    '*'               1   *      -1  -1
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    10 April 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *S, the string containing the
//    data to be read.  Reading will begin at position 1 and
//    terminate at the end of the string, or when no more
//    characters can be read.
//
//    Output, int *R, the repetition factor, which defaults to 1.
//
//    Output, char *CODE, the format code.
//
//    Output, int *W, the field width.
//
//    Output, int *M, the mantissa width.
//
{
  char c;
  int d;
  bool debug = true;
  int LEFT = 1;
  int paren_sum;
  int pos;
  int RIGHT = -1;
  int s_length;
  int state;

  state = 0;
  paren_sum = 0;
  pos = 0;
  s_length = s_len_trim ( s );

  *r = 0;
  *w = 0;
  *code = '?';
  *m = 0;

  while ( pos < s_length )
  {
    c = s[pos];
    pos = pos + 1;
//
//  BLANK character:
//
    if ( c == ' ' )
    {
      if ( state == 4 )
      {
        state = 5;
      }
      else if ( state == 6 )
      {
        state = 7;
      }
      else if ( state == 10 )
      {
        state = 11;
      }
      else if ( state == 12 )
      {
        state = 13;
      }
    }
//
//  LEFT PAREN
//
    else if ( c == '(' )
    {
      if ( state < 2 )
      {
        paren_sum = paren_sum + LEFT;
      }
      else
      {
        if ( debug )
        {
          cerr << "\n";
          cerr << "S_TO_FORMAT - Fatal error!\n";
          cerr << "  Current state = " << state << "\n";
          cerr << "  Input character = '" << c << "'.\n";
        }
        state = -1;
        break;
      }
    }
//
//  DIGIT (R, F, or W)
//
    else if ( ch_is_digit ( c ) )
    {
      if ( state <= 3 )
      {
        state = 4;
        *r = ch_to_digit ( c );
      }
      else if ( state == 4 )
      {
        d = ch_to_digit ( c );
        *r = 10 * (*r) + d;
      }
      else if ( state == 6 || state == 7 )
      {
        if ( *code == '*' )
        {
          if ( debug )
          {
            cerr << "\n";
            cerr << "S_TO_FORMAT - Fatal error!\n";
            cerr << "  Current state = " << state << "\n";
            cerr << "  Current code = '" << *code << "'.\n";
            cerr << "  Input character = '" << c << "'.\n";
          }
          state = -1;
          break;
        }
        state = 8;
        *w = ch_to_digit ( c );
      }
      else if ( state == 8 )
      {
        d = ch_to_digit ( c );
        *w = 10 * (*w) + d;
      }
      else if ( state == 9 )
      {
        state = 10;
        *m = ch_to_digit ( c );
      }
      else if ( state == 10 )
      {
        d = ch_to_digit ( c );
        *m = 10 * (*m) + d;
      }
      else
      {
        if ( debug )
        {
          cerr << "\n";
          cerr << "S_TO_FORMAT - Fatal error!\n";
          cerr << "  Current state = " << state << "\n";
          cerr << "  Input character = '" << c << "'.\n";
        }
        state = -1;
        break;
      }
    }
//
//  DECIMAL POINT
//
    else if ( c == '.' )
    {
      if ( state == 8 )
      {
        state = 9;
      }
      else
      {
        if ( debug )
        {
          cerr << "\n";
          cerr << "S_TO_FORMAT - Fatal error!\n";
          cerr << "  Current state = " << state << "\n";
          cerr << "  Input character = '" << c << "'.\n";
        }
        state = -1;
        break;
      }
    }
//
//  RIGHT PAREN
//
    else if ( c == ')' )
    {
      paren_sum = paren_sum + RIGHT;

      if ( paren_sum != 0 )
      {
        if ( debug )
        {
          cerr << "\n";
          cerr << "S_TO_FORMAT - Fatal error!\n";
          cerr << "  Current paren sum = " << paren_sum << "\n";
          cerr << "  Input character = '" << c << "'.\n";
        }
        state = -1;
        break;
      }

      if ( state == 6 && *code == '*' )
      {
        state = 12;
      }
      else if ( 6 <= state )
      {
        state = 12;
      }
      else
      {
        if ( debug )
        {
          cerr << "\n";
          cerr << "S_TO_FORMAT - Fatal error!\n";
          cerr << "  Current state = " << state << "\n";
          cerr << "  Input character = '" << c << "'.\n";
        }
        state = -1;
        break;
      }
    }
//
//  Code
//
    else if ( ch_is_format_code ( c ) )
    {
      if ( state < 6 )
      {
        state = 6;
        *code = c;
      }
      else
      {
        if ( debug )
        {
          cerr << "\n";
          cerr << "S_TO_FORMAT - Fatal error!\n";
          cerr << "  Current state = " << state << "\n";
          cerr << "  Input character = '" << c << "'.\n";
        }
        state = -1;
        break;
      }
    }
//
//  Unexpected character
//
    else
    {
      if ( debug )
      {
        cerr << "\n";
        cerr << "S_TO_FORMAT - Fatal error!\n";
        cerr << "  Current state = " << state << "\n";
        cerr << "  Input character = '" << c << "'.\n";
      }
      state = -1;
      break;
    }
  }

  if ( paren_sum != 0 )
  {
    cerr << "\n";
    cerr << "S_TO_FORMAT - Fatal error!\n";
    cerr << "  Parentheses mismatch.\n";
    exit ( 1 );
  }

  if ( state < 0 )
  {
    cerr << "\n";
    cerr << "S_TO_FORMAT - Fatal error!\n";
    cerr << "  Parsing error.\n";
    exit ( 1 );
  }

  if ( *r == 0 )
  {
    *r = 1;
  }

  return;
}
//****************************************************************************80

void s_trim ( char *s )

//****************************************************************************80
//
//  Purpose:
//
//    S_TRIM promotes the final null forward through trailing blanks.
//
//  Discussion:
//
//    What we're trying to say is that we reposition the null character
//    so that trailing blanks are no longer visible.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    10 April 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input/output, char *S, the string to be trimmed.
//
{
  char c;
  int n;
  char *t;

  n = strlen ( s );
  t = s + strlen ( s ) - 1;

  while ( 0 < n ) 
  {
    if ( *t != ' ' )
    {
      return;
    }
    c      = *t;
    *t     = *(t+1);
    *(t+1) = c;
    t--;
    n--;
  }

  return;
}
//**********************************************************************

void timestamp ( void )

//**********************************************************************
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    24 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}