iohb.c

很好用的库

    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, 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(char* message) 
{
   fprintf(stderr,message);
   exit(1);
}