src/glphbm.c
changeset 1 c445c931472f
equal deleted inserted replaced
-1:000000000000 0:9830d3ae7a85
       
     1 /* glphbm.c */
       
     2 
       
     3 /***********************************************************************
       
     4 *  This code is part of GLPK (GNU Linear Programming Kit).
       
     5 *
       
     6 *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
       
     7 *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
       
     8 *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
       
     9 *  E-mail: <mao@gnu.org>.
       
    10 *
       
    11 *  GLPK is free software: you can redistribute it and/or modify it
       
    12 *  under the terms of the GNU General Public License as published by
       
    13 *  the Free Software Foundation, either version 3 of the License, or
       
    14 *  (at your option) any later version.
       
    15 *
       
    16 *  GLPK is distributed in the hope that it will be useful, but WITHOUT
       
    17 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
       
    18 *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
       
    19 *  License for more details.
       
    20 *
       
    21 *  You should have received a copy of the GNU General Public License
       
    22 *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
       
    23 ***********************************************************************/
       
    24 
       
    25 #define _GLPSTD_ERRNO
       
    26 #define _GLPSTD_STDIO
       
    27 #include "glphbm.h"
       
    28 #include "glpenv.h"
       
    29 
       
    30 /***********************************************************************
       
    31 *  NAME
       
    32 *
       
    33 *  hbm_read_mat - read sparse matrix in Harwell-Boeing format
       
    34 *
       
    35 *  SYNOPSIS
       
    36 *
       
    37 *  #include "glphbm.h"
       
    38 *  HBM *hbm_read_mat(const char *fname);
       
    39 *
       
    40 *  DESCRIPTION
       
    41 *
       
    42 *  The routine hbm_read_mat reads a sparse matrix in the Harwell-Boeing
       
    43 *  format from a text file whose name is the character string fname.
       
    44 *
       
    45 *  Detailed description of the Harwell-Boeing format recognised by this
       
    46 *  routine is given in the following report:
       
    47 *
       
    48 *  I.S.Duff, R.G.Grimes, J.G.Lewis. User's Guide for the Harwell-Boeing
       
    49 *  Sparse Matrix Collection (Release I), TR/PA/92/86, October 1992.
       
    50 *
       
    51 *  RETURNS
       
    52 *
       
    53 *  If no error occured, the routine hbm_read_mat returns a pointer to
       
    54 *  a data structure containing the matrix. In case of error the routine
       
    55 *  prints an appropriate error message and returns NULL. */
       
    56 
       
    57 struct dsa
       
    58 {     /* working area used by routine hbm_read_mat */
       
    59       const char *fname;
       
    60       /* name of input text file */
       
    61       FILE *fp;
       
    62       /* stream assigned to input text file */
       
    63       int seqn;
       
    64       /* card sequential number */
       
    65       char card[80+1];
       
    66       /* card image buffer */
       
    67       int fmt_p;
       
    68       /* scale factor */
       
    69       int fmt_k;
       
    70       /* iterator */
       
    71       int fmt_f;
       
    72       /* format code */
       
    73       int fmt_w;
       
    74       /* field width */
       
    75       int fmt_d;
       
    76       /* number of decimal places after point */
       
    77 };
       
    78 
       
    79 /***********************************************************************
       
    80 *  read_card - read next data card
       
    81 *
       
    82 *  This routine reads the next 80-column card from the input text file
       
    83 *  and stores its image into the character string card. If the card was
       
    84 *  read successfully, the routine returns zero, otherwise non-zero. */
       
    85 
       
    86 static int read_card(struct dsa *dsa)
       
    87 {     int k, c;
       
    88       dsa->seqn++;
       
    89       memset(dsa->card, ' ', 80), dsa->card[80] = '\0';
       
    90       k = 0;
       
    91       for (;;)
       
    92       {  c = fgetc(dsa->fp);
       
    93          if (ferror(dsa->fp))
       
    94          {  xprintf("%s:%d: read error - %s\n", dsa->fname, dsa->seqn,
       
    95                strerror(errno));
       
    96             return 1;
       
    97          }
       
    98          if (feof(dsa->fp))
       
    99          {  if (k == 0)
       
   100                xprintf("%s:%d: unexpected EOF\n", dsa->fname,
       
   101                   dsa->seqn);
       
   102             else
       
   103                xprintf("%s:%d: missing final LF\n", dsa->fname,
       
   104                   dsa->seqn);
       
   105             return 1;
       
   106          }
       
   107          if (c == '\r') continue;
       
   108          if (c == '\n') break;
       
   109          if (iscntrl(c))
       
   110          {  xprintf("%s:%d: invalid control character 0x%02X\n",
       
   111                dsa->fname, dsa->seqn, c);
       
   112             return 1;
       
   113          }
       
   114          if (k == 80)
       
   115          {  xprintf("%s:%d: card image too long\n", dsa->fname,
       
   116                dsa->seqn);
       
   117             return 1;
       
   118          }
       
   119          dsa->card[k++] = (char)c;
       
   120       }
       
   121       return 0;
       
   122 }
       
   123 
       
   124 /***********************************************************************
       
   125 *  scan_int - scan integer value from the current card
       
   126 *
       
   127 *  This routine scans an integer value from the current card, where fld
       
   128 *  is the name of the field, pos is the position of the field, width is
       
   129 *  the width of the field, val points to a location to which the scanned
       
   130 *  value should be stored. If the value was scanned successfully, the
       
   131 *  routine returns zero, otherwise non-zero. */
       
   132 
       
   133 static int scan_int(struct dsa *dsa, char *fld, int pos, int width,
       
   134       int *val)
       
   135 {     char str[80+1];
       
   136       xassert(1 <= width && width <= 80);
       
   137       memcpy(str, dsa->card + pos, width), str[width] = '\0';
       
   138       if (str2int(strspx(str), val))
       
   139       {  xprintf("%s:%d: field `%s' contains invalid value `%s'\n",
       
   140             dsa->fname, dsa->seqn, fld, str);
       
   141          return 1;
       
   142       }
       
   143       return 0;
       
   144 }
       
   145 
       
   146 /***********************************************************************
       
   147 *  parse_fmt - parse Fortran format specification
       
   148 *
       
   149 *  This routine parses the Fortran format specification represented as
       
   150 *  character string which fmt points to and stores format elements into
       
   151 *  appropriate static locations. Should note that not all valid Fortran
       
   152 *  format specifications may be recognised. If the format specification
       
   153 *  was recognised, the routine returns zero, otherwise non-zero. */
       
   154 
       
   155 static int parse_fmt(struct dsa *dsa, char *fmt)
       
   156 {     int k, s, val;
       
   157       char str[80+1];
       
   158       /* first character should be left parenthesis */
       
   159       if (fmt[0] != '(')
       
   160 fail: {  xprintf("hbm_read_mat: format `%s' not recognised\n", fmt);
       
   161          return 1;
       
   162       }
       
   163       k = 1;
       
   164       /* optional scale factor */
       
   165       dsa->fmt_p = 0;
       
   166       if (isdigit((unsigned char)fmt[k]))
       
   167       {  s = 0;
       
   168          while (isdigit((unsigned char)fmt[k]))
       
   169          {  if (s == 80) goto fail;
       
   170             str[s++] = fmt[k++];
       
   171          }
       
   172          str[s] = '\0';
       
   173          if (str2int(str, &val)) goto fail;
       
   174          if (toupper((unsigned char)fmt[k]) != 'P') goto iter;
       
   175          dsa->fmt_p = val, k++;
       
   176          if (!(0 <= dsa->fmt_p && dsa->fmt_p <= 255)) goto fail;
       
   177          /* optional comma may follow scale factor */
       
   178          if (fmt[k] == ',') k++;
       
   179       }
       
   180       /* optional iterator */
       
   181       dsa->fmt_k = 1;
       
   182       if (isdigit((unsigned char)fmt[k]))
       
   183       {  s = 0;
       
   184          while (isdigit((unsigned char)fmt[k]))
       
   185          {  if (s == 80) goto fail;
       
   186             str[s++] = fmt[k++];
       
   187          }
       
   188          str[s] = '\0';
       
   189          if (str2int(str, &val)) goto fail;
       
   190 iter:    dsa->fmt_k = val;
       
   191          if (!(1 <= dsa->fmt_k && dsa->fmt_k <= 255)) goto fail;
       
   192       }
       
   193       /* format code */
       
   194       dsa->fmt_f = toupper((unsigned char)fmt[k++]);
       
   195       if (!(dsa->fmt_f == 'D' || dsa->fmt_f == 'E' ||
       
   196             dsa->fmt_f == 'F' || dsa->fmt_f == 'G' ||
       
   197             dsa->fmt_f == 'I')) goto fail;
       
   198       /* field width */
       
   199       if (!isdigit((unsigned char)fmt[k])) goto fail;
       
   200       s = 0;
       
   201       while (isdigit((unsigned char)fmt[k]))
       
   202       {  if (s == 80) goto fail;
       
   203          str[s++] = fmt[k++];
       
   204       }
       
   205       str[s] = '\0';
       
   206       if (str2int(str, &dsa->fmt_w)) goto fail;
       
   207       if (!(1 <= dsa->fmt_w && dsa->fmt_w <= 255)) goto fail;
       
   208       /* optional number of decimal places after point */
       
   209       dsa->fmt_d = 0;
       
   210       if (fmt[k] == '.')
       
   211       {  k++;
       
   212          if (!isdigit((unsigned char)fmt[k])) goto fail;
       
   213          s = 0;
       
   214          while (isdigit((unsigned char)fmt[k]))
       
   215          {  if (s == 80) goto fail;
       
   216             str[s++] = fmt[k++];
       
   217          }
       
   218          str[s] = '\0';
       
   219          if (str2int(str, &dsa->fmt_d)) goto fail;
       
   220          if (!(0 <= dsa->fmt_d && dsa->fmt_d <= 255)) goto fail;
       
   221       }
       
   222       /* last character should be right parenthesis */
       
   223       if (!(fmt[k] == ')' && fmt[k+1] == '\0')) goto fail;
       
   224       return 0;
       
   225 }
       
   226 
       
   227 /***********************************************************************
       
   228 *  read_int_array - read array of integer type
       
   229 *
       
   230 *  This routine reads an integer array from the input text file, where
       
   231 *  name is array name, fmt is Fortran format specification that controls
       
   232 *  reading, n is number of array elements, val is array of integer type.
       
   233 *  If the array was read successful, the routine returns zero, otherwise
       
   234 *  non-zero. */
       
   235 
       
   236 static int read_int_array(struct dsa *dsa, char *name, char *fmt,
       
   237       int n, int val[])
       
   238 {     int k, pos;
       
   239       char str[80+1];
       
   240       if (parse_fmt(dsa, fmt)) return 1;
       
   241       if (!(dsa->fmt_f == 'I' && dsa->fmt_w <= 80 &&
       
   242             dsa->fmt_k * dsa->fmt_w <= 80))
       
   243       {  xprintf(
       
   244             "%s:%d: can't read array `%s' - invalid format `%s'\n",
       
   245             dsa->fname, dsa->seqn, name, fmt);
       
   246          return 1;
       
   247       }
       
   248       for (k = 1, pos = INT_MAX; k <= n; k++, pos++)
       
   249       {  if (pos >= dsa->fmt_k)
       
   250          {  if (read_card(dsa)) return 1;
       
   251             pos = 0;
       
   252          }
       
   253          memcpy(str, dsa->card + dsa->fmt_w * pos, dsa->fmt_w);
       
   254          str[dsa->fmt_w] = '\0';
       
   255          strspx(str);
       
   256          if (str2int(str, &val[k]))
       
   257          {  xprintf(
       
   258                "%s:%d: can't read array `%s' - invalid value `%s'\n",
       
   259                dsa->fname, dsa->seqn, name, str);
       
   260             return 1;
       
   261          }
       
   262       }
       
   263       return 0;
       
   264 }
       
   265 
       
   266 /***********************************************************************
       
   267 *  read_real_array - read array of real type
       
   268 *
       
   269 *  This routine reads a real array from the input text file, where name
       
   270 *  is array name, fmt is Fortran format specification that controls
       
   271 *  reading, n is number of array elements, val is array of real type.
       
   272 *  If the array was read successful, the routine returns zero, otherwise
       
   273 *  non-zero. */
       
   274 
       
   275 static int read_real_array(struct dsa *dsa, char *name, char *fmt,
       
   276       int n, double val[])
       
   277 {     int k, pos;
       
   278       char str[80+1], *ptr;
       
   279       if (parse_fmt(dsa, fmt)) return 1;
       
   280       if (!(dsa->fmt_f != 'I' && dsa->fmt_w <= 80 &&
       
   281             dsa->fmt_k * dsa->fmt_w <= 80))
       
   282       {  xprintf(
       
   283             "%s:%d: can't read array `%s' - invalid format `%s'\n",
       
   284             dsa->fname, dsa->seqn, name, fmt);
       
   285          return 1;
       
   286       }
       
   287       for (k = 1, pos = INT_MAX; k <= n; k++, pos++)
       
   288       {  if (pos >= dsa->fmt_k)
       
   289          {  if (read_card(dsa)) return 1;
       
   290             pos = 0;
       
   291          }
       
   292          memcpy(str, dsa->card + dsa->fmt_w * pos, dsa->fmt_w);
       
   293          str[dsa->fmt_w] = '\0';
       
   294          strspx(str);
       
   295          if (strchr(str, '.') == NULL && strcmp(str, "0"))
       
   296          {  xprintf("%s(%d): can't read array `%s' - value `%s' has no "
       
   297                "decimal point\n", dsa->fname, dsa->seqn, name, str);
       
   298             return 1;
       
   299          }
       
   300          /* sometimes lower case letters appear */
       
   301          for (ptr = str; *ptr; ptr++)
       
   302             *ptr = (char)toupper((unsigned char)*ptr);
       
   303          ptr = strchr(str, 'D');
       
   304          if (ptr != NULL) *ptr = 'E';
       
   305          /* value may appear with decimal exponent but without letters
       
   306             E or D (for example, -123.456-012), so missing letter should
       
   307             be inserted */
       
   308          ptr = strchr(str+1, '+');
       
   309          if (ptr == NULL) ptr = strchr(str+1, '-');
       
   310          if (ptr != NULL && *(ptr-1) != 'E')
       
   311          {  xassert(strlen(str) < 80);
       
   312             memmove(ptr+1, ptr, strlen(ptr)+1);
       
   313             *ptr = 'E';
       
   314          }
       
   315          if (str2num(str, &val[k]))
       
   316          {  xprintf(
       
   317                "%s:%d: can't read array `%s' - invalid value `%s'\n",
       
   318                dsa->fname, dsa->seqn, name, str);
       
   319             return 1;
       
   320          }
       
   321       }
       
   322       return 0;
       
   323 }
       
   324 
       
   325 HBM *hbm_read_mat(const char *fname)
       
   326 {     struct dsa _dsa, *dsa = &_dsa;
       
   327       HBM *hbm = NULL;
       
   328       dsa->fname = fname;
       
   329       xprintf("hbm_read_mat: reading matrix from `%s'...\n",
       
   330          dsa->fname);
       
   331       dsa->fp = fopen(dsa->fname, "r");
       
   332       if (dsa->fp == NULL)
       
   333       {  xprintf("hbm_read_mat: unable to open `%s' - %s\n",
       
   334             dsa->fname, strerror(errno));
       
   335          goto fail;
       
   336       }
       
   337       dsa->seqn = 0;
       
   338       hbm = xmalloc(sizeof(HBM));
       
   339       memset(hbm, 0, sizeof(HBM));
       
   340       /* read the first heading card */
       
   341       if (read_card(dsa)) goto fail;
       
   342       memcpy(hbm->title, dsa->card, 72), hbm->title[72] = '\0';
       
   343       strtrim(hbm->title);
       
   344       xprintf("%s\n", hbm->title);
       
   345       memcpy(hbm->key, dsa->card+72, 8), hbm->key[8] = '\0';
       
   346       strspx(hbm->key);
       
   347       xprintf("key = %s\n", hbm->key);
       
   348       /* read the second heading card */
       
   349       if (read_card(dsa)) goto fail;
       
   350       if (scan_int(dsa, "totcrd",  0, 14, &hbm->totcrd)) goto fail;
       
   351       if (scan_int(dsa, "ptrcrd", 14, 14, &hbm->ptrcrd)) goto fail;
       
   352       if (scan_int(dsa, "indcrd", 28, 14, &hbm->indcrd)) goto fail;
       
   353       if (scan_int(dsa, "valcrd", 42, 14, &hbm->valcrd)) goto fail;
       
   354       if (scan_int(dsa, "rhscrd", 56, 14, &hbm->rhscrd)) goto fail;
       
   355       xprintf("totcrd = %d; ptrcrd = %d; indcrd = %d; valcrd = %d; rhsc"
       
   356          "rd = %d\n", hbm->totcrd, hbm->ptrcrd, hbm->indcrd,
       
   357          hbm->valcrd, hbm->rhscrd);
       
   358       /* read the third heading card */
       
   359       if (read_card(dsa)) goto fail;
       
   360       memcpy(hbm->mxtype, dsa->card, 3), hbm->mxtype[3] = '\0';
       
   361       if (strchr("RCP",   hbm->mxtype[0]) == NULL ||
       
   362           strchr("SUHZR", hbm->mxtype[1]) == NULL ||
       
   363           strchr("AE",    hbm->mxtype[2]) == NULL)
       
   364       {  xprintf("%s:%d: matrix type `%s' not recognised\n",
       
   365             dsa->fname, dsa->seqn, hbm->mxtype);
       
   366          goto fail;
       
   367       }
       
   368       if (scan_int(dsa, "nrow", 14, 14, &hbm->nrow)) goto fail;
       
   369       if (scan_int(dsa, "ncol", 28, 14, &hbm->ncol)) goto fail;
       
   370       if (scan_int(dsa, "nnzero", 42, 14, &hbm->nnzero)) goto fail;
       
   371       if (scan_int(dsa, "neltvl", 56, 14, &hbm->neltvl)) goto fail;
       
   372       xprintf("mxtype = %s; nrow = %d; ncol = %d; nnzero = %d; neltvl ="
       
   373          " %d\n", hbm->mxtype, hbm->nrow, hbm->ncol, hbm->nnzero,
       
   374          hbm->neltvl);
       
   375       /* read the fourth heading card */
       
   376       if (read_card(dsa)) goto fail;
       
   377       memcpy(hbm->ptrfmt, dsa->card, 16), hbm->ptrfmt[16] = '\0';
       
   378       strspx(hbm->ptrfmt);
       
   379       memcpy(hbm->indfmt, dsa->card+16, 16), hbm->indfmt[16] = '\0';
       
   380       strspx(hbm->indfmt);
       
   381       memcpy(hbm->valfmt, dsa->card+32, 20), hbm->valfmt[20] = '\0';
       
   382       strspx(hbm->valfmt);
       
   383       memcpy(hbm->rhsfmt, dsa->card+52, 20), hbm->rhsfmt[20] = '\0';
       
   384       strspx(hbm->rhsfmt);
       
   385       xprintf("ptrfmt = %s; indfmt = %s; valfmt = %s; rhsfmt = %s\n",
       
   386          hbm->ptrfmt, hbm->indfmt, hbm->valfmt, hbm->rhsfmt);
       
   387       /* read the fifth heading card (optional) */
       
   388       if (hbm->rhscrd <= 0)
       
   389       {  strcpy(hbm->rhstyp, "???");
       
   390          hbm->nrhs = 0;
       
   391          hbm->nrhsix = 0;
       
   392       }
       
   393       else
       
   394       {  if (read_card(dsa)) goto fail;
       
   395          memcpy(hbm->rhstyp, dsa->card, 3), hbm->rhstyp[3] = '\0';
       
   396          if (scan_int(dsa, "nrhs", 14, 14, &hbm->nrhs)) goto fail;
       
   397          if (scan_int(dsa, "nrhsix", 28, 14, &hbm->nrhsix)) goto fail;
       
   398          xprintf("rhstyp = `%s'; nrhs = %d; nrhsix = %d\n",
       
   399             hbm->rhstyp, hbm->nrhs, hbm->nrhsix);
       
   400       }
       
   401       /* read matrix structure */
       
   402       hbm->colptr = xcalloc(1+hbm->ncol+1, sizeof(int));
       
   403       if (read_int_array(dsa, "colptr", hbm->ptrfmt, hbm->ncol+1,
       
   404          hbm->colptr)) goto fail;
       
   405       hbm->rowind = xcalloc(1+hbm->nnzero, sizeof(int));
       
   406       if (read_int_array(dsa, "rowind", hbm->indfmt, hbm->nnzero,
       
   407          hbm->rowind)) goto fail;
       
   408       /* read matrix values */
       
   409       if (hbm->valcrd <= 0) goto done;
       
   410       if (hbm->mxtype[2] == 'A')
       
   411       {  /* assembled matrix */
       
   412          hbm->values = xcalloc(1+hbm->nnzero, sizeof(double));
       
   413          if (read_real_array(dsa, "values", hbm->valfmt, hbm->nnzero,
       
   414             hbm->values)) goto fail;
       
   415       }
       
   416       else
       
   417       {  /* elemental (unassembled) matrix */
       
   418          hbm->values = xcalloc(1+hbm->neltvl, sizeof(double));
       
   419          if (read_real_array(dsa, "values", hbm->valfmt, hbm->neltvl,
       
   420             hbm->values)) goto fail;
       
   421       }
       
   422       /* read right-hand sides */
       
   423       if (hbm->nrhs <= 0) goto done;
       
   424       if (hbm->rhstyp[0] == 'F')
       
   425       {  /* dense format */
       
   426          hbm->nrhsvl = hbm->nrow * hbm->nrhs;
       
   427          hbm->rhsval = xcalloc(1+hbm->nrhsvl, sizeof(double));
       
   428          if (read_real_array(dsa, "rhsval", hbm->rhsfmt, hbm->nrhsvl,
       
   429             hbm->rhsval)) goto fail;
       
   430       }
       
   431       else if (hbm->rhstyp[0] == 'M' && hbm->mxtype[2] == 'A')
       
   432       {  /* sparse format */
       
   433          /* read pointers */
       
   434          hbm->rhsptr = xcalloc(1+hbm->nrhs+1, sizeof(int));
       
   435          if (read_int_array(dsa, "rhsptr", hbm->ptrfmt, hbm->nrhs+1,
       
   436             hbm->rhsptr)) goto fail;
       
   437          /* read sparsity pattern */
       
   438          hbm->rhsind = xcalloc(1+hbm->nrhsix, sizeof(int));
       
   439          if (read_int_array(dsa, "rhsind", hbm->indfmt, hbm->nrhsix,
       
   440             hbm->rhsind)) goto fail;
       
   441          /* read values */
       
   442          hbm->rhsval = xcalloc(1+hbm->nrhsix, sizeof(double));
       
   443          if (read_real_array(dsa, "rhsval", hbm->rhsfmt, hbm->nrhsix,
       
   444             hbm->rhsval)) goto fail;
       
   445       }
       
   446       else if (hbm->rhstyp[0] == 'M' && hbm->mxtype[2] == 'E')
       
   447       {  /* elemental format */
       
   448          hbm->rhsval = xcalloc(1+hbm->nrhsvl, sizeof(double));
       
   449          if (read_real_array(dsa, "rhsval", hbm->rhsfmt, hbm->nrhsvl,
       
   450             hbm->rhsval)) goto fail;
       
   451       }
       
   452       else
       
   453       {  xprintf("%s:%d: right-hand side type `%c' not recognised\n",
       
   454             dsa->fname, dsa->seqn, hbm->rhstyp[0]);
       
   455          goto fail;
       
   456       }
       
   457       /* read starting guesses */
       
   458       if (hbm->rhstyp[1] == 'G')
       
   459       {  hbm->nguess = hbm->nrow * hbm->nrhs;
       
   460          hbm->sguess = xcalloc(1+hbm->nguess, sizeof(double));
       
   461          if (read_real_array(dsa, "sguess", hbm->rhsfmt, hbm->nguess,
       
   462             hbm->sguess)) goto fail;
       
   463       }
       
   464       /* read solution vectors */
       
   465       if (hbm->rhstyp[2] == 'X')
       
   466       {  hbm->nexact = hbm->nrow * hbm->nrhs;
       
   467          hbm->xexact = xcalloc(1+hbm->nexact, sizeof(double));
       
   468          if (read_real_array(dsa, "xexact", hbm->rhsfmt, hbm->nexact,
       
   469             hbm->xexact)) goto fail;
       
   470       }
       
   471 done: /* reading has been completed */
       
   472       xprintf("hbm_read_mat: %d cards were read\n", dsa->seqn);
       
   473       fclose(dsa->fp);
       
   474       return hbm;
       
   475 fail: /* something wrong in Danish kingdom */
       
   476       if (hbm != NULL)
       
   477       {  if (hbm->colptr != NULL) xfree(hbm->colptr);
       
   478          if (hbm->rowind != NULL) xfree(hbm->rowind);
       
   479          if (hbm->rhsptr != NULL) xfree(hbm->rhsptr);
       
   480          if (hbm->rhsind != NULL) xfree(hbm->rhsind);
       
   481          if (hbm->values != NULL) xfree(hbm->values);
       
   482          if (hbm->rhsval != NULL) xfree(hbm->rhsval);
       
   483          if (hbm->sguess != NULL) xfree(hbm->sguess);
       
   484          if (hbm->xexact != NULL) xfree(hbm->xexact);
       
   485          xfree(hbm);
       
   486       }
       
   487       if (dsa->fp != NULL) fclose(dsa->fp);
       
   488       return NULL;
       
   489 }
       
   490 
       
   491 /***********************************************************************
       
   492 *  NAME
       
   493 *
       
   494 *  hbm_free_mat - free sparse matrix in Harwell-Boeing format
       
   495 *
       
   496 *  SYNOPSIS
       
   497 *
       
   498 *  #include "glphbm.h"
       
   499 *  void hbm_free_mat(HBM *hbm);
       
   500 *
       
   501 *  DESCRIPTION
       
   502 *
       
   503 *  The hbm_free_mat routine frees all the memory allocated to the data
       
   504 *  structure containing a sparse matrix in the Harwell-Boeing format. */
       
   505 
       
   506 void hbm_free_mat(HBM *hbm)
       
   507 {     if (hbm->colptr != NULL) xfree(hbm->colptr);
       
   508       if (hbm->rowind != NULL) xfree(hbm->rowind);
       
   509       if (hbm->rhsptr != NULL) xfree(hbm->rhsptr);
       
   510       if (hbm->rhsind != NULL) xfree(hbm->rhsind);
       
   511       if (hbm->values != NULL) xfree(hbm->values);
       
   512       if (hbm->rhsval != NULL) xfree(hbm->rhsval);
       
   513       if (hbm->sguess != NULL) xfree(hbm->sguess);
       
   514       if (hbm->xexact != NULL) xfree(hbm->xexact);
       
   515       xfree(hbm);
       
   516       return;
       
   517 }
       
   518 
       
   519 /* eof */