iohb.c

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

       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_rhs_double(filename, b, j)
char *filename; double b[]; int j;
{
/****************************************************************************/
/*  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.                          */
/*                                                                          */
/*    ----------                                                            */
/*    **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.                 **        */
/*                                                                          */
/****************************************************************************/
    int i, n, null_entries, lines_left;
    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;
    int Rhsperline, Rhswidth;
    int Rhsflag;
    char buffer[BUFSIZ];
    int ind, col;
    int count;
    char *line_ptr;
    char line[80];
    char *ThisElement;

    if ((in_file = fopen( filename, "r")) == 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);

    if (Nrhs <= 0)
    {
      printf("Attempt to read rhs when none is present.\n");
      exit(1);
      return;
    }
    if (Rhstype[0] != 'F')
    {
      printf("Attempt to read rhs which is not stored in Full form.\n");
      printf("Rhs must be specified as full. \n");
      exit(1);
      return;
    }

    ParseRfmt(Rhsfmt, &Rhsperline, &Rhswidth, &Rhsflag);

/*  Lines to skip before starting to read RHS values... */
    n = Ptrcrd + Indcrd + Valcrd + j*Nrow/Rhsperline;

/*  number of entries on the line of interest to skip before */
/*  reading RHS values...                                    */
    null_entries = j*Nrow%Rhsperline;
    lines_left = (int) ( .5 + (double) (Nrow-Rhsperline+null_entries)/ 
                                            (double) Rhsperline );

    for (i = 0; i < n; i++)
      fgets(buffer, BUFSIZ, in_file);

    count = 0;

/*  Handle first line separately, in case j != 0 and the right-hand-side */
/*  of interest starts mid-line (after null_entries entries)...          */

    line_ptr = fgets(line, 82, in_file);
    nullchk(line_ptr);
    col = 0;
    for (ind = 0; ind < Rhsperline; ind++)
    {
      if (ind < null_entries)
        col += Rhswidth;
      else
      {
        if (count > Nrow-1)
          break;
        ThisElement = substr(line, col, Rhswidth);
        b[count] = atof(ThisElement);
        count++; col += Rhswidth;
      }
    }

/*  Handle subsequent lines...              */
    for (i = 0; i < lines_left; i++)
    {
      line_ptr = fgets(line, 82, in_file);
      nullchk(line_ptr);
      col = 0;
      for (ind = 0; ind < Rhsperline; ind++)
      {
         if (count > Nrow-1)
           break;
         ThisElement = substr(line, col, Rhswidth);
         b[count] = atof(ThisElement);
         count++; col += Rhswidth;
      }
    }

    return;
}

void readHB_rhs_float(filename, b, j)
char *filename; float b[]; int j;
{
/****************************************************************************/
/*  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.                          */
/*                                                                          */
/*    ----------                                                            */
/*    **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.                 **        */
/*                                                                          */
/****************************************************************************/
    int i, n, null_entries, lines_left;
    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;
    int Rhsperline, Rhswidth;
    int Rhsflag;
    char buffer[BUFSIZ];
    int ind, col;
    int count;
    char *line_ptr;
    char line[80];
    char *ThisElement;

    if ((in_file = fopen( filename, "r")) == 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);

    if (Nrhs <= 0)
    {
      printf("Attempt to read rhs when none is present.\n");
      exit(1);
      return;
    }
    if (Rhstype[0] != 'F')
    {
      printf("Attempt to read rhs which is not stored in Full form.\n");
      printf("Rhs must be specified as full. \n");
      exit(1);
      return;
    }

    ParseRfmt(Rhsfmt, &Rhsperline, &Rhswidth, &Rhsflag);

/*  Lines to skip before starting to read RHS values... */
    n = Ptrcrd + Indcrd + Valcrd + j*Nrow/Rhsperline;

/*  number of entries on the line of interest to skip before */
/*  reading RHS values...                                    */
    null_entries = j*Nrow%Rhsperline;
    lines_left = (int)( .5 + (double) (Nrow-Rhsperline+null_entries)/
                                             (double) Rhsperline);

    for (i = 0; i < n; i++)
      fgets(buffer, BUFSIZ, in_file);

    count = 0;

/*  Handle first line separately, in case j != 0 and the right-hand-side */
/*  of interest starts mid-line (after null_entries entries)...          */

    line_ptr = fgets(line, 82, in_file);
    nullchk(line_ptr);
    col = 0;
    for (ind = 0; ind < Rhsperline; ind++)
    {
      if (ind < null_entries)
        col += Rhswidth;
      else
      {
        if (count > Nrow-1)
          break;
        ThisElement = substr(line, col, Rhswidth);
        b[count] = atof(ThisElement);
        count++; col += Rhswidth;
      }
    }

/*  Handle subsequent lines...              */
    for (i = 0; i < lines_left; i++)
    {
      line_ptr = fgets(line, 82, in_file);
      nullchk(line_ptr);
      col = 0;
      for (ind = 0; ind < Rhsperline; ind++)
      {
         if (count > Nrow-1)
           break;
         ThisElement = substr(line, col, Rhswidth);
         b[count] = atof(ThisElement);
         count++; col += Rhswidth;
      }
    }

    return;
}


