iohb.c

http://gams.cam.nist.gov/acmd/Staff/RPozo/sparselib++.html

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*             ********   ***                                 SparseLib++    */
/*          *******  **  ***       ***      ***               v. 1.5c        */
/*           *****      ***     ******** ********                            */
/*            *****    ***     ******** ********              R. Pozo        */
/*       **  *******  ***   **   ***      ***                 K. Remington   */
/*        ********   ********                                 A. Lumsdaine   */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*                                                                           */
/*                                                                           */
/*                     SparseLib++ : Sparse Matrix Library                   */
/*                                                                           */
/*               National Institute of Standards and Technology              */
/*                        University of Notre Dame                           */
/*              Authors: R. Pozo, K. Remington, A. Lumsdaine                 */
/*                                                                           */
/*                                 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 notice appear in all copies and supporting        */
/* documentation.                                                            */
/*                                                                           */
/* Neither the Institutions (National Institute of Standards and Technology, */
/* University of Notre Dame) nor the Authors make any representations about  */
/* the suitability of this software for any purpose.  This software is       */
/* provided ``as is'' without expressed or implied warranty.                 */
/*                                                                           */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/


/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*      I/O for Harwell-Boeing files                                        */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*                                                                          */
/*  Several user I/O functions for Harwell-Boeing files are provided,       */
/*  together with various utility functions.                                */
/*  For a description of the Harwell-Boeing standard, see:                  */
/*                                                                          */
/*            Duff, et al.,  ACM TOMS Vol.15, No.1, March 1989              */
/*                                                                          */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Function:                                                               */
/*  void readHB_newmat_double(char *filename, int *M, int *N, *int nz,      */
/*        int **colptr, int **rowind,  double**val)                         */
/*                                                                          */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Function:                                                               */
/*  void readHB_info(char *filename, int *M, int *N, int *nz, 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), and the number of right-hand-sides stored along with   */
/*  the matrix (nrhs).                                                      */
/*  The optional verbose parameter, if set to 1,  triggers the output of    */
/*  more detailed information from the header, including title, etc.        */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Function:                                                               */
/*  void readHB_mat_double(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 and creating compressed column storage scheme vectors to hold    */
/*  the index and nonzero value information.                                */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Function:                                                               */
/*  void readHB_rhs_double(char *filename, double *b, int j)                */
/*                                                                          */
/*  Description:                                                            */
/*                                                                          */
/*  This function opens and reads the specified file, returning a right-    */
/*  hand-side vector b.  If the file specifies multiple right-hand-sides    */
/*  (that is, a right-hand-side matrix B), then the optional argument j     */
/*  may be used to select the (j+1)st column of B.                          */
/*                                                                          */
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/*  Function:                                                               */
/*  void writeHB_mat_double(char *filename, int M, int N, int nz,           */
/*         *colptr,  int *rowind, double *val, int nrhs, double *rhs,       */
/*         char *Title, char *Key)                                          */
/*                                                                          */
/*  Description:                                                            */
/*                                                                          */
/*  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.                                                                 */
/*                                                                          */
/****************************************************************************/


#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<malloc.h>

void nullchk();
void ParseIfmt();
void ParseRfmt();
void convertDtoE();
void readHB_header();
void readHB_mat_double();
void readHB_rhs_double();
void writeHB_mat_double();
char * substr();
char * substr_after();
char * substr_before();
char * substr_through();
void upcase();
int my_index();

void readHB_info(filename, M, N, nz, nrhs)
char *filename; int *M; int *N; int *nz; 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).                                                      */
/*  The optional verbose parameter, if set to 1,  triggers the output of    */
/*  more detailed information from the header, including title, etc.        */
/*                                                                          */
/*  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.                 **        */
/*                                                                          */
/****************************************************************************/
    int verbose = 0;
    FILE *in_file;
    char Title[73], Key[9], Type[4], Rhstype[4];
    char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
    int Ptrcrd, Indcrd, Valcrd, Rhscrd; 
    int Nrow, Ncol, Nnzero;
    int Nrhs;
    
    in_file = fopen( filename, "r");
    if (in_file == NULL)
    {
       printf("Cannot open file: %s\n",filename);
       exit(1);
    }

    readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
                  Ptrfmt, Indfmt, Valfmt, Rhsfmt, 
                  &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);

    *M    = Nrow;
    *N    = Ncol;
    *nz   = Nnzero;
    if (Rhscrd != 0) {*nrhs = Nrhs;}
                else {*nrhs = 0;}

/*  In verbose mode, print some of the header information:   */
    if (verbose == 1)
    {
        printf("Reading from Harwell-Boeing file %s...\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) are stored.\n",*nrhs);
    }
 
    return;

}


/*************************************************************************/
/*  Read header information from the named H/B file...                   */
/*************************************************************************/

