1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glphbm.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,519 @@
     1.4 +/* glphbm.c */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#define _GLPSTD_ERRNO
    1.29 +#define _GLPSTD_STDIO
    1.30 +#include "glphbm.h"
    1.31 +#include "glpenv.h"
    1.32 +
    1.33 +/***********************************************************************
    1.34 +*  NAME
    1.35 +*
    1.36 +*  hbm_read_mat - read sparse matrix in Harwell-Boeing format
    1.37 +*
    1.38 +*  SYNOPSIS
    1.39 +*
    1.40 +*  #include "glphbm.h"
    1.41 +*  HBM *hbm_read_mat(const char *fname);
    1.42 +*
    1.43 +*  DESCRIPTION
    1.44 +*
    1.45 +*  The routine hbm_read_mat reads a sparse matrix in the Harwell-Boeing
    1.46 +*  format from a text file whose name is the character string fname.
    1.47 +*
    1.48 +*  Detailed description of the Harwell-Boeing format recognised by this
    1.49 +*  routine is given in the following report:
    1.50 +*
    1.51 +*  I.S.Duff, R.G.Grimes, J.G.Lewis. User's Guide for the Harwell-Boeing
    1.52 +*  Sparse Matrix Collection (Release I), TR/PA/92/86, October 1992.
    1.53 +*
    1.54 +*  RETURNS
    1.55 +*
    1.56 +*  If no error occured, the routine hbm_read_mat returns a pointer to
    1.57 +*  a data structure containing the matrix. In case of error the routine
    1.58 +*  prints an appropriate error message and returns NULL. */
    1.59 +
    1.60 +struct dsa
    1.61 +{     /* working area used by routine hbm_read_mat */
    1.62 +      const char *fname;
    1.63 +      /* name of input text file */
    1.64 +      FILE *fp;
    1.65 +      /* stream assigned to input text file */
    1.66 +      int seqn;
    1.67 +      /* card sequential number */
    1.68 +      char card[80+1];
    1.69 +      /* card image buffer */
    1.70 +      int fmt_p;
    1.71 +      /* scale factor */
    1.72 +      int fmt_k;
    1.73 +      /* iterator */
    1.74 +      int fmt_f;
    1.75 +      /* format code */
    1.76 +      int fmt_w;
    1.77 +      /* field width */
    1.78 +      int fmt_d;
    1.79 +      /* number of decimal places after point */
    1.80 +};
    1.81 +
    1.82 +/***********************************************************************
    1.83 +*  read_card - read next data card
    1.84 +*
    1.85 +*  This routine reads the next 80-column card from the input text file
    1.86 +*  and stores its image into the character string card. If the card was
    1.87 +*  read successfully, the routine returns zero, otherwise non-zero. */
    1.88 +
    1.89 +static int read_card(struct dsa *dsa)
    1.90 +{     int k, c;
    1.91 +      dsa->seqn++;
    1.92 +      memset(dsa->card, ' ', 80), dsa->card[80] = '\0';
    1.93 +      k = 0;
    1.94 +      for (;;)
    1.95 +      {  c = fgetc(dsa->fp);
    1.96 +         if (ferror(dsa->fp))
    1.97 +         {  xprintf("%s:%d: read error - %s\n", dsa->fname, dsa->seqn,
    1.98 +               strerror(errno));
    1.99 +            return 1;
   1.100 +         }
   1.101 +         if (feof(dsa->fp))
   1.102 +         {  if (k == 0)
   1.103 +               xprintf("%s:%d: unexpected EOF\n", dsa->fname,
   1.104 +                  dsa->seqn);
   1.105 +            else
   1.106 +               xprintf("%s:%d: missing final LF\n", dsa->fname,
   1.107 +                  dsa->seqn);
   1.108 +            return 1;
   1.109 +         }
   1.110 +         if (c == '\r') continue;
   1.111 +         if (c == '\n') break;
   1.112 +         if (iscntrl(c))
   1.113 +         {  xprintf("%s:%d: invalid control character 0x%02X\n",
   1.114 +               dsa->fname, dsa->seqn, c);
   1.115 +            return 1;
   1.116 +         }
   1.117 +         if (k == 80)
   1.118 +         {  xprintf("%s:%d: card image too long\n", dsa->fname,
   1.119 +               dsa->seqn);
   1.120 +            return 1;
   1.121 +         }
   1.122 +         dsa->card[k++] = (char)c;
   1.123 +      }
   1.124 +      return 0;
   1.125 +}
   1.126 +
   1.127 +/***********************************************************************
   1.128 +*  scan_int - scan integer value from the current card
   1.129 +*
   1.130 +*  This routine scans an integer value from the current card, where fld
   1.131 +*  is the name of the field, pos is the position of the field, width is
   1.132 +*  the width of the field, val points to a location to which the scanned
   1.133 +*  value should be stored. If the value was scanned successfully, the
   1.134 +*  routine returns zero, otherwise non-zero. */
   1.135 +
   1.136 +static int scan_int(struct dsa *dsa, char *fld, int pos, int width,
   1.137 +      int *val)
   1.138 +{     char str[80+1];
   1.139 +      xassert(1 <= width && width <= 80);
   1.140 +      memcpy(str, dsa->card + pos, width), str[width] = '\0';
   1.141 +      if (str2int(strspx(str), val))
   1.142 +      {  xprintf("%s:%d: field `%s' contains invalid value `%s'\n",
   1.143 +            dsa->fname, dsa->seqn, fld, str);
   1.144 +         return 1;
   1.145 +      }
   1.146 +      return 0;
   1.147 +}
   1.148 +
   1.149 +/***********************************************************************
   1.150 +*  parse_fmt - parse Fortran format specification
   1.151 +*
   1.152 +*  This routine parses the Fortran format specification represented as
   1.153 +*  character string which fmt points to and stores format elements into
   1.154 +*  appropriate static locations. Should note that not all valid Fortran
   1.155 +*  format specifications may be recognised. If the format specification
   1.156 +*  was recognised, the routine returns zero, otherwise non-zero. */
   1.157 +
   1.158 +static int parse_fmt(struct dsa *dsa, char *fmt)
   1.159 +{     int k, s, val;
   1.160 +      char str[80+1];
   1.161 +      /* first character should be left parenthesis */
   1.162 +      if (fmt[0] != '(')
   1.163 +fail: {  xprintf("hbm_read_mat: format `%s' not recognised\n", fmt);
   1.164 +         return 1;
   1.165 +      }
   1.166 +      k = 1;
   1.167 +      /* optional scale factor */
   1.168 +      dsa->fmt_p = 0;
   1.169 +      if (isdigit((unsigned char)fmt[k]))
   1.170 +      {  s = 0;
   1.171 +         while (isdigit((unsigned char)fmt[k]))
   1.172 +         {  if (s == 80) goto fail;
   1.173 +            str[s++] = fmt[k++];
   1.174 +         }
   1.175 +         str[s] = '\0';
   1.176 +         if (str2int(str, &val)) goto fail;
   1.177 +         if (toupper((unsigned char)fmt[k]) != 'P') goto iter;
   1.178 +         dsa->fmt_p = val, k++;
   1.179 +         if (!(0 <= dsa->fmt_p && dsa->fmt_p <= 255)) goto fail;
   1.180 +         /* optional comma may follow scale factor */
   1.181 +         if (fmt[k] == ',') k++;
   1.182 +      }
   1.183 +      /* optional iterator */
   1.184 +      dsa->fmt_k = 1;
   1.185 +      if (isdigit((unsigned char)fmt[k]))
   1.186 +      {  s = 0;
   1.187 +         while (isdigit((unsigned char)fmt[k]))
   1.188 +         {  if (s == 80) goto fail;
   1.189 +            str[s++] = fmt[k++];
   1.190 +         }
   1.191 +         str[s] = '\0';
   1.192 +         if (str2int(str, &val)) goto fail;
   1.193 +iter:    dsa->fmt_k = val;
   1.194 +         if (!(1 <= dsa->fmt_k && dsa->fmt_k <= 255)) goto fail;
   1.195 +      }
   1.196 +      /* format code */
   1.197 +      dsa->fmt_f = toupper((unsigned char)fmt[k++]);
   1.198 +      if (!(dsa->fmt_f == 'D' || dsa->fmt_f == 'E' ||
   1.199 +            dsa->fmt_f == 'F' || dsa->fmt_f == 'G' ||
   1.200 +            dsa->fmt_f == 'I')) goto fail;
   1.201 +      /* field width */
   1.202 +      if (!isdigit((unsigned char)fmt[k])) goto fail;
   1.203 +      s = 0;
   1.204 +      while (isdigit((unsigned char)fmt[k]))
   1.205 +      {  if (s == 80) goto fail;
   1.206 +         str[s++] = fmt[k++];
   1.207 +      }
   1.208 +      str[s] = '\0';
   1.209 +      if (str2int(str, &dsa->fmt_w)) goto fail;
   1.210 +      if (!(1 <= dsa->fmt_w && dsa->fmt_w <= 255)) goto fail;
   1.211 +      /* optional number of decimal places after point */
   1.212 +      dsa->fmt_d = 0;
   1.213 +      if (fmt[k] == '.')
   1.214 +      {  k++;
   1.215 +         if (!isdigit((unsigned char)fmt[k])) goto fail;
   1.216 +         s = 0;
   1.217 +         while (isdigit((unsigned char)fmt[k]))
   1.218 +         {  if (s == 80) goto fail;
   1.219 +            str[s++] = fmt[k++];
   1.220 +         }
   1.221 +         str[s] = '\0';
   1.222 +         if (str2int(str, &dsa->fmt_d)) goto fail;
   1.223 +         if (!(0 <= dsa->fmt_d && dsa->fmt_d <= 255)) goto fail;
   1.224 +      }
   1.225 +      /* last character should be right parenthesis */
   1.226 +      if (!(fmt[k] == ')' && fmt[k+1] == '\0')) goto fail;
   1.227 +      return 0;
   1.228 +}
   1.229 +
   1.230 +/***********************************************************************
   1.231 +*  read_int_array - read array of integer type
   1.232 +*
   1.233 +*  This routine reads an integer array from the input text file, where
   1.234 +*  name is array name, fmt is Fortran format specification that controls
   1.235 +*  reading, n is number of array elements, val is array of integer type.
   1.236 +*  If the array was read successful, the routine returns zero, otherwise
   1.237 +*  non-zero. */
   1.238 +
   1.239 +static int read_int_array(struct dsa *dsa, char *name, char *fmt,
   1.240 +      int n, int val[])
   1.241 +{     int k, pos;
   1.242 +      char str[80+1];
   1.243 +      if (parse_fmt(dsa, fmt)) return 1;
   1.244 +      if (!(dsa->fmt_f == 'I' && dsa->fmt_w <= 80 &&
   1.245 +            dsa->fmt_k * dsa->fmt_w <= 80))
   1.246 +      {  xprintf(
   1.247 +            "%s:%d: can't read array `%s' - invalid format `%s'\n",
   1.248 +            dsa->fname, dsa->seqn, name, fmt);
   1.249 +         return 1;
   1.250 +      }
   1.251 +      for (k = 1, pos = INT_MAX; k <= n; k++, pos++)
   1.252 +      {  if (pos >= dsa->fmt_k)
   1.253 +         {  if (read_card(dsa)) return 1;
   1.254 +            pos = 0;
   1.255 +         }
   1.256 +         memcpy(str, dsa->card + dsa->fmt_w * pos, dsa->fmt_w);
   1.257 +         str[dsa->fmt_w] = '\0';
   1.258 +         strspx(str);
   1.259 +         if (str2int(str, &val[k]))
   1.260 +         {  xprintf(
   1.261 +               "%s:%d: can't read array `%s' - invalid value `%s'\n",
   1.262 +               dsa->fname, dsa->seqn, name, str);
   1.263 +            return 1;
   1.264 +         }
   1.265 +      }
   1.266 +      return 0;
   1.267 +}
   1.268 +
   1.269 +/***********************************************************************
   1.270 +*  read_real_array - read array of real type
   1.271 +*
   1.272 +*  This routine reads a real array from the input text file, where name
   1.273 +*  is array name, fmt is Fortran format specification that controls
   1.274 +*  reading, n is number of array elements, val is array of real type.
   1.275 +*  If the array was read successful, the routine returns zero, otherwise
   1.276 +*  non-zero. */
   1.277 +
   1.278 +static int read_real_array(struct dsa *dsa, char *name, char *fmt,
   1.279 +      int n, double val[])
   1.280 +{     int k, pos;
   1.281 +      char str[80+1], *ptr;
   1.282 +      if (parse_fmt(dsa, fmt)) return 1;
   1.283 +      if (!(dsa->fmt_f != 'I' && dsa->fmt_w <= 80 &&
   1.284 +            dsa->fmt_k * dsa->fmt_w <= 80))
   1.285 +      {  xprintf(
   1.286 +            "%s:%d: can't read array `%s' - invalid format `%s'\n",
   1.287 +            dsa->fname, dsa->seqn, name, fmt);
   1.288 +         return 1;
   1.289 +      }
   1.290 +      for (k = 1, pos = INT_MAX; k <= n; k++, pos++)
   1.291 +      {  if (pos >= dsa->fmt_k)
   1.292 +         {  if (read_card(dsa)) return 1;
   1.293 +            pos = 0;
   1.294 +         }
   1.295 +         memcpy(str, dsa->card + dsa->fmt_w * pos, dsa->fmt_w);
   1.296 +         str[dsa->fmt_w] = '\0';
   1.297 +         strspx(str);
   1.298 +         if (strchr(str, '.') == NULL && strcmp(str, "0"))
   1.299 +         {  xprintf("%s(%d): can't read array `%s' - value `%s' has no "
   1.300 +               "decimal point\n", dsa->fname, dsa->seqn, name, str);
   1.301 +            return 1;
   1.302 +         }
   1.303 +         /* sometimes lower case letters appear */
   1.304 +         for (ptr = str; *ptr; ptr++)
   1.305 +            *ptr = (char)toupper((unsigned char)*ptr);
   1.306 +         ptr = strchr(str, 'D');
   1.307 +         if (ptr != NULL) *ptr = 'E';
   1.308 +         /* value may appear with decimal exponent but without letters
   1.309 +            E or D (for example, -123.456-012), so missing letter should
   1.310 +            be inserted */
   1.311 +         ptr = strchr(str+1, '+');
   1.312 +         if (ptr == NULL) ptr = strchr(str+1, '-');
   1.313 +         if (ptr != NULL && *(ptr-1) != 'E')
   1.314 +         {  xassert(strlen(str) < 80);
   1.315 +            memmove(ptr+1, ptr, strlen(ptr)+1);
   1.316 +            *ptr = 'E';
   1.317 +         }
   1.318 +         if (str2num(str, &val[k]))
   1.319 +         {  xprintf(
   1.320 +               "%s:%d: can't read array `%s' - invalid value `%s'\n",
   1.321 +               dsa->fname, dsa->seqn, name, str);
   1.322 +            return 1;
   1.323 +         }
   1.324 +      }
   1.325 +      return 0;
   1.326 +}
   1.327 +
   1.328 +HBM *hbm_read_mat(const char *fname)
   1.329 +{     struct dsa _dsa, *dsa = &_dsa;
   1.330 +      HBM *hbm = NULL;
   1.331 +      dsa->fname = fname;
   1.332 +      xprintf("hbm_read_mat: reading matrix from `%s'...\n",
   1.333 +         dsa->fname);
   1.334 +      dsa->fp = fopen(dsa->fname, "r");
   1.335 +      if (dsa->fp == NULL)
   1.336 +      {  xprintf("hbm_read_mat: unable to open `%s' - %s\n",
   1.337 +            dsa->fname, strerror(errno));
   1.338 +         goto fail;
   1.339 +      }
   1.340 +      dsa->seqn = 0;
   1.341 +      hbm = xmalloc(sizeof(HBM));
   1.342 +      memset(hbm, 0, sizeof(HBM));
   1.343 +      /* read the first heading card */
   1.344 +      if (read_card(dsa)) goto fail;
   1.345 +      memcpy(hbm->title, dsa->card, 72), hbm->title[72] = '\0';
   1.346 +      strtrim(hbm->title);
   1.347 +      xprintf("%s\n", hbm->title);
   1.348 +      memcpy(hbm->key, dsa->card+72, 8), hbm->key[8] = '\0';
   1.349 +      strspx(hbm->key);
   1.350 +      xprintf("key = %s\n", hbm->key);
   1.351 +      /* read the second heading card */
   1.352 +      if (read_card(dsa)) goto fail;
   1.353 +      if (scan_int(dsa, "totcrd",  0, 14, &hbm->totcrd)) goto fail;
   1.354 +      if (scan_int(dsa, "ptrcrd", 14, 14, &hbm->ptrcrd)) goto fail;
   1.355 +      if (scan_int(dsa, "indcrd", 28, 14, &hbm->indcrd)) goto fail;
   1.356 +      if (scan_int(dsa, "valcrd", 42, 14, &hbm->valcrd)) goto fail;
   1.357 +      if (scan_int(dsa, "rhscrd", 56, 14, &hbm->rhscrd)) goto fail;
   1.358 +      xprintf("totcrd = %d; ptrcrd = %d; indcrd = %d; valcrd = %d; rhsc"
   1.359 +         "rd = %d\n", hbm->totcrd, hbm->ptrcrd, hbm->indcrd,
   1.360 +         hbm->valcrd, hbm->rhscrd);
   1.361 +      /* read the third heading card */
   1.362 +      if (read_card(dsa)) goto fail;
   1.363 +      memcpy(hbm->mxtype, dsa->card, 3), hbm->mxtype[3] = '\0';
   1.364 +      if (strchr("RCP",   hbm->mxtype[0]) == NULL ||
   1.365 +          strchr("SUHZR", hbm->mxtype[1]) == NULL ||
   1.366 +          strchr("AE",    hbm->mxtype[2]) == NULL)
   1.367 +      {  xprintf("%s:%d: matrix type `%s' not recognised\n",
   1.368 +            dsa->fname, dsa->seqn, hbm->mxtype);
   1.369 +         goto fail;
   1.370 +      }
   1.371 +      if (scan_int(dsa, "nrow", 14, 14, &hbm->nrow)) goto fail;
   1.372 +      if (scan_int(dsa, "ncol", 28, 14, &hbm->ncol)) goto fail;
   1.373 +      if (scan_int(dsa, "nnzero", 42, 14, &hbm->nnzero)) goto fail;
   1.374 +      if (scan_int(dsa, "neltvl", 56, 14, &hbm->neltvl)) goto fail;
   1.375 +      xprintf("mxtype = %s; nrow = %d; ncol = %d; nnzero = %d; neltvl ="
   1.376 +         " %d\n", hbm->mxtype, hbm->nrow, hbm->ncol, hbm->nnzero,
   1.377 +         hbm->neltvl);
   1.378 +      /* read the fourth heading card */
   1.379 +      if (read_card(dsa)) goto fail;
   1.380 +      memcpy(hbm->ptrfmt, dsa->card, 16), hbm->ptrfmt[16] = '\0';
   1.381 +      strspx(hbm->ptrfmt);
   1.382 +      memcpy(hbm->indfmt, dsa->card+16, 16), hbm->indfmt[16] = '\0';
   1.383 +      strspx(hbm->indfmt);
   1.384 +      memcpy(hbm->valfmt, dsa->card+32, 20), hbm->valfmt[20] = '\0';
   1.385 +      strspx(hbm->valfmt);
   1.386 +      memcpy(hbm->rhsfmt, dsa->card+52, 20), hbm->rhsfmt[20] = '\0';
   1.387 +      strspx(hbm->rhsfmt);
   1.388 +      xprintf("ptrfmt = %s; indfmt = %s; valfmt = %s; rhsfmt = %s\n",
   1.389 +         hbm->ptrfmt, hbm->indfmt, hbm->valfmt, hbm->rhsfmt);
   1.390 +      /* read the fifth heading card (optional) */
   1.391 +      if (hbm->rhscrd <= 0)
   1.392 +      {  strcpy(hbm->rhstyp, "???");
   1.393 +         hbm->nrhs = 0;
   1.394 +         hbm->nrhsix = 0;
   1.395 +      }
   1.396 +      else
   1.397 +      {  if (read_card(dsa)) goto fail;
   1.398 +         memcpy(hbm->rhstyp, dsa->card, 3), hbm->rhstyp[3] = '\0';
   1.399 +         if (scan_int(dsa, "nrhs", 14, 14, &hbm->nrhs)) goto fail;
   1.400 +         if (scan_int(dsa, "nrhsix", 28, 14, &hbm->nrhsix)) goto fail;
   1.401 +         xprintf("rhstyp = `%s'; nrhs = %d; nrhsix = %d\n",
   1.402 +            hbm->rhstyp, hbm->nrhs, hbm->nrhsix);
   1.403 +      }
   1.404 +      /* read matrix structure */
   1.405 +      hbm->colptr = xcalloc(1+hbm->ncol+1, sizeof(int));
   1.406 +      if (read_int_array(dsa, "colptr", hbm->ptrfmt, hbm->ncol+1,
   1.407 +         hbm->colptr)) goto fail;
   1.408 +      hbm->rowind = xcalloc(1+hbm->nnzero, sizeof(int));
   1.409 +      if (read_int_array(dsa, "rowind", hbm->indfmt, hbm->nnzero,
   1.410 +         hbm->rowind)) goto fail;
   1.411 +      /* read matrix values */
   1.412 +      if (hbm->valcrd <= 0) goto done;
   1.413 +      if (hbm->mxtype[2] == 'A')
   1.414 +      {  /* assembled matrix */
   1.415 +         hbm->values = xcalloc(1+hbm->nnzero, sizeof(double));
   1.416 +         if (read_real_array(dsa, "values", hbm->valfmt, hbm->nnzero,
   1.417 +            hbm->values)) goto fail;
   1.418 +      }
   1.419 +      else
   1.420 +      {  /* elemental (unassembled) matrix */
   1.421 +         hbm->values = xcalloc(1+hbm->neltvl, sizeof(double));
   1.422 +         if (read_real_array(dsa, "values", hbm->valfmt, hbm->neltvl,
   1.423 +            hbm->values)) goto fail;
   1.424 +      }
   1.425 +      /* read right-hand sides */
   1.426 +      if (hbm->nrhs <= 0) goto done;
   1.427 +      if (hbm->rhstyp[0] == 'F')
   1.428 +      {  /* dense format */
   1.429 +         hbm->nrhsvl = hbm->nrow * hbm->nrhs;
   1.430 +         hbm->rhsval = xcalloc(1+hbm->nrhsvl, sizeof(double));
   1.431 +         if (read_real_array(dsa, "rhsval", hbm->rhsfmt, hbm->nrhsvl,
   1.432 +            hbm->rhsval)) goto fail;
   1.433 +      }
   1.434 +      else if (hbm->rhstyp[0] == 'M' && hbm->mxtype[2] == 'A')
   1.435 +      {  /* sparse format */
   1.436 +         /* read pointers */
   1.437 +         hbm->rhsptr = xcalloc(1+hbm->nrhs+1, sizeof(int));
   1.438 +         if (read_int_array(dsa, "rhsptr", hbm->ptrfmt, hbm->nrhs+1,
   1.439 +            hbm->rhsptr)) goto fail;
   1.440 +         /* read sparsity pattern */
   1.441 +         hbm->rhsind = xcalloc(1+hbm->nrhsix, sizeof(int));
   1.442 +         if (read_int_array(dsa, "rhsind", hbm->indfmt, hbm->nrhsix,
   1.443 +            hbm->rhsind)) goto fail;
   1.444 +         /* read values */
   1.445 +         hbm->rhsval = xcalloc(1+hbm->nrhsix, sizeof(double));
   1.446 +         if (read_real_array(dsa, "rhsval", hbm->rhsfmt, hbm->nrhsix,
   1.447 +            hbm->rhsval)) goto fail;
   1.448 +      }
   1.449 +      else if (hbm->rhstyp[0] == 'M' && hbm->mxtype[2] == 'E')
   1.450 +      {  /* elemental format */
   1.451 +         hbm->rhsval = xcalloc(1+hbm->nrhsvl, sizeof(double));
   1.452 +         if (read_real_array(dsa, "rhsval", hbm->rhsfmt, hbm->nrhsvl,
   1.453 +            hbm->rhsval)) goto fail;
   1.454 +      }
   1.455 +      else
   1.456 +      {  xprintf("%s:%d: right-hand side type `%c' not recognised\n",
   1.457 +            dsa->fname, dsa->seqn, hbm->rhstyp[0]);
   1.458 +         goto fail;
   1.459 +      }
   1.460 +      /* read starting guesses */
   1.461 +      if (hbm->rhstyp[1] == 'G')
   1.462 +      {  hbm->nguess = hbm->nrow * hbm->nrhs;
   1.463 +         hbm->sguess = xcalloc(1+hbm->nguess, sizeof(double));
   1.464 +         if (read_real_array(dsa, "sguess", hbm->rhsfmt, hbm->nguess,
   1.465 +            hbm->sguess)) goto fail;
   1.466 +      }
   1.467 +      /* read solution vectors */
   1.468 +      if (hbm->rhstyp[2] == 'X')
   1.469 +      {  hbm->nexact = hbm->nrow * hbm->nrhs;
   1.470 +         hbm->xexact = xcalloc(1+hbm->nexact, sizeof(double));
   1.471 +         if (read_real_array(dsa, "xexact", hbm->rhsfmt, hbm->nexact,
   1.472 +            hbm->xexact)) goto fail;
   1.473 +      }
   1.474 +done: /* reading has been completed */
   1.475 +      xprintf("hbm_read_mat: %d cards were read\n", dsa->seqn);
   1.476 +      fclose(dsa->fp);
   1.477 +      return hbm;
   1.478 +fail: /* something wrong in Danish kingdom */
   1.479 +      if (hbm != NULL)
   1.480 +      {  if (hbm->colptr != NULL) xfree(hbm->colptr);
   1.481 +         if (hbm->rowind != NULL) xfree(hbm->rowind);
   1.482 +         if (hbm->rhsptr != NULL) xfree(hbm->rhsptr);
   1.483 +         if (hbm->rhsind != NULL) xfree(hbm->rhsind);
   1.484 +         if (hbm->values != NULL) xfree(hbm->values);
   1.485 +         if (hbm->rhsval != NULL) xfree(hbm->rhsval);
   1.486 +         if (hbm->sguess != NULL) xfree(hbm->sguess);
   1.487 +         if (hbm->xexact != NULL) xfree(hbm->xexact);
   1.488 +         xfree(hbm);
   1.489 +      }
   1.490 +      if (dsa->fp != NULL) fclose(dsa->fp);
   1.491 +      return NULL;
   1.492 +}
   1.493 +
   1.494 +/***********************************************************************
   1.495 +*  NAME
   1.496 +*
   1.497 +*  hbm_free_mat - free sparse matrix in Harwell-Boeing format
   1.498 +*
   1.499 +*  SYNOPSIS
   1.500 +*
   1.501 +*  #include "glphbm.h"
   1.502 +*  void hbm_free_mat(HBM *hbm);
   1.503 +*
   1.504 +*  DESCRIPTION
   1.505 +*
   1.506 +*  The hbm_free_mat routine frees all the memory allocated to the data
   1.507 +*  structure containing a sparse matrix in the Harwell-Boeing format. */
   1.508 +
   1.509 +void hbm_free_mat(HBM *hbm)
   1.510 +{     if (hbm->colptr != NULL) xfree(hbm->colptr);
   1.511 +      if (hbm->rowind != NULL) xfree(hbm->rowind);
   1.512 +      if (hbm->rhsptr != NULL) xfree(hbm->rhsptr);
   1.513 +      if (hbm->rhsind != NULL) xfree(hbm->rhsind);
   1.514 +      if (hbm->values != NULL) xfree(hbm->values);
   1.515 +      if (hbm->rhsval != NULL) xfree(hbm->rhsval);
   1.516 +      if (hbm->sguess != NULL) xfree(hbm->sguess);
   1.517 +      if (hbm->xexact != NULL) xfree(hbm->xexact);
   1.518 +      xfree(hbm);
   1.519 +      return;
   1.520 +}
   1.521 +
   1.522 +/* eof */