iohb.c

来自「Boost provides free peer-reviewed portab」· C语言 代码 · 共 1,611 行 · 第 1/5 页

    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 */