void  readHB_header(in_file, Title, Key, Type, 
                    Nrow, Ncol, Nnzero, Nrhs,
                    Ptrfmt, Indfmt, Valfmt, Rhsfmt, 
                    Ptrcrd, Indcrd, Valcrd, Rhscrd, 
                    Rhstype)
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;
{
    char line[82];
    char *line_ptr;
    int Totcrd;
    int Neltvl, Nrhsix;

/*  First line:   */
    line_ptr = fgets(line, 82, in_file);
    nullchk(line_ptr);
    (void) sscanf(line, "%72c %8c", Title, Key);

/*  Second line:  */
    line_ptr = fgets(line, 82, in_file);
    nullchk(line_ptr);
    (void) sscanf(line, "%i %i %i %i %i", &Totcrd, Ptrcrd, Indcrd, 
                                          Valcrd, Rhscrd);

/*  Third line:   */
    line_ptr = fgets(line, 82, in_file);
    nullchk(line_ptr);
    (void) sscanf(line, "%3c %i %i %i %i", Type, Nrow, Ncol, 
                                          Nnzero, &Neltvl);

/*  Fourth line:  */
    line_ptr = fgets(line, 82, in_file);
    nullchk(line_ptr);
    (void) sscanf(line, "%s %s %s %s", Ptrfmt, Indfmt, Valfmt, Rhsfmt);
   
/*  (Optional) Fifth line: */
    if (*Rhscrd != 0 )
    { 
       line_ptr = fgets(line, 82, in_file);
       nullchk(line_ptr);
       (void) sscanf(line, "%3c %i %i", Rhstype, Nrhs, &Nrhsix);
    }

    return;
}



void readHB_mat_double(filename, colptr, rowind, val)
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;
    char Title[73], Key[9], Type[4], Rhstype[4];
    char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
    int Ptrcrd, Indcrd, Valcrd, Rhscrd;
    int Nrow, Ncol, Nnzero;
    int Nrhs;
    char line[82];
    char* line_ptr;
    int Ptrperline, Ptrwidth, Indperline, Indwidth;
    int Valperline, Valwidth;
    int Valflag;           /* Indicates 'E','D', or 'F' float format */
    int i, ind, col;
    int count;
    char* ThisElement;

    in_file = fopen( filename, "r");
    if (in_file == NULL)
    {
       printf("Cannot open file: %s\n",filename);
       exit(1);
    }

    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);
    ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valflag);

/*  Read column pointer array:   */

    count=0;
    for (i=0;i<Ptrcrd;i++)
    {
       line_ptr = fgets(line, 82, in_file);
       nullchk(line_ptr);
       col =  0;
       for (ind = 0;ind<Ptrperline;ind++)
       {
          if (count > Ncol) break;
          ThisElement = substr(line,col,Ptrwidth);
          colptr[count] = atoi(ThisElement)-1;
          count++; col += Ptrwidth;
       }
    }

/*  Read row index array:  */

    count = 0;
    for (i=0;i<Indcrd;i++)
    {
       line_ptr = fgets(line, 82, in_file);
       nullchk(line_ptr);
       col =  0;
       for (ind = 0;ind<Indperline;ind++)
       {
          if (count == Nnzero) break;
          ThisElement = substr(line,col,Indwidth);
          rowind[count] = atoi(ThisElement)-1;
          count++; col += Indwidth;
       }
    }

/*  Read array of values:  */

    count = 0;
    for (i=0;i<Valcrd;i++)
    {
       line_ptr = fgets(line, 82, in_file);
       nullchk(line_ptr);
       if (Valflag == 'D') convertDtoE(line);
       col =  0;
       for (ind = 0;ind<Valperline;ind++)
       {
          if (count == Nnzero) break;
          ThisElement = substr(line,col,Valwidth);
          val[count] = atof(ThisElement);
          count++; col += Valwidth;
       }
    }

    return;
}


void readHB_mat_float(filename, colptr, rowind, val)
char *filename; int *colptr; int *rowind; float *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;
    char Title[73], Key[9], Type[4], Rhstype[4];
    char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
    int Ptrcrd, Indcrd, Valcrd, Rhscrd;
    int Nrow, Ncol, Nnzero;
    int Nrhs;
    char line[82];
    char* line_ptr;
    int Ptrperline, Ptrwidth, Indperline, Indwidth;
    int Valperline, Valwidth;
    int Valflag;           /* Indicates 'E','D', or 'F' float format */
    int i, ind, col;
    int count;
    char* ThisElement;

    in_file = fopen( filename, "r");
    if (in_file == NULL)
    {
       printf("Cannot open file: %s\n",filename);
       exit(1);
    }

    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);
    ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valflag);

/*  Read column pointer array:   */

    count=0;
    for (i=0;i<Ptrcrd;i++)
    {
       line_ptr = fgets(line, 82, in_file);