void writeHB_mat_double(filename, M, N, nz, colptr, rowind, val, nrhs, rhs, 
                                                                  Title, Key)
char *filename; int M, N, nz, colptr[], rowind[]; double val[]; int nrhs;
double rhs[]; char *Title; char *Key;
{
/****************************************************************************/
/*  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 totcrd, ptrcrd, indcrd, valcrd, rhscrd;
    char *ptrfmt, *indfmt, *valfmt, *rhsfmt;
    char *Type;
    char *Titlefill, *Keyfill;
    int filllen, i;
    int entry, finfo;

    out_file = fopen( filename, "w");

    ptrcrd = (N+1)/8;
    if ( (N+1)%8 != 0) ptrcrd++;

    indcrd = nz/8;
    if ( nz%8 != 0) indcrd++;

    valcrd = nz/4;
    if ( nz%4 != 0) valcrd++;

    rhscrd = nrhs*M/4; 
    if ( nrhs*M%4 != 0) rhscrd++;

    totcrd = 4+ptrcrd+indcrd+valcrd+rhscrd;

    ptrfmt = "(8I10)          ";
    indfmt = ptrfmt;
    valfmt = "(4E20.16)           ";
    rhsfmt = valfmt;

    Type = "RUA";

/*  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%20s\n", ptrfmt, indfmt, valfmt, rhsfmt);
    if ( nrhs != 0 ) {
/*     Print optional fifth header line for right-hand-side information : */
       fprintf(out_file,"F             %d\n", nrhs);
    }

/*  Print column pointers:   */
    for (i=0;i<N+1;i++)
    {
       entry = colptr[i]+1;
       fprintf(out_file,"%10d",entry);
       if ( (i+1)%8 == 0 ) fprintf(out_file,"\n");
    }

   if ( (N+1) % 8 != 0 ) fprintf(out_file,"\n");

/*  Print row indices:       */
    for (i=0;i<nz;i++)
    {
       entry = rowind[i]+1;
       fprintf(out_file,"%10d",entry);
       if ( (i+1)%8 == 0 ) fprintf(out_file,"\n");
    }

   if ( nz % 8 != 0 ) fprintf(out_file,"\n");

/*  Print values:            */
    for (i=0;i<nz;i++)
    {
       fprintf(out_file,"% 20.12E",val[i]);
       if ( (i+1)%4 == 0 ) fprintf(out_file,"\n");
    }

    if ( nz % 4 != 0 ) fprintf(out_file,"\n");

/*  Print right hand sides:  */
    if ( nrhs > 0 ) {
       for (i=0;i<nrhs*M;i++)
       {
          fprintf(out_file,"% 20.12E",rhs[i]);
          if ( (i+1)%4 == 0 ) fprintf(out_file,"\n");
       }
    }

    finfo = fclose(out_file);
    if (finfo != 0) printf("Error closing file in writeHB_mat_double().\n");

    return;
}


void writeHB_mat_float(filename, M, N, nz, colptr, rowind, val, nrhs, rhs, 
                                                                  Title, Key)
char *filename; int M, N, nz, colptr[], rowind[]; float val[]; int nrhs;
float rhs[]; char *Title; char *Key;
{
/****************************************************************************/
/*  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              */
/*                                                                          */