Adding upstream version 16.2.11+ds.upstream/16.2.11+dsupstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 1610 insertions, 0 deletions

diff --git a/src/boost/libs/graph/example/iohb.c b/src/boost/libs/graph/example/iohb.c
new file mode 100644
index 000000000..8a216358d
--- /dev/null
+++ b/src/boost/libs/graph/example/iohb.c
@@ -0,0 +1,1610 @@
+//  (C) Copyright Jeremy Siek 2004
+//  Distributed under the Boost Software License, Version 1.0. (See
+//  accompanying file LICENSE_1_0.txt or copy at
+//  http://www.boost.org/LICENSE_1_0.txt)
+
+/*
+Fri Aug 15 16:29:47 EDT 1997
+
+                       Harwell-Boeing File I/O in C
+                                V. 1.0
+
+           National Institute of Standards and Technology, MD.
+                             K.A. Remington
+
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+                                NOTICE
+
+ Permission to use, copy, modify, and distribute this software and
+ its documentation for any purpose and without fee is hereby granted
+ provided that the above copyright notice appear in all copies and
+ that both the copyright notice and this permission notice appear in
+ supporting documentation.
+
+ Neither the Author nor the Institution (National Institute of Standards
+ and Technology) make any representations about the suitability of this
+ software for any purpose. This software is provided "as is" without
+ expressed or implied warranty.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+                         ---------------------
+                         INTERFACE DESCRIPTION
+                         ---------------------
+  ---------------
+  QUERY FUNCTIONS
+  ---------------
+
+  FUNCTION:
+
+  int readHB_info(const char *filename, int *M, int *N, int *nz,
+  char **Type, int *Nrhs)
+
+  DESCRIPTION:
+
+  The readHB_info function opens and reads the header information from
+  the specified Harwell-Boeing file, and reports back the number of rows
+  and columns in the stored matrix (M and N), the number of nonzeros in
+  the matrix (nz), the 3-character matrix type(Type), and the number of
+  right-hand-sides stored along with the matrix (Nrhs).  This function
+  is designed to retrieve basic size information which can be used to
+  allocate arrays.
+
+  FUNCTION:
+
+  int  readHB_header(FILE* in_file, char* Title, char* Key, char* Type,
+                    int* Nrow, int* Ncol, int* Nnzero, int* Nrhs,
+                    char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
+                    int* Ptrcrd, int* Indcrd, int* Valcrd, int* Rhscrd,
+                    char *Rhstype)
+
+  DESCRIPTION:
+
+  More detailed than the readHB_info function, readHB_header() reads from
+  the specified Harwell-Boeing file all of the header information.
+
+
+  ------------------------------
+  DOUBLE PRECISION I/O FUNCTIONS
+  ------------------------------
+  FUNCTION:
+
+  int readHB_newmat_double(const char *filename, int *M, int *N, *int nz,
+  int **colptr, int **rowind,  double**val)
+
+  int readHB_mat_double(const char *filename, int *colptr, int *rowind,
+  double*val)
+
+
+  DESCRIPTION:
+
+  This function opens and reads the specified file, interpreting its
+  contents as a sparse matrix stored in the Harwell/Boeing standard
+  format.  (See readHB_aux_double to read auxillary vectors.)
+        -- Values are interpreted as double precision numbers. --
+
+  The "mat" function uses _pre-allocated_ vectors to hold the index and
+  nonzero value information.
+
+  The "newmat" function allocates vectors to hold the index and nonzero
+  value information, and returns pointers to these vectors along with
+  matrix dimension and number of nonzeros.
+
+  FUNCTION:
+
+  int readHB_aux_double(const char* filename, const char AuxType, double b[])
+
+  int readHB_newaux_double(const char* filename, const char AuxType, double** b)
+
+  DESCRIPTION:
+
+  This function opens and reads from the specified file auxillary vector(s).
+  The char argument Auxtype determines which type of auxillary vector(s)
+  will be read (if present in the file).
+
+                  AuxType = 'F'   right-hand-side
+                  AuxType = 'G'   initial estimate (Guess)
+                  AuxType = 'X'   eXact solution
+
+  If Nrhs > 1, all of the Nrhs vectors of the given type are read and
+  stored in column-major order in the vector b.
+
+  The "newaux" function allocates a vector to hold the values retrieved.
+  The "mat" function uses a _pre-allocated_ vector to hold the values.
+
+  FUNCTION:
+
+  int writeHB_mat_double(const char* filename, int M, int N,
+                        int nz, const int colptr[], const int rowind[],
+                        const double val[], int Nrhs, const double rhs[],
+                        const double guess[], const double exact[],
+                        const char* Title, const char* Key, const char* Type,
+                        char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
+                        const char* Rhstype)
+
+  DESCRIPTION:
+
+  The writeHB_mat_double function opens the named file and writes the specified
+  matrix and optional auxillary vector(s) to that file in Harwell-Boeing
+  format.  The format arguments (Ptrfmt,Indfmt,Valfmt, and Rhsfmt) are
+  character strings specifying "Fortran-style" output formats -- as they
+  would appear in a Harwell-Boeing file.  They are used to produce output
+  which is as close as possible to what would be produced by Fortran code,
+  but note that "D" and "P" edit descriptors are not supported.
+  If NULL, the following defaults will be used:
+                    Ptrfmt = Indfmt = "(8I10)"
+                    Valfmt = Rhsfmt = "(4E20.13)"
+
+  -----------------------
+  CHARACTER I/O FUNCTIONS
+  -----------------------
+  FUNCTION:
+
+  int readHB_mat_char(const char* filename, int colptr[], int rowind[],
+                                           char val[], char* Valfmt)
+  int readHB_newmat_char(const char* filename, int* M, int* N, int* nonzeros,
+                          int** colptr, int** rowind, char** val, char** Valfmt)
+
+  DESCRIPTION:
+
+  This function opens and reads the specified file, interpreting its
+  contents as a sparse matrix stored in the Harwell/Boeing standard
+  format.  (See readHB_aux_char to read auxillary vectors.)
+              -- Values are interpreted as char strings.     --
+  (Used to translate exact values from the file into a new storage format.)
+
+  The "mat" function uses _pre-allocated_ arrays to hold the index and
+  nonzero value information.
+
+  The "newmat" function allocates char arrays to hold the index
+  and nonzero value information, and returns pointers to these arrays
+  along with matrix dimension and number of nonzeros.
+
+  FUNCTION:
+
+  int readHB_aux_char(const char* filename, const char AuxType, char b[])
+  int readHB_newaux_char(const char* filename, const char AuxType, char** b,
+                         char** Rhsfmt)
+
+  DESCRIPTION:
+
+  This function opens and reads from the specified file auxillary vector(s).
+  The char argument Auxtype determines which type of auxillary vector(s)
+  will be read (if present in the file).
+
+                  AuxType = 'F'   right-hand-side
+                  AuxType = 'G'   initial estimate (Guess)
+                  AuxType = 'X'   eXact solution
+
+  If Nrhs > 1, all of the Nrhs vectors of the given type are read and
+  stored in column-major order in the vector b.
+
+  The "newaux" function allocates a character array to hold the values
+                retrieved.
+  The "mat" function uses a _pre-allocated_ array to hold the values.
+
+  FUNCTION:
+
+  int writeHB_mat_char(const char* filename, int M, int N,
+                        int nz, const int colptr[], const int rowind[],
+                        const char val[], int Nrhs, const char rhs[],
+                        const char guess[], const char exact[],
+                        const char* Title, const char* Key, const char* Type,
+                        char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
+                        const char* Rhstype)
+
+  DESCRIPTION:
+
+  The writeHB_mat_char function opens the named file and writes the specified
+  matrix and optional auxillary vector(s) to that file in Harwell-Boeing
+  format.  The format arguments (Ptrfmt,Indfmt,Valfmt, and Rhsfmt) are
+  character strings specifying "Fortran-style" output formats -- as they
+  would appear in a Harwell-Boeing file.  Valfmt and Rhsfmt must accurately
+  represent the character representation of the values stored in val[]
+  and rhs[].
+
+  If NULL, the following defaults will be used for the integer vectors:
+                    Ptrfmt = Indfmt = "(8I10)"
+                    Valfmt = Rhsfmt = "(4E20.13)"
+
+
+*/
+
+/*---------------------------------------------------------------------*/
+/* If zero-based indexing is desired, _SP_base should be set to 0      */
+/* This will cause indices read from H-B files to be decremented by 1  */
+/*             and indices written to H-B files to be incremented by 1 */
+/*            <<<  Standard usage is _SP_base = 1  >>>                 */
+#ifndef _SP_base
+#define _SP_base 1
+#endif
+/*---------------------------------------------------------------------*/
+
+#include "iohb.h"
+#include<stdio.h>
+#include<stdlib.h>
+#include<string.h>
+#include<math.h>
+
+char* substr(const char* S, const int pos, const int len);
+void upcase(char* S);
+void IOHBTerminate(const char* message);
+
+int readHB_info(const char* filename, int* M, int* N, int* nz, char** Type,
+                                                      int* Nrhs)
+{
+/****************************************************************************/
+/*  The readHB_info function opens and reads the header information from    */
+/*  the specified Harwell-Boeing file, and reports back the number of rows  */
+/*  and columns in the stored matrix (M and N), the number of nonzeros in   */
+/*  the matrix (nz), and the number of right-hand-sides stored along with   */
+/*  the matrix (Nrhs).                                                      */
+/*                                                                          */
+/*  For a description of the Harwell Boeing standard, see:                  */
+/*            Duff, et al.,  ACM TOMS Vol.15, No.1, March 1989              */
+/*                                                                          */
+/*    ----------                                                            */
+/*    **CAVEAT**                                                            */
+/*    ----------                                                            */
+/*  **  If the input file does not adhere to the H/B format, the  **        */
+/*  **             results will be unpredictable.                 **        */
+/*                                                                          */
+/****************************************************************************/
+    FILE *in_file;
+    int Ptrcrd, Indcrd, Valcrd, Rhscrd;
+    int Nrow, Ncol, Nnzero;
+    char* mat_type;
+    char Title[73], Key[9], Rhstype[4];
+    char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
+
+    mat_type = *Type;
+    if ( mat_type == NULL ) IOHBTerminate("Insufficient memory for mat_typen");
+
+    if ( (in_file = fopen( filename, "r")) == NULL ) {
+       fprintf(stderr,"Error: Cannot open file: %s\n",filename);
+       return 0;
+    }
+
+    readHB_header(in_file, Title, Key, mat_type, &Nrow, &Ncol, &Nnzero, Nrhs,
+                  Ptrfmt, Indfmt, Valfmt, Rhsfmt,
+                  &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
+    fclose(in_file);
+    *Type = mat_type;
+    *(*Type+3) = (char) NULL;
+    *M    = Nrow;
+    *N    = Ncol;
+    *nz   = Nnzero;
+    if (Rhscrd == 0) {*Nrhs = 0;}
+
+/*  In verbose mode, print some of the header information:   */
+/*
+    if (verbose == 1)
+    {
+        printf("Reading from Harwell-Boeing file %s (verbose on)...\n",filename);
+        printf("  Title: %s\n",Title);
+        printf("  Key:   %s\n",Key);
+        printf("  The stored matrix is %i by %i with %i nonzeros.\n",
+                *M, *N, *nz );
+        printf("  %i right-hand--side(s) stored.\n",*Nrhs);
+    }
+*/
+
+    return 1;
+
+}
+
+
+
+int readHB_header(FILE* in_file, char* Title, char* Key, char* Type,
+                    int* Nrow, int* Ncol, int* Nnzero, int* Nrhs,
+                    char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
+                    int* Ptrcrd, int* Indcrd, int* Valcrd, int* Rhscrd,
+                    char *Rhstype)
+{
+/*************************************************************************/
+/*  Read header information from the named H/B file...                   */
+/*************************************************************************/
+    int Totcrd,Neltvl,Nrhsix;
+    char line[BUFSIZ];
+
+/*  First line:   */
+    fgets(line, BUFSIZ, in_file);
+    if ( sscanf(line,"%*s") < 0 )
+        IOHBTerminate("iohb.c: Null (or blank) first line of HB file.\n");
+    (void) sscanf(line, "%72c%8[^\n]", Title, Key);
+    *(Key+8) = (char) NULL;
+    *(Title+72) = (char) NULL;
+
+/*  Second line:  */
+    fgets(line, BUFSIZ, in_file);
+    if ( sscanf(line,"%*s") < 0 )
+        IOHBTerminate("iohb.c: Null (or blank) second line of HB file.\n");
+    if ( sscanf(line,"%i",&Totcrd) != 1) Totcrd = 0;
+    if ( sscanf(line,"%*i%i",Ptrcrd) != 1) *Ptrcrd = 0;
+    if ( sscanf(line,"%*i%*i%i",Indcrd) != 1) *Indcrd = 0;
+    if ( sscanf(line,"%*i%*i%*i%i",Valcrd) != 1) *Valcrd = 0;
+    if ( sscanf(line,"%*i%*i%*i%*i%i",Rhscrd) != 1) *Rhscrd = 0;
+
+/*  Third line:   */
+    fgets(line, BUFSIZ, in_file);
+    if ( sscanf(line,"%*s") < 0 )
+        IOHBTerminate("iohb.c: Null (or blank) third line of HB file.\n");
+    if ( sscanf(line, "%3c", Type) != 1)
+        IOHBTerminate("iohb.c: Invalid Type info, line 3 of Harwell-Boeing file.\n");
+    upcase(Type);
+    if ( sscanf(line,"%*3c%i",Nrow) != 1) *Nrow = 0 ;
+    if ( sscanf(line,"%*3c%*i%i",Ncol) != 1) *Ncol = 0 ;
+    if ( sscanf(line,"%*3c%*i%*i%i",Nnzero) != 1) *Nnzero = 0 ;
+    if ( sscanf(line,"%*3c%*i%*i%*i%i",&Neltvl) != 1) Neltvl = 0 ;
+
+/*  Fourth line:  */
+    fgets(line, BUFSIZ, in_file);
+    if ( sscanf(line,"%*s") < 0 )
+        IOHBTerminate("iohb.c: Null (or blank) fourth line of HB file.\n");
+    if ( sscanf(line, "%16c",Ptrfmt) != 1)
+        IOHBTerminate("iohb.c: Invalid format info, line 4 of Harwell-Boeing file.\n");
+    if ( sscanf(line, "%*16c%16c",Indfmt) != 1)
+        IOHBTerminate("iohb.c: Invalid format info, line 4 of Harwell-Boeing file.\n");
+    if ( sscanf(line, "%*16c%*16c%20c",Valfmt) != 1)
+        IOHBTerminate("iohb.c: Invalid format info, line 4 of Harwell-Boeing file.\n");
+    sscanf(line, "%*16c%*16c%*20c%20c",Rhsfmt);
+    *(Ptrfmt+16) = (char) NULL;
+    *(Indfmt+16) = (char) NULL;
+    *(Valfmt+20) = (char) NULL;
+    *(Rhsfmt+20) = (char) NULL;
+
+/*  (Optional) Fifth line: */
+    if (*Rhscrd != 0 )
+    {
+       fgets(line, BUFSIZ, in_file);
+       if ( sscanf(line,"%*s") < 0 )
+           IOHBTerminate("iohb.c: Null (or blank) fifth line of HB file.\n");
+       if ( sscanf(line, "%3c", Rhstype) != 1)
+         IOHBTerminate("iohb.c: Invalid RHS type information, line 5 of Harwell-Boeing file.\n");
+       if ( sscanf(line, "%*3c%i", Nrhs) != 1) *Nrhs = 0;
+       if ( sscanf(line, "%*3c%*i%i", &Nrhsix) != 1) Nrhsix = 0;
+    }
+    return 1;
+}
+
+
+int readHB_mat_double(const char* filename, int colptr[], int rowind[],
+                                                                 double val[])
+{
+/****************************************************************************/
+/*  This function opens and reads the specified file, interpreting its      */
+/*  contents as a sparse matrix stored in the Harwell/Boeing standard       */
+/*  format and creating compressed column storage scheme vectors to hold    */
+/*  the index and nonzero value information.                                */
+/*                                                                          */
+/*    ----------                                                            */
+/*    **CAVEAT**                                                            */
+/*    ----------                                                            */
+/*  Parsing real formats from Fortran is tricky, and this file reader       */
+/*  does not claim to be foolproof.   It has been tested for cases when     */
+/*  the real values are printed consistently and evenly spaced on each      */
+/*  line, with Fixed (F), and Exponential (E or D) formats.                 */
+/*                                                                          */
+/*  **  If the input file does not adhere to the H/B format, the  **        */
+/*  **             results will be unpredictable.                 **        */
+/*                                                                          */
+/****************************************************************************/
+    FILE *in_file;
+    int i,j,ind,col,offset,count,last,Nrhs;
+    int Ptrcrd, Indcrd, Valcrd, Rhscrd;
+    int Nrow, Ncol, Nnzero, Nentries;
+    int Ptrperline, Ptrwidth, Indperline, Indwidth;
+    int Valperline, Valwidth, Valprec;
+    int Valflag;           /* Indicates 'E','D', or 'F' float format */
+    char* ThisElement;
+    char Title[73], Key[8], Type[4], Rhstype[4];
+    char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
+    char line[BUFSIZ];
+
+    if ( (in_file = fopen( filename, "r")) == NULL ) {
+       fprintf(stderr,"Error: Cannot open file: %s\n",filename);
+       return 0;
+    }
+
+    readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
+                  Ptrfmt, Indfmt, Valfmt, Rhsfmt,
+                  &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
+
+/*  Parse the array input formats from Line 3 of HB file  */
+    ParseIfmt(Ptrfmt,&Ptrperline,&Ptrwidth);
+    ParseIfmt(Indfmt,&Indperline,&Indwidth);
+    if ( Type[0] != 'P' ) {          /* Skip if pattern only  */
+    ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
+    }
+
+/*  Read column pointer array:   */
+
+    offset = 1-_SP_base;  /* if base 0 storage is declared (via macro definition), */
+                          /* then storage entries are offset by 1                  */
+
+    ThisElement = (char *) malloc(Ptrwidth+1);
+    if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
+    *(ThisElement+Ptrwidth) = (char) NULL;
+    count=0;
+    for (i=0;i<Ptrcrd;i++)
+    {
+       fgets(line, BUFSIZ, in_file);
+       if ( sscanf(line,"%*s") < 0 )
+         IOHBTerminate("iohb.c: Null (or blank) line in pointer data region of HB file.\n");
+       col =  0;
+       for (ind = 0;ind<Ptrperline;ind++)
+       {
+          if (count > Ncol) break;
+          strncpy(ThisElement,line+col,Ptrwidth);
+  /* ThisElement = substr(line,col,Ptrwidth); */
+          colptr[count] = atoi(ThisElement)-offset;
+          count++; col += Ptrwidth;
+       }
+    }
+    free(ThisElement);
+
+/*  Read row index array:  */
+
+    ThisElement = (char *) malloc(Indwidth+1);
+    if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
+    *(ThisElement+Indwidth) = (char) NULL;
+    count = 0;
+    for (i=0;i<Indcrd;i++)
+    {
+       fgets(line, BUFSIZ, in_file);
+       if ( sscanf(line,"%*s") < 0 )
+         IOHBTerminate("iohb.c: Null (or blank) line in index data region of HB file.\n");
+       col =  0;
+       for (ind = 0;ind<Indperline;ind++)
+       {
+          if (count == Nnzero) break;
+          strncpy(ThisElement,line+col,Indwidth);
+/*        ThisElement = substr(line,col,Indwidth); */
+          rowind[count] = atoi(ThisElement)-offset;
+          count++; col += Indwidth;
+       }
+    }
+    free(ThisElement);
+
+/*  Read array of values:  */
+
+    if ( Type[0] != 'P' ) {          /* Skip if pattern only  */
+
+       if ( Type[0] == 'C' ) Nentries = 2*Nnzero;
+           else Nentries = Nnzero;
+
+    ThisElement = (char *) malloc(Valwidth+1);
+    if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
+    *(ThisElement+Valwidth) = (char) NULL;
+    count = 0;
+    for (i=0;i<Valcrd;i++)
+    {
+       fgets(line, BUFSIZ, in_file);
+       if ( sscanf(line,"%*s") < 0 )
+         IOHBTerminate("iohb.c: Null (or blank) line in value data region of HB file.\n");
+       if (Valflag == 'D')  {
+          while( strchr(line,'D') ) *strchr(line,'D') = 'E';
+/*           *strchr(Valfmt,'D') = 'E'; */
+       }
+       col =  0;
+       for (ind = 0;ind<Valperline;ind++)
+       {
+          if (count == Nentries) break;
+          strncpy(ThisElement,line+col,Valwidth);
+          /*ThisElement = substr(line,col,Valwidth);*/
+          if ( Valflag != 'F' && strchr(ThisElement,'E') == NULL ) {
+             /* insert a char prefix for exp */
+             last = strlen(ThisElement);
+             for (j=last+1;j>=0;j--) {
+                ThisElement[j] = ThisElement[j-1];
+                if ( ThisElement[j] == '+' || ThisElement[j] == '-' ) {
+                   ThisElement[j-1] = Valflag;
+                   break;
+                }
+             }
+          }
+          val[count] = atof(ThisElement);
+          count++; col += Valwidth;
+       }
+    }
+    free(ThisElement);
+    }
+
+    fclose(in_file);
+    return 1;
+}
+
+int readHB_newmat_double(const char* filename, int* M, int* N, int* nonzeros,
+                         int** colptr, int** rowind, double** val)
+{
+        int Nrhs;
+        char *Type;
+
+        readHB_info(filename, M, N, nonzeros, &Type, &Nrhs);
+
+        *colptr = (int *)malloc((*N+1)*sizeof(int));
+        if ( *colptr == NULL ) IOHBTerminate("Insufficient memory for colptr.\n");
+        *rowind = (int *)malloc(*nonzeros*sizeof(int));
+        if ( *rowind == NULL ) IOHBTerminate("Insufficient memory for rowind.\n");
+        if ( Type[0] == 'C' ) {
+/*
+   fprintf(stderr, "Warning: Reading complex data from HB file %s.\n",filename);
+   fprintf(stderr, "         Real and imaginary parts will be interlaced in val[].\n");
+*/
+           /* Malloc enough space for real AND imaginary parts of val[] */
+           *val = (double *)malloc(*nonzeros*sizeof(double)*2);
+           if ( *val == NULL ) IOHBTerminate("Insufficient memory for val.\n");
+        } else {
+           if ( Type[0] != 'P' ) {
+             /* Malloc enough space for real array val[] */
+             *val = (double *)malloc(*nonzeros*sizeof(double));
+             if ( *val == NULL ) IOHBTerminate("Insufficient memory for val.\n");
+           }
+        }  /* No val[] space needed if pattern only */
+        return readHB_mat_double(filename, *colptr, *rowind, *val);
+
+}
+
+int readHB_aux_double(const char* filename, const char AuxType, double b[])
+{
+/****************************************************************************/
+/*  This function opens and reads the specified file, placing auxillary     */
+/*  vector(s) of the given type (if available) in b.                        */
+/*  Return value is the number of vectors successfully read.                */
+/*                                                                          */
+/*                AuxType = 'F'   full right-hand-side vector(s)            */
+/*                AuxType = 'G'   initial Guess vector(s)                   */
+/*                AuxType = 'X'   eXact solution vector(s)                  */
+/*                                                                          */
+/*    ----------                                                            */
+/*    **CAVEAT**                                                            */
+/*    ----------                                                            */
+/*  Parsing real formats from Fortran is tricky, and this file reader       */
+/*  does not claim to be foolproof.   It has been tested for cases when     */
+/*  the real values are printed consistently and evenly spaced on each      */
+/*  line, with Fixed (F), and Exponential (E or D) formats.                 */
+/*                                                                          */
+/*  **  If the input file does not adhere to the H/B format, the  **        */
+/*  **             results will be unpredictable.                 **        */
+/*                                                                          */
+/****************************************************************************/
+    FILE *in_file;
+    int i,j,n,maxcol,start,stride,col,last,linel;
+    int Ptrcrd, Indcrd, Valcrd, Rhscrd;
+    int Nrow, Ncol, Nnzero, Nentries;
+    int Nrhs, nvecs, rhsi;
+    int Rhsperline, Rhswidth, Rhsprec;
+    int Rhsflag;
+    char *ThisElement;
+    char Title[73], Key[9], Type[4], Rhstype[4];
+    char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
+    char line[BUFSIZ];
+
+    if ((in_file = fopen( filename, "r")) == NULL) {
+      fprintf(stderr,"Error: Cannot open file: %s\n",filename);
+      return 0;
+     }
+
+    readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
+                  Ptrfmt, Indfmt, Valfmt, Rhsfmt,
+                  &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
+
+    if (Nrhs <= 0)
+    {
+      fprintf(stderr, "Warn: Attempt to read auxillary vector(s) when none are present.\n");
+      return 0;
+    }
+    if (Rhstype[0] != 'F' )
+    {
+      fprintf(stderr,"Warn: Attempt to read auxillary vector(s) which are not stored in Full form.\n");
+      fprintf(stderr,"       Rhs must be specified as full. \n");
+      return 0;
+    }
+
+/* If reading complex data, allow for interleaved real and imaginary values. */
+    if ( Type[0] == 'C' ) {
+       Nentries = 2*Nrow;
+     } else {
+       Nentries = Nrow;
+    }
+
+    nvecs = 1;
+
+    if ( Rhstype[1] == 'G' ) nvecs++;
+    if ( Rhstype[2] == 'X' ) nvecs++;
+
+    if ( AuxType == 'G' && Rhstype[1] != 'G' ) {
+      fprintf(stderr, "Warn: Attempt to read auxillary Guess vector(s) when none are present.\n");
+      return 0;
+    }
+    if ( AuxType == 'X' && Rhstype[2] != 'X' ) {
+      fprintf(stderr, "Warn: Attempt to read auxillary eXact solution vector(s) when none are present.\n");
+      return 0;
+    }
+
+    ParseRfmt(Rhsfmt, &Rhsperline, &Rhswidth, &Rhsprec,&Rhsflag);
+    maxcol = Rhsperline*Rhswidth;
+
+/*  Lines to skip before starting to read RHS values... */
+    n = Ptrcrd + Indcrd + Valcrd;
+
+    for (i = 0; i < n; i++)
+      fgets(line, BUFSIZ, in_file);
+
+/*  start  - number of initial aux vector entries to skip   */
+/*           to reach first  vector requested               */
+/*  stride - number of aux vector entries to skip between   */
+/*           requested vectors                              */
+    if ( AuxType == 'F' ) start = 0;
+    else if ( AuxType == 'G' ) start = Nentries;
+    else start = (nvecs-1)*Nentries;
+    stride = (nvecs-1)*Nentries;
+
+    fgets(line, BUFSIZ, in_file);
+    linel= strchr(line,'\n')-line;
+    col = 0;
+/*  Skip to initial offset */
+
+    for (i=0;i<start;i++) {
+       if ( col >=  ( maxcol<linel?maxcol:linel ) ) {
+           fgets(line, BUFSIZ, in_file);
+           linel= strchr(line,'\n')-line;
+           col = 0;
+       }
+       col += Rhswidth;
+    }
+    if (Rhsflag == 'D')  {
+       while( strchr(line,'D') ) *strchr(line,'D') = 'E';
+    }
+
+/*  Read a vector of desired type, then skip to next */
+/*  repeating to fill Nrhs vectors                   */
+
+  ThisElement = (char *) malloc(Rhswidth+1);
+  if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
+  *(ThisElement+Rhswidth) = (char) NULL;
+  for (rhsi=0;rhsi<Nrhs;rhsi++) {
+
+    for (i=0;i<Nentries;i++) {
+       if ( col >= ( maxcol<linel?maxcol:linel ) ) {
+           fgets(line, BUFSIZ, in_file);
+           linel= strchr(line,'\n')-line;
+           if (Rhsflag == 'D')  {
+              while( strchr(line,'D') ) *strchr(line,'D') = 'E';
+           }
+           col = 0;
+       }
+       strncpy(ThisElement,line+col,Rhswidth);
+       /*ThisElement = substr(line, col, Rhswidth);*/
+          if ( Rhsflag != 'F' && strchr(ThisElement,'E') == NULL ) {
+             /* insert a char prefix for exp */
+             last = strlen(ThisElement);
+             for (j=last+1;j>=0;j--) {
+                ThisElement[j] = ThisElement[j-1];
+                if ( ThisElement[j] == '+' || ThisElement[j] == '-' ) {
+                   ThisElement[j-1] = Rhsflag;
+                   break;
+                }
+             }
+          }
+       b[i] = atof(ThisElement);
+       col += Rhswidth;
+    }
+
+/*  Skip any interleaved Guess/eXact vectors */
+
+    for (i=0;i<stride;i++) {
+       if ( col >= ( maxcol<linel?maxcol:linel ) ) {
+           fgets(line, BUFSIZ, in_file);
+           linel= strchr(line,'\n')-line;
+           col = 0;
+       }
+       col += Rhswidth;
+    }
+
+  }
+  free(ThisElement);
+
+
+    fclose(in_file);
+    return Nrhs;
+}
+
+int readHB_newaux_double(const char* filename, const char AuxType, double** b)
+{
+        int Nrhs,M,N,nonzeros;
+        char *Type;
+
+        readHB_info(filename, &M, &N, &nonzeros, &Type, &Nrhs);
+        if ( Nrhs <= 0 ) {
+          fprintf(stderr,"Warn: Requested read of aux vector(s) when none are present.\n");
+          return 0;
+        } else {
+          if ( Type[0] == 'C' ) {
+            fprintf(stderr, "Warning: Reading complex aux vector(s) from HB file %s.",filename);
+            fprintf(stderr, "         Real and imaginary parts will be interlaced in b[].");
+            *b = (double *)malloc(M*Nrhs*sizeof(double)*2);
+            if ( *b == NULL ) IOHBTerminate("Insufficient memory for rhs.\n");
+            return readHB_aux_double(filename, AuxType, *b);
+          } else {
+            *b = (double *)malloc(M*Nrhs*sizeof(double));
+            if ( *b == NULL ) IOHBTerminate("Insufficient memory for rhs.\n");
+            return readHB_aux_double(filename, AuxType, *b);
+          }
+        }
+}
+
+int writeHB_mat_double(const char* filename, int M, int N,
+                        int nz, const int colptr[], const int rowind[],
+                        const double val[], int Nrhs, const double rhs[],
+                        const double guess[], const double exact[],
+                        const char* Title, const char* Key, const char* Type,
+                        char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
+                        const char* Rhstype)
+{
+/****************************************************************************/
+/*  The writeHB function opens the named file and writes the specified      */
+/*  matrix and optional right-hand-side(s) to that file in Harwell-Boeing   */
+/*  format.                                                                 */
+/*                                                                          */
+/*  For a description of the Harwell Boeing standard, see:                  */
+/*            Duff, et al.,  ACM TOMS Vol.15, No.1, March 1989              */
+/*                                                                          */
+/****************************************************************************/
+    FILE *out_file;
+    int i,j,entry,offset,acount,linemod;
+    int totcrd, ptrcrd, indcrd, valcrd, rhscrd;
+    int nvalentries, nrhsentries;
+    int Ptrperline, Ptrwidth, Indperline, Indwidth;
+    int Rhsperline, Rhswidth, Rhsprec;
+    int Rhsflag;
+    int Valperline, Valwidth, Valprec;
+    int Valflag;           /* Indicates 'E','D', or 'F' float format */
+    char pformat[16],iformat[16],vformat[19],rformat[19];
+
+    if ( Type[0] == 'C' ) {
+         nvalentries = 2*nz;
+         nrhsentries = 2*M;
+    } else {
+         nvalentries = nz;
+         nrhsentries = M;
+    }
+
+    if ( filename != NULL ) {
+       if ( (out_file = fopen( filename, "w")) == NULL ) {
+         fprintf(stderr,"Error: Cannot open file: %s\n",filename);
+         return 0;
+       }
+    } else out_file = stdout;
+
+    if ( Ptrfmt == NULL ) Ptrfmt = "(8I10)";
+    ParseIfmt(Ptrfmt,&Ptrperline,&Ptrwidth);
+    sprintf(pformat,"%%%dd",Ptrwidth);
+    ptrcrd = (N+1)/Ptrperline;
+    if ( (N+1)%Ptrperline != 0) ptrcrd++;
+
+    if ( Indfmt == NULL ) Indfmt =  Ptrfmt;
+    ParseIfmt(Indfmt,&Indperline,&Indwidth);
+    sprintf(iformat,"%%%dd",Indwidth);
+    indcrd = nz/Indperline;
+    if ( nz%Indperline != 0) indcrd++;
+
+    if ( Type[0] != 'P' ) {          /* Skip if pattern only  */
+      if ( Valfmt == NULL ) Valfmt = "(4E20.13)";
+      ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
+      if (Valflag == 'D') *strchr(Valfmt,'D') = 'E';
+      if (Valflag == 'F')
+         sprintf(vformat,"%% %d.%df",Valwidth,Valprec);
+      else
+         sprintf(vformat,"%% %d.%dE",Valwidth,Valprec);
+      valcrd = nvalentries/Valperline;
+      if ( nvalentries%Valperline != 0) valcrd++;
+    } else valcrd = 0;
+
+    if ( Nrhs > 0 ) {
+       if ( Rhsfmt == NULL ) Rhsfmt = Valfmt;
+       ParseRfmt(Rhsfmt,&Rhsperline,&Rhswidth,&Rhsprec, &Rhsflag);
+       if (Rhsflag == 'F')
+          sprintf(rformat,"%% %d.%df",Rhswidth,Rhsprec);
+       else
+          sprintf(rformat,"%% %d.%dE",Rhswidth,Rhsprec);
+       if (Rhsflag == 'D') *strchr(Rhsfmt,'D') = 'E';
+       rhscrd = nrhsentries/Rhsperline;
+       if ( nrhsentries%Rhsperline != 0) rhscrd++;
+       if ( Rhstype[1] == 'G' ) rhscrd+=rhscrd;
+       if ( Rhstype[2] == 'X' ) rhscrd+=rhscrd;
+       rhscrd*=Nrhs;
+    } else rhscrd = 0;
+
+    totcrd = 4+ptrcrd+indcrd+valcrd+rhscrd;
+
+
+/*  Print header information:  */
+
+    fprintf(out_file,"%-72s%-8s\n%14d%14d%14d%14d%14d\n",Title, Key, totcrd,
+            ptrcrd, indcrd, valcrd, rhscrd);
+    fprintf(out_file,"%3s%11s%14d%14d%14d\n",Type,"          ", M, N, nz);
+    fprintf(out_file,"%-16s%-16s%-20s", Ptrfmt, Indfmt, Valfmt);
+    if ( Nrhs != 0 ) {
+/*     Print Rhsfmt on fourth line and                                    */
+/*           optional fifth header line for auxillary vector information: */
+       fprintf(out_file,"%-20s\n%-14s%d\n",Rhsfmt,Rhstype,Nrhs);
+    } else fprintf(out_file,"\n");
+
+    offset = 1-_SP_base;  /* if base 0 storage is declared (via macro definition), */
+                          /* then storage entries are offset by 1                  */
+
+/*  Print column pointers:   */
+    for (i=0;i<N+1;i++)
+    {
+       entry = colptr[i]+offset;
+       fprintf(out_file,pformat,entry);
+       if ( (i+1)%Ptrperline == 0 ) fprintf(out_file,"\n");
+    }
+
+   if ( (N+1) % Ptrperline != 0 ) fprintf(out_file,"\n");
+
+/*  Print row indices:       */
+    for (i=0;i<nz;i++)
+    {
+       entry = rowind[i]+offset;
+       fprintf(out_file,iformat,entry);
+       if ( (i+1)%Indperline == 0 ) fprintf(out_file,"\n");
+    }
+
+   if ( nz % Indperline != 0 ) fprintf(out_file,"\n");
+
+/*  Print values:            */
+
+    if ( Type[0] != 'P' ) {          /* Skip if pattern only  */
+
+    for (i=0;i<nvalentries;i++)
+    {
+       fprintf(out_file,vformat,val[i]);
+       if ( (i+1)%Valperline == 0 ) fprintf(out_file,"\n");
+    }
+
+    if ( nvalentries % Valperline != 0 ) fprintf(out_file,"\n");
+
+/*  If available,  print right hand sides,
+           guess vectors and exact solution vectors:  */
+    acount = 1;
+    linemod = 0;
+    if ( Nrhs > 0 ) {
+       for (i=0;i<Nrhs;i++)
+       {
+          for ( j=0;j<nrhsentries;j++ ) {
+            fprintf(out_file,rformat,rhs[j]);
+            if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
+          }
+          if ( acount%Rhsperline != linemod ) {
+            fprintf(out_file,"\n");
+            linemod = (acount-1)%Rhsperline;
+          }
+          rhs += nrhsentries;
+          if ( Rhstype[1] == 'G' ) {
+            for ( j=0;j<nrhsentries;j++ ) {
+              fprintf(out_file,rformat,guess[j]);
+              if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
+            }
+            if ( acount%Rhsperline != linemod ) {
+              fprintf(out_file,"\n");
+              linemod = (acount-1)%Rhsperline;
+            }
+            guess += nrhsentries;
+          }
+          if ( Rhstype[2] == 'X' ) {
+            for ( j=0;j<nrhsentries;j++ ) {
+              fprintf(out_file,rformat,exact[j]);
+              if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
+            }
+            if ( acount%Rhsperline != linemod ) {
+              fprintf(out_file,"\n");
+              linemod = (acount-1)%Rhsperline;
+            }
+            exact += nrhsentries;
+          }
+       }
+    }
+
+    }
+
+    if ( fclose(out_file) != 0){
+      fprintf(stderr,"Error closing file in writeHB_mat_double().\n");
+      return 0;
+    } else return 1;
+
+}
+
+int readHB_mat_char(const char* filename, int colptr[], int rowind[],
+                                           char val[], char* Valfmt)
+{
+/****************************************************************************/
+/*  This function opens and reads the specified file, interpreting its      */
+/*  contents as a sparse matrix stored in the Harwell/Boeing standard       */
+/*  format and creating compressed column storage scheme vectors to hold    */
+/*  the index and nonzero value information.                                */
+/*                                                                          */
+/*    ----------                                                            */
+/*    **CAVEAT**                                                            */
+/*    ----------                                                            */
+/*  Parsing real formats from Fortran is tricky, and this file reader       */
+/*  does not claim to be foolproof.   It has been tested for cases when     */
+/*  the real values are printed consistently and evenly spaced on each      */
+/*  line, with Fixed (F), and Exponential (E or D) formats.                 */
+/*                                                                          */
+/*  **  If the input file does not adhere to the H/B format, the  **        */
+/*  **             results will be unpredictable.                 **        */
+/*                                                                          */
+/****************************************************************************/
+    FILE *in_file;
+    int i,j,ind,col,offset,count,last;
+    int Nrow,Ncol,Nnzero,Nentries,Nrhs;
+    int Ptrcrd, Indcrd, Valcrd, Rhscrd;
+    int Ptrperline, Ptrwidth, Indperline, Indwidth;
+    int Valperline, Valwidth, Valprec;
+    int Valflag;           /* Indicates 'E','D', or 'F' float format */
+    char* ThisElement;
+    char line[BUFSIZ];
+    char Title[73], Key[8], Type[4], Rhstype[4];
+    char Ptrfmt[17], Indfmt[17], Rhsfmt[21];
+
+    if ( (in_file = fopen( filename, "r")) == NULL ) {
+       fprintf(stderr,"Error: Cannot open file: %s\n",filename);
+       return 0;
+    }
+
+    readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
+                  Ptrfmt, Indfmt, Valfmt, Rhsfmt,
+                  &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
+
+/*  Parse the array input formats from Line 3 of HB file  */
+    ParseIfmt(Ptrfmt,&Ptrperline,&Ptrwidth);
+    ParseIfmt(Indfmt,&Indperline,&Indwidth);
+    if ( Type[0] != 'P' ) {          /* Skip if pattern only  */
+       ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
+       if (Valflag == 'D') {
+          *strchr(Valfmt,'D') = 'E';
+       }
+    }
+
+/*  Read column pointer array:   */
+
+    offset = 1-_SP_base;  /* if base 0 storage is declared (via macro definition), */
+                          /* then storage entries are offset by 1                  */
+
+    ThisElement = (char *) malloc(Ptrwidth+1);
+    if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
+    *(ThisElement+Ptrwidth) = (char) NULL;
+    count=0;
+    for (i=0;i<Ptrcrd;i++)
+    {
+       fgets(line, BUFSIZ, in_file);
+       if ( sscanf(line,"%*s") < 0 )
+         IOHBTerminate("iohb.c: Null (or blank) line in pointer data region of HB file.\n");
+       col =  0;
+       for (ind = 0;ind<Ptrperline;ind++)
+       {
+          if (count > Ncol) break;
+          strncpy(ThisElement,line+col,Ptrwidth);
+          /*ThisElement = substr(line,col,Ptrwidth);*/
+          colptr[count] = atoi(ThisElement)-offset;
+          count++; col += Ptrwidth;
+       }
+    }
+    free(ThisElement);
+
+/*  Read row index array:  */
+
+    ThisElement = (char *) malloc(Indwidth+1);
+    if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
+    *(ThisElement+Indwidth) = (char) NULL;
+    count = 0;
+    for (i=0;i<Indcrd;i++)
+    {
+       fgets(line, BUFSIZ, in_file);
+       if ( sscanf(line,"%*s") < 0 )
+         IOHBTerminate("iohb.c: Null (or blank) line in index data region of HB file.\n");
+       col =  0;
+       for (ind = 0;ind<Indperline;ind++)
+       {
+          if (count == Nnzero) break;
+          strncpy(ThisElement,line+col,Indwidth);
+          /*ThisElement = substr(line,col,Indwidth);*/
+          rowind[count] = atoi(ThisElement)-offset;
+          count++; col += Indwidth;
+       }
+    }
+    free(ThisElement);
+
+/*  Read array of values:  AS CHARACTERS*/
+
+    if ( Type[0] != 'P' ) {          /* Skip if pattern only  */
+
+       if ( Type[0] == 'C' ) Nentries = 2*Nnzero;
+           else Nentries = Nnzero;
+
+    ThisElement = (char *) malloc(Valwidth+1);
+    if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
+    *(ThisElement+Valwidth) = (char) NULL;
+    count = 0;
+    for (i=0;i<Valcrd;i++)
+    {
+       fgets(line, BUFSIZ, in_file);
+       if ( sscanf(line,"%*s") < 0 )
+         IOHBTerminate("iohb.c: Null (or blank) line in value data region of HB file.\n");
+       if (Valflag == 'D') {
+          while( strchr(line,'D') ) *strchr(line,'D') = 'E';
+       }
+       col =  0;
+       for (ind = 0;ind<Valperline;ind++)
+       {
+          if (count == Nentries) break;
+          ThisElement = &val[count*Valwidth];
+          strncpy(ThisElement,line+col,Valwidth);
+          /*strncpy(ThisElement,substr(line,col,Valwidth),Valwidth);*/
+          if ( Valflag != 'F' && strchr(ThisElement,'E') == NULL ) {
+             /* insert a char prefix for exp */
+             last = strlen(ThisElement);
+             for (j=last+1;j>=0;j--) {
+                ThisElement[j] = ThisElement[j-1];
+                if ( ThisElement[j] == '+' || ThisElement[j] == '-' ) {
+                   ThisElement[j-1] = Valflag;
+                   break;
+                }
+             }
+          }
+          count++; col += Valwidth;
+       }
+    }
+    }
+
+    return 1;
+}
+
+int readHB_newmat_char(const char* filename, int* M, int* N, int* nonzeros, int** colptr,
+                          int** rowind, char** val, char** Valfmt)
+{
+    FILE *in_file;
+    int Nrhs;
+    int Ptrcrd, Indcrd, Valcrd, Rhscrd;
+    int Valperline, Valwidth, Valprec;
+    int Valflag;           /* Indicates 'E','D', or 'F' float format */
+    char Title[73], Key[9], Type[4], Rhstype[4];
+    char Ptrfmt[17], Indfmt[17], Rhsfmt[21];
+
+    if ((in_file = fopen( filename, "r")) == NULL) {
+      fprintf(stderr,"Error: Cannot open file: %s\n",filename);
+      return 0;
+     }
+
+    *Valfmt = (char *)malloc(21*sizeof(char));
+    if ( *Valfmt == NULL ) IOHBTerminate("Insufficient memory for Valfmt.");
+    readHB_header(in_file, Title, Key, Type, M, N, nonzeros, &Nrhs,
+                  Ptrfmt, Indfmt, (*Valfmt), Rhsfmt,
+                  &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
+    fclose(in_file);
+    ParseRfmt(*Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
+
+        *colptr = (int *)malloc((*N+1)*sizeof(int));
+        if ( *colptr == NULL ) IOHBTerminate("Insufficient memory for colptr.\n");
+        *rowind = (int *)malloc(*nonzeros*sizeof(int));
+        if ( *rowind == NULL ) IOHBTerminate("Insufficient memory for rowind.\n");
+        if ( Type[0] == 'C' ) {
+/*
+   fprintf(stderr, "Warning: Reading complex data from HB file %s.\n",filename);
+   fprintf(stderr, "         Real and imaginary parts will be interlaced in val[].\n");
+*/
+           /* Malloc enough space for real AND imaginary parts of val[] */
+           *val = (char *)malloc(*nonzeros*Valwidth*sizeof(char)*2);
+           if ( *val == NULL ) IOHBTerminate("Insufficient memory for val.\n");
+        } else {
+           if ( Type[0] != 'P' ) {
+             /* Malloc enough space for real array val[] */
+             *val = (char *)malloc(*nonzeros*Valwidth*sizeof(char));
+             if ( *val == NULL ) IOHBTerminate("Insufficient memory for val.\n");
+           }
+        }  /* No val[] space needed if pattern only */
+        return readHB_mat_char(filename, *colptr, *rowind, *val, *Valfmt);
+
+}
+
+int readHB_aux_char(const char* filename, const char AuxType, char b[])
+{
+/****************************************************************************/
+/*  This function opens and reads the specified file, placing auxilary      */
+/*  vector(s) of the given type (if available) in b :                       */
+/*  Return value is the number of vectors successfully read.                */
+/*                                                                          */
+/*                AuxType = 'F'   full right-hand-side vector(s)            */
+/*                AuxType = 'G'   initial Guess vector(s)                   */
+/*                AuxType = 'X'   eXact solution vector(s)                  */
+/*                                                                          */
+/*    ----------                                                            */
+/*    **CAVEAT**                                                            */
+/*    ----------                                                            */
+/*  Parsing real formats from Fortran is tricky, and this file reader       */
+/*  does not claim to be foolproof.   It has been tested for cases when     */
+/*  the real values are printed consistently and evenly spaced on each      */
+/*  line, with Fixed (F), and Exponential (E or D) formats.                 */
+/*                                                                          */
+/*  **  If the input file does not adhere to the H/B format, the  **        */
+/*  **             results will be unpredictable.                 **        */
+/*                                                                          */
+/****************************************************************************/
+    FILE *in_file;
+    int i,j,n,maxcol,start,stride,col,last,linel,nvecs,rhsi;
+    int Nrow, Ncol, Nnzero, Nentries,Nrhs;
+    int Ptrcrd, Indcrd, Valcrd, Rhscrd;
+    int Rhsperline, Rhswidth, Rhsprec;
+    int Rhsflag;
+    char Title[73], Key[9], Type[4], Rhstype[4];
+    char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
+    char line[BUFSIZ];
+    char *ThisElement;
+
+    if ((in_file = fopen( filename, "r")) == NULL) {
+      fprintf(stderr,"Error: Cannot open file: %s\n",filename);
+      return 0;
+     }
+
+    readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
+                  Ptrfmt, Indfmt, Valfmt, Rhsfmt,
+                  &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
+
+    if (Nrhs <= 0)
+    {
+      fprintf(stderr, "Warn: Attempt to read auxillary vector(s) when none are present.\n");
+      return 0;
+    }
+    if (Rhstype[0] != 'F' )
+    {
+      fprintf(stderr,"Warn: Attempt to read auxillary vector(s) which are not stored in Full form.\n");
+      fprintf(stderr,"       Rhs must be specified as full. \n");
+      return 0;
+    }
+
+/* If reading complex data, allow for interleaved real and imaginary values. */
+    if ( Type[0] == 'C' ) {
+       Nentries = 2*Nrow;
+     } else {
+       Nentries = Nrow;
+    }
+
+    nvecs = 1;
+
+    if ( Rhstype[1] == 'G' ) nvecs++;
+    if ( Rhstype[2] == 'X' ) nvecs++;
+
+    if ( AuxType == 'G' && Rhstype[1] != 'G' ) {
+      fprintf(stderr, "Warn: Attempt to read auxillary Guess vector(s) when none are present.\n");
+      return 0;
+    }
+    if ( AuxType == 'X' && Rhstype[2] != 'X' ) {
+      fprintf(stderr, "Warn: Attempt to read auxillary eXact solution vector(s) when none are present.\n");
+      return 0;
+    }
+
+    ParseRfmt(Rhsfmt, &Rhsperline, &Rhswidth, &Rhsprec,&Rhsflag);
+    maxcol = Rhsperline*Rhswidth;
+
+/*  Lines to skip before starting to read RHS values... */
+    n = Ptrcrd + Indcrd + Valcrd;
+
+    for (i = 0; i < n; i++)
+      fgets(line, BUFSIZ, in_file);
+
+/*  start  - number of initial aux vector entries to skip   */
+/*           to reach first  vector requested               */
+/*  stride - number of aux vector entries to skip between   */
+/*           requested vectors                              */
+    if ( AuxType == 'F' ) start = 0;
+    else if ( AuxType == 'G' ) start = Nentries;
+    else start = (nvecs-1)*Nentries;
+    stride = (nvecs-1)*Nentries;
+
+    fgets(line, BUFSIZ, in_file);
+    linel= strchr(line,'\n')-line;
+    if ( sscanf(line,"%*s") < 0 )
+       IOHBTerminate("iohb.c: Null (or blank) line in auxillary vector data region of HB file.\n");
+    col = 0;
+/*  Skip to initial offset */
+
+    for (i=0;i<start;i++) {
+       col += Rhswidth;
+       if ( col >= ( maxcol<linel?maxcol:linel ) ) {
+           fgets(line, BUFSIZ, in_file);
+           linel= strchr(line,'\n')-line;
+       if ( sscanf(line,"%*s") < 0 )
+       IOHBTerminate("iohb.c: Null (or blank) line in auxillary vector data region of HB file.\n");
+           col = 0;
+       }
+    }
+
+    if (Rhsflag == 'D')  {
+      while( strchr(line,'D') ) *strchr(line,'D') = 'E';
+    }
+/*  Read a vector of desired type, then skip to next */
+/*  repeating to fill Nrhs vectors                   */
+
+  for (rhsi=0;rhsi<Nrhs;rhsi++) {
+
+    for (i=0;i<Nentries;i++) {
+       if ( col >= ( maxcol<linel?maxcol:linel ) ) {
+           fgets(line, BUFSIZ, in_file);
+           linel= strchr(line,'\n')-line;
+       if ( sscanf(line,"%*s") < 0 )
+       IOHBTerminate("iohb.c: Null (or blank) line in auxillary vector data region of HB file.\n");
+           if (Rhsflag == 'D')  {
+              while( strchr(line,'D') ) *strchr(line,'D') = 'E';
+           }
+           col = 0;
+       }
+       ThisElement = &b[i*Rhswidth];
+       strncpy(ThisElement,line+col,Rhswidth);
+          if ( Rhsflag != 'F' && strchr(ThisElement,'E') == NULL ) {
+             /* insert a char prefix for exp */
+             last = strlen(ThisElement);
+             for (j=last+1;j>=0;j--) {
+                ThisElement[j] = ThisElement[j-1];
+                if ( ThisElement[j] == '+' || ThisElement[j] == '-' ) {
+                   ThisElement[j-1] = Rhsflag;
+                   break;
+                }
+             }
+          }
+       col += Rhswidth;
+    }
+    b+=Nentries*Rhswidth;
+
+/*  Skip any interleaved Guess/eXact vectors */
+
+    for (i=0;i<stride;i++) {
+       col += Rhswidth;
+       if ( col >= ( maxcol<linel?maxcol:linel ) ) {
+           fgets(line, BUFSIZ, in_file);
+           linel= strchr(line,'\n')-line;
+       if ( sscanf(line,"%*s") < 0 )
+       IOHBTerminate("iohb.c: Null (or blank) line in auxillary vector data region of HB file.\n");
+           col = 0;
+       }
+    }
+
+  }
+
+
+    fclose(in_file);
+    return Nrhs;
+}
+
+int readHB_newaux_char(const char* filename, const char AuxType, char** b, char** Rhsfmt)
+{
+    FILE *in_file;
+    int Ptrcrd, Indcrd, Valcrd, Rhscrd;
+    int Nrow,Ncol,Nnzero,Nrhs;
+    int Rhsperline, Rhswidth, Rhsprec;
+    int Rhsflag;
+    char Title[73], Key[9], Type[4], Rhstype[4];
+    char Ptrfmt[17], Indfmt[17], Valfmt[21];
+
+    if ((in_file = fopen( filename, "r")) == NULL) {
+      fprintf(stderr,"Error: Cannot open file: %s\n",filename);
+      return 0;
+     }
+
+    *Rhsfmt = (char *)malloc(21*sizeof(char));
+    if ( *Rhsfmt == NULL ) IOHBTerminate("Insufficient memory for Rhsfmt.");
+    readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
+                  Ptrfmt, Indfmt, Valfmt, (*Rhsfmt),
+                  &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
+     fclose(in_file);
+        if ( Nrhs == 0 ) {
+          fprintf(stderr,"Warn: Requested read of aux vector(s) when none are present.\n");
+          return 0;
+        } else {
+          ParseRfmt(*Rhsfmt,&Rhsperline,&Rhswidth,&Rhsprec,&Rhsflag);
+          if ( Type[0] == 'C' ) {
+            fprintf(stderr, "Warning: Reading complex aux vector(s) from HB file %s.",filename);
+            fprintf(stderr, "         Real and imaginary parts will be interlaced in b[].");
+            *b = (char *)malloc(Nrow*Nrhs*Rhswidth*sizeof(char)*2);
+            if ( *b == NULL ) IOHBTerminate("Insufficient memory for rhs.\n");
+            return readHB_aux_char(filename, AuxType, *b);
+          } else {
+            *b = (char *)malloc(Nrow*Nrhs*Rhswidth*sizeof(char));
+            if ( *b == NULL ) IOHBTerminate("Insufficient memory for rhs.\n");
+            return readHB_aux_char(filename, AuxType, *b);
+          }
+        }
+}
+
+int writeHB_mat_char(const char* filename, int M, int N,
+                        int nz, const int colptr[], const int rowind[],
+                        const char val[], int Nrhs, const char rhs[],
+                        const char guess[], const char exact[],
+                        const char* Title, const char* Key, const char* Type,
+                        char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
+                        const char* Rhstype)
+{
+/****************************************************************************/
+/*  The writeHB function opens the named file and writes the specified      */
+/*  matrix and optional right-hand-side(s) to that file in Harwell-Boeing   */
+/*  format.                                                                 */
+/*                                                                          */
+/*  For a description of the Harwell Boeing standard, see:                  */
+/*            Duff, et al.,  ACM TOMS Vol.15, No.1, March 1989              */
+/*                                                                          */
+/****************************************************************************/
+    FILE *out_file;
+    int i,j,acount,linemod,entry,offset;
+    int totcrd, ptrcrd, indcrd, valcrd, rhscrd;
+    int nvalentries, nrhsentries;
+    int Ptrperline, Ptrwidth, Indperline, Indwidth;
+    int Rhsperline, Rhswidth, Rhsprec;
+    int Rhsflag;
+    int Valperline = 1, Valwidth, Valprec;
+    int Valflag;           /* Indicates 'E','D', or 'F' float format */
+    char pformat[16],iformat[16],vformat[19],rformat[19];
+
+    if ( Type[0] == 'C' ) {
+         nvalentries = 2*nz;
+         nrhsentries = 2*M;
+    } else {
+         nvalentries = nz;
+         nrhsentries = M;
+    }
+
+    if ( filename != NULL ) {
+       if ( (out_file = fopen( filename, "w")) == NULL ) {
+         fprintf(stderr,"Error: Cannot open file: %s\n",filename);
+         return 0;
+       }
+    } else out_file = stdout;
+
+    if ( Ptrfmt == NULL ) Ptrfmt = "(8I10)";
+    ParseIfmt(Ptrfmt,&Ptrperline,&Ptrwidth);
+    sprintf(pformat,"%%%dd",Ptrwidth);
+
+    if ( Indfmt == NULL ) Indfmt =  Ptrfmt;
+    ParseIfmt(Indfmt,&Indperline,&Indwidth);
+    sprintf(iformat,"%%%dd",Indwidth);
+
+    if ( Type[0] != 'P' ) {          /* Skip if pattern only  */
+      if ( Valfmt == NULL ) Valfmt = "(4E20.13)";
+      ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
+      sprintf(vformat,"%%%ds",Valwidth);
+    }
+
+    ptrcrd = (N+1)/Ptrperline;
+    if ( (N+1)%Ptrperline != 0) ptrcrd++;
+
+    indcrd = nz/Indperline;
+    if ( nz%Indperline != 0) indcrd++;
+
+    valcrd = nvalentries/Valperline;
+    if ( nvalentries%Valperline != 0) valcrd++;
+
+    if ( Nrhs > 0 ) {
+       if ( Rhsfmt == NULL ) Rhsfmt = Valfmt;
+       ParseRfmt(Rhsfmt,&Rhsperline,&Rhswidth,&Rhsprec, &Rhsflag);
+       sprintf(rformat,"%%%ds",Rhswidth);
+       rhscrd = nrhsentries/Rhsperline;
+       if ( nrhsentries%Rhsperline != 0) rhscrd++;
+       if ( Rhstype[1] == 'G' ) rhscrd+=rhscrd;
+       if ( Rhstype[2] == 'X' ) rhscrd+=rhscrd;
+       rhscrd*=Nrhs;
+    } else rhscrd = 0;
+
+    totcrd = 4+ptrcrd+indcrd+valcrd+rhscrd;
+
+
+/*  Print header information:  */
+
+    fprintf(out_file,"%-72s%-8s\n%14d%14d%14d%14d%14d\n",Title, Key, totcrd,
+            ptrcrd, indcrd, valcrd, rhscrd);
+    fprintf(out_file,"%3s%11s%14d%14d%14d\n",Type,"          ", M, N, nz);
+    fprintf(out_file,"%-16s%-16s%-20s", Ptrfmt, Indfmt, Valfmt);
+    if ( Nrhs != 0 ) {
+/*     Print Rhsfmt on fourth line and                                    */
+/*           optional fifth header line for auxillary vector information: */
+       fprintf(out_file,"%-20s\n%-14s%d\n",Rhsfmt,Rhstype,Nrhs);
+    } else fprintf(out_file,"\n");
+
+    offset = 1-_SP_base;  /* if base 0 storage is declared (via macro definition), */
+                          /* then storage entries are offset by 1                  */
+
+/*  Print column pointers:   */
+    for (i=0;i<N+1;i++)
+    {
+       entry = colptr[i]+offset;
+       fprintf(out_file,pformat,entry);
+       if ( (i+1)%Ptrperline == 0 ) fprintf(out_file,"\n");
+    }
+
+   if ( (N+1) % Ptrperline != 0 ) fprintf(out_file,"\n");
+
+/*  Print row indices:       */
+    for (i=0;i<nz;i++)
+    {
+       entry = rowind[i]+offset;
+       fprintf(out_file,iformat,entry);
+       if ( (i+1)%Indperline == 0 ) fprintf(out_file,"\n");
+    }
+
+   if ( nz % Indperline != 0 ) fprintf(out_file,"\n");
+
+/*  Print values:            */
+
+    if ( Type[0] != 'P' ) {          /* Skip if pattern only  */
+    for (i=0;i<nvalentries;i++)
+    {
+       fprintf(out_file,vformat,val+i*Valwidth);
+       if ( (i+1)%Valperline == 0 ) fprintf(out_file,"\n");
+    }
+
+    if ( nvalentries % Valperline != 0 ) fprintf(out_file,"\n");
+
+/*  Print right hand sides:  */
+    acount = 1;
+    linemod=0;
+    if ( Nrhs > 0 ) {
+      for (j=0;j<Nrhs;j++) {
+       for (i=0;i<nrhsentries;i++)
+       {
+          fprintf(out_file,rformat,rhs+i*Rhswidth);
+          if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
+       }
+       if ( acount%Rhsperline != linemod ) {
+          fprintf(out_file,"\n");
+          linemod = (acount-1)%Rhsperline;
+       }
+       if ( Rhstype[1] == 'G' ) {
+         for (i=0;i<nrhsentries;i++)
+         {
+           fprintf(out_file,rformat,guess+i*Rhswidth);
+           if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
+         }
+         if ( acount%Rhsperline != linemod ) {
+            fprintf(out_file,"\n");
+            linemod = (acount-1)%Rhsperline;
+         }
+       }
+       if ( Rhstype[2] == 'X' ) {
+         for (i=0;i<nrhsentries;i++)
+         {
+           fprintf(out_file,rformat,exact+i*Rhswidth);
+           if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
+         }
+         if ( acount%Rhsperline != linemod ) {
+            fprintf(out_file,"\n");
+            linemod = (acount-1)%Rhsperline;
+         }
+       }
+      }
+    }
+
+    }
+
+    if ( fclose(out_file) != 0){
+      fprintf(stderr,"Error closing file in writeHB_mat_char().\n");
+      return 0;
+    } else return 1;
+
+}
+
+int ParseIfmt(char* fmt, int* perline, int* width)
+{
+/*************************************************/
+/*  Parse an *integer* format field to determine */
+/*  width and number of elements per line.       */
+/*************************************************/
+    char *tmp;
+    if (fmt == NULL ) {
+      *perline = 0; *width = 0; return 0;
+    }
+    upcase(fmt);
+    tmp = strchr(fmt,'(');
+    tmp = substr(fmt,tmp - fmt + 1, strchr(fmt,'I') - tmp - 1);
+    *perline = atoi(tmp);
+    tmp = strchr(fmt,'I');
+    tmp = substr(fmt,tmp - fmt + 1, strchr(fmt,')') - tmp - 1);
+    return *width = atoi(tmp);
+}
+
+int ParseRfmt(char* fmt, int* perline, int* width, int* prec, int* flag)
+{
+/*************************************************/
+/*  Parse a *real* format field to determine     */
+/*  width and number of elements per line.       */
+/*  Also sets flag indicating 'E' 'F' 'P' or 'D' */
+/*  format.                                      */
+/*************************************************/
+    char* tmp;
+    char* tmp2;
+    char* tmp3;
+    int len;
+
+    if (fmt == NULL ) {
+      *perline = 0;
+      *width = 0;
+      flag = NULL;
+      return 0;
+    }
+
+    upcase(fmt);
+    if (strchr(fmt,'(') != NULL)  fmt = strchr(fmt,'(');
+    if (strchr(fmt,')') != NULL)  {
+       tmp2 = strchr(fmt,')');
+       while ( strchr(tmp2+1,')') != NULL ) {
+          tmp2 = strchr(tmp2+1,')');
+       }
+       *(tmp2+1) = (int) NULL;
+    }
+    if (strchr(fmt,'P') != NULL)  /* Remove any scaling factor, which */
+    {                             /* affects output only, not input */
+      if (strchr(fmt,'(') != NULL) {
+        tmp = strchr(fmt,'P');
+        if ( *(++tmp) == ',' ) tmp++;
+        tmp3 = strchr(fmt,'(')+1;
+        len = tmp-tmp3;
+        tmp2 = tmp3;
+        while ( *(tmp2+len) != (int) NULL ) {
+           *tmp2=*(tmp2+len);
+           tmp2++;
+        }
+        *(strchr(fmt,')')+1) = (int) NULL;
+      }
+    }
+    if (strchr(fmt,'E') != NULL) {
+       *flag = 'E';
+    } else if (strchr(fmt,'D') != NULL) {
+       *flag = 'D';
+    } else if (strchr(fmt,'F') != NULL) {
+       *flag = 'F';
+    } else {
+      fprintf(stderr,"Real format %s in H/B file not supported.\n",fmt);
+      return 0;
+    }
+    tmp = strchr(fmt,'(');
+    tmp = substr(fmt,tmp - fmt + 1, strchr(fmt,*flag) - tmp - 1);
+    *perline = atoi(tmp);
+    tmp = strchr(fmt,*flag);
+    if ( strchr(fmt,'.') ) {
+      *prec = atoi( substr( fmt, strchr(fmt,'.') - fmt + 1, strchr(fmt,')') - strchr(fmt,'.')-1) );
+      tmp = substr(fmt,tmp - fmt + 1, strchr(fmt,'.') - tmp - 1);
+    } else {
+      tmp = substr(fmt,tmp - fmt + 1, strchr(fmt,')') - tmp - 1);
+    }
+    return *width = atoi(tmp);
+}
+
+char* substr(const char* S, const int pos, const int len)
+{
+    int i;
+    char *SubS;
+    if ( pos+len <= strlen(S)) {
+    SubS = (char *)malloc(len+1);
+    if ( SubS == NULL ) IOHBTerminate("Insufficient memory for SubS.");
+    for (i=0;i<len;i++) SubS[i] = S[pos+i];
+    SubS[len] = (char) NULL;
+    } else {
+      SubS = NULL;
+    }
+    return SubS;
+}
+
+#include<ctype.h>
+void upcase(char* S)
+{
+/*  Convert S to uppercase     */
+    int i,len;
+    if ( S == NULL ) return;
+    len = strlen(S);
+    for (i=0;i< len;i++)
+       S[i] = toupper(S[i]);
+}
+
+void IOHBTerminate(const char* message)
+{
+   fprintf(stderr,"%s",message);
+   exit(1);
+}
+