/* glphbm.c */

/*----------------------------------------------------------------------
-- This code is part of GNU Linear Programming Kit (GLPK).
--
-- Copyright (C) 2000, 01, 02, 03, 04, 05, 06 Andrew Makhorin,
-- Department for Applied Informatics, Moscow Aviation Institute,
-- Moscow, Russia. All rights reserved. E-mail: <mao@mai2.rcnet.ru>.
--
-- GLPK is free software; you can redistribute it and/or modify it
-- under the terms of the GNU General Public License as published by
-- the Free Software Foundation; either version 2, or (at your option)
-- any later version.
--
-- GLPK is distributed in the hope that it will be useful, but WITHOUT
-- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
-- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
-- License for more details.
--
-- You should have received a copy of the GNU General Public License
-- along with GLPK; see the file COPYING. If not, write to the Free
-- Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
-- 02110-1301, USA.
----------------------------------------------------------------------*/

#include <ctype.h>
#include <errno.h>
#include <limits.h>
#include <stdio.h>
#include <string.h>
#include "glphbm.h"
#include "glplib.h"

/*----------------------------------------------------------------------
-- hbm_read_mat - read sparse matrix in Harwell-Boeing format.
--
-- *Synopsis*
--
-- #include "glphbm.h"
-- HBM *hbm_read_mat(char *fname);
--
-- *Description*
--
-- The routine hbm_read_mat reads a sparse matrix in the Harwell-Boeing
-- format from a text file whose name is the character string fname.
--
-- Detailed description of the Harwell-Boeing format recognised by this
-- routine is given in the following report:
--
-- I.S.Duff, R.G.Grimes, J.G.Lewis. User's Guide for the Harwell-Boeing
-- Sparse Matrix Collection (Release I), TR/PA/92/86, October 1992.
--
-- *Returns*
--
-- If no error occured, the routine hbm_read_mat returns a pointer to
-- a data structure containing the matrix. In case of error the routine
-- prints an appropriate error message and returns NULL. */

struct dsa
{     /* working area used by routine hbm_read_mat */
      char *fname;
      /* name of input text file */
      FILE *fp;
      /* stream assigned to input text file */
      int seqn;
      /* card sequential number */
      char card[80+1];
      /* card image buffer */
      int fmt_p;
      /* scale factor */
      int fmt_k;
      /* iterator */
      int fmt_f;
      /* format code */
      int fmt_w;
      /* field width */
      int fmt_d;
      /* number of decimal places after point */
};

/*----------------------------------------------------------------------
-- read_card - read next data card.
--
-- This routine reads the next 80-column card from the input text file
-- and stores its image into the character string card. If the card was
-- read successfully, the routine returns zero, otherwise non-zero. */

static int read_card(struct dsa *dsa)
{     int k, c;
      dsa->seqn++;
      memset(dsa->card, ' ', 80), dsa->card[80] = '\0';
      k = 0;
      for (;;)
      {  c = fgetc(dsa->fp);
         if (ferror(dsa->fp))
         {  print("%s:%d: read error - %s", dsa->fname, dsa->seqn,
               strerror(errno));
            return 1;
         }
         if (feof(dsa->fp))
         {  if (k == 0)
               print("%s:%d: unexpected EOF", dsa->fname, dsa->seqn);
            else
               print("%s:%d: missing final LF", dsa->fname, dsa->seqn);
            return 1;
         }
         if (c == '\r') continue;
         if (c == '\n') break;
         if (iscntrl(c))
         {  print("%s:%d: invalid control character 0x%02X", dsa->fname,
               dsa->seqn, c);
            return 1;
         }
         if (k == 80)
         {  print("%s:%d: card image too long", dsa->fname, dsa->seqn);
            return 1;
         }
         dsa->card[k++] = (char)c;
      }
      return 0;
}

/*----------------------------------------------------------------------
-- scan_int - scan integer value from the current card.
--
-- This routine scans an integer value from the current card, where fld
-- is the name of the field, pos is the position of the field, width is
-- the width of the field, val points to a location to which the scanned
-- value should be stored. If the value was scanned successfully, the
-- routine returns zero, otherwise non-zero. */

static int scan_int(struct dsa *dsa, char *fld, int pos, int width,
      int *val)
{     char str[80+1];
      insist(1 <= width && width <= 80);
      memcpy(str, dsa->card + pos, width), str[width] = '\0';
      if (str2int(strspx(str), val))
      {  print("%s:%d: field `%s' contains invalid value `%s'",
            dsa->fname, dsa->seqn, fld, str);
         return 1;
      }
      return 0;
}

/*----------------------------------------------------------------------
-- parse_fmt - parse Fortran format specification.
--
-- This routine parses the Fortran format specification represented as
-- character string which fmt points to and stores format elements into
-- appropriate static locations. Should note that not all valid Fortran
-- format specifications may be recognised. If the format specification
-- was recognised, the routine returns zero, otherwise non-zero. */

static int parse_fmt(struct dsa *dsa, char *fmt)
{     int k, s, val;
      char str[80+1];
      /* first character should be left parenthesis */
      if (fmt[0] != '(')
fail: {  print("hbm_read_mat: format `%s' not recognised", fmt);
         return 1;
      }
      k = 1;
      /* optional scale factor */
      dsa->fmt_p = 0;
      if (isdigit((unsigned char)fmt[k]))
      {  s = 0;
         while (isdigit((unsigned char)fmt[k]))
         {  if (s == 80) goto fail;
            str[s++] = fmt[k++];
         }
         str[s] = '\0';
         if (str2int(str, &val)) goto fail;
         if (toupper((unsigned char)fmt[k]) != 'P') goto iter;
         dsa->fmt_p = val, k++;
         if (!(0 <= dsa->fmt_p && dsa->fmt_p <= 255)) goto fail;
         /* optional comma may follow scale factor */
         if (fmt[k] == ',') k++;
      }
      /* optional iterator */
      dsa->fmt_k = 1;
      if (isdigit((unsigned char)fmt[k]))
      {  s = 0;
         while (isdigit((unsigned char)fmt[k]))
         {  if (s == 80) goto fail;
            str[s++] = fmt[k++];
         }
         str[s] = '\0';
         if (str2int(str, &val)) goto fail;
iter:    dsa->fmt_k = val;
         if (!(1 <= dsa->fmt_k && dsa->fmt_k <= 255)) goto fail;
      }
      /* format code */
      dsa->fmt_f = toupper((unsigned char)fmt[k++]);
      if (!(dsa->fmt_f == 'D' || dsa->fmt_f == 'E' ||
            dsa->fmt_f == 'F' || dsa->fmt_f == 'G' ||
            dsa->fmt_f == 'I')) goto fail;
      /* field width */
      if (!isdigit((unsigned char)fmt[k])) goto fail;
      s = 0;
      while (isdigit((unsigned char)fmt[k]))
      {  if (s == 80) goto fail;
         str[s++] = fmt[k++];
      }
      str[s] = '\0';
      if (str2int(str, &dsa->fmt_w)) goto fail;
      if (!(1 <= dsa->fmt_w && dsa->fmt_w <= 255)) goto fail;
      /* optional number of decimal places after point */
      dsa->fmt_d = 0;
      if (fmt[k] == '.')
      {  k++;
         if (!isdigit((unsigned char)fmt[k])) goto fail;
         s = 0;
         while (isdigit((unsigned char)fmt[k]))
         {  if (s == 80) goto fail;
            str[s++] = fmt[k++];
         }
         str[s] = '\0';
         if (str2int(str, &dsa->fmt_d)) goto fail;
         if (!(0 <= dsa->fmt_d && dsa->fmt_d <= 255)) goto fail;
      }
      /* last character should be right parenthesis */
      if (!(fmt[k] == ')' && fmt[k+1] == '\0')) goto fail;
      return 0;
}

/*----------------------------------------------------------------------
-- read_int_array - read array of integer type.
--
-- This routine reads an integer array from the input text file, where
-- name is array name, fmt is Fortran format specification that controls
-- reading, n is number of array elements, val is array of integer type.
-- If the array was read successful, the routine returns zero, otherwise
-- non-zero. */

static int read_int_array(struct dsa *dsa, char *name, char *fmt,
      int n, int val[])
{     int k, pos;
      char str[80+1];
      if (parse_fmt(dsa, fmt)) return 1;
      if (!(dsa->fmt_f == 'I' && dsa->fmt_w <= 80 &&
            dsa->fmt_k * dsa->fmt_w <= 80))
      {  print("%s:%d: can't read array `%s' - invalid format `%s'",
            dsa->fname, dsa->seqn, name, fmt);
         return 1;
      }
      for (k = 1, pos = INT_MAX; k <= n; k++, pos++)
      {  if (pos >= dsa->fmt_k)
         {  if (read_card(dsa)) return 1;
            pos = 0;
         }
         memcpy(str, dsa->card + dsa->fmt_w * pos, dsa->fmt_w);
         str[dsa->fmt_w] = '\0';
         strspx(str);
         if (str2int(str, &val[k]))
         {  print("%s:%d: can't read array `%s' - invalid value `%s'",
               dsa->fname, dsa->seqn, name, str);
            return 1;
         }
      }
      return 0;
}

/*----------------------------------------------------------------------
-- read_real_array - read array of real type.
--
-- This routine reads a real array from the input text file, where name
-- is array name, fmt is Fortran format specification that controls
-- reading, n is number of array elements, val is array of real type.
-- If the array was read successful, the routine returns zero, otherwise
-- non-zero. */

static int read_real_array(struct dsa *dsa, char *name, char *fmt,
      int n, double val[])
{     int k, pos;
      char str[80+1], *ptr;
      if (parse_fmt(dsa, fmt)) return 1;
      if (!(dsa->fmt_f != 'I' && dsa->fmt_w <= 80 &&
            dsa->fmt_k * dsa->fmt_w <= 80))
      {  print("%s:%d: can't read array `%s' - invalid format `%s'",
            dsa->fname, dsa->seqn, name, fmt);
         return 1;
      }
      for (k = 1, pos = INT_MAX; k <= n; k++, pos++)
      {  if (pos >= dsa->fmt_k)
         {  if (read_card(dsa)) return 1;
            pos = 0;
         }
         memcpy(str, dsa->card + dsa->fmt_w * pos, dsa->fmt_w);
         str[dsa->fmt_w] = '\0';
         strspx(str);
         if (strchr(str, '.') == NULL && strcmp(str, "0"))
         {  print("%s(%d): can't read array `%s' - value `%s' has no de"
               "cimal point", dsa->fname, dsa->seqn, name, str);
            return 1;
         }
         /* sometimes lower case letters appear */
         for (ptr = str; *ptr; ptr++)
            *ptr = (char)toupper((unsigned char)*ptr);
         ptr = strchr(str, 'D');
         if (ptr != NULL) *ptr = 'E';
         /* value may appear with decimal exponent but without letters
            E or D (for example, -123.456-012), so missing letter should
            be inserted */
         ptr = strchr(str+1, '+');
         if (ptr == NULL) ptr = strchr(str+1, '-');
         if (ptr != NULL && *(ptr-1) != 'E')
         {  insist(strlen(str) < 80);
            memmove(ptr+1, ptr, strlen(ptr)+1);
            *ptr = 'E';
         }
         if (str2dbl(str, &val[k]))
         {  print("%s:%d: can't read array `%s' - invalid value `%s'",
               dsa->fname, dsa->seqn, name, str);
            return 1;
         }
      }
      return 0;
}

HBM *hbm_read_mat(char *fname)
{     struct dsa _dsa, *dsa = &_dsa;
      HBM *hbm = NULL;
      dsa->fname = fname;
      print("hbm_read_mat: reading matrix from `%s'...", dsa->fname);
      dsa->fp = ufopen(dsa->fname, "r");
      if (dsa->fp == NULL)
      {  print("hbm_read_mat: unable to open `%s' - %s", dsa->fname,
            strerror(errno));
         goto fail;
      }
      dsa->seqn = 0;
      hbm = umalloc(sizeof(HBM));
      memset(hbm, 0, sizeof(HBM));
      /* read the first heading card */
      if (read_card(dsa)) goto fail;
      memcpy(hbm->title, dsa->card, 72), hbm->title[72] = '\0';
      strtrim(hbm->title);
      print("%s", hbm->title);
      memcpy(hbm->key, dsa->card+72, 8), hbm->key[8] = '\0';
      strspx(hbm->key);
      print("key = %s", hbm->key);
      /* read the second heading card */
      if (read_card(dsa)) goto fail;
      if (scan_int(dsa, "totcrd",  0, 14, &hbm->totcrd)) goto fail;
      if (scan_int(dsa, "ptrcrd", 14, 14, &hbm->ptrcrd)) goto fail;
      if (scan_int(dsa, "indcrd", 28, 14, &hbm->indcrd)) goto fail;
      if (scan_int(dsa, "valcrd", 42, 14, &hbm->valcrd)) goto fail;
      if (scan_int(dsa, "rhscrd", 56, 14, &hbm->rhscrd)) goto fail;
      print("totcrd = %d; ptrcrd = %d; indcrd = %d; valcrd = %d; rhscrd"
         " = %d", hbm->totcrd, hbm->ptrcrd, hbm->indcrd, hbm->valcrd,
         hbm->rhscrd);
      /* read the third heading card */
      if (read_card(dsa)) goto fail;
      memcpy(hbm->mxtype, dsa->card, 3), hbm->mxtype[3] = '\0';
      if (strchr("RCP",   hbm->mxtype[0]) == NULL ||
          strchr("SUHZR", hbm->mxtype[1]) == NULL ||
          strchr("AE",    hbm->mxtype[2]) == NULL)
      {  print("%s:%d: matrix type `%s' not recognised", dsa->fname,
            dsa->seqn, hbm->mxtype);
         goto fail;
      }
      if (scan_int(dsa, "nrow", 14, 14, &hbm->nrow)) goto fail;
      if (scan_int(dsa, "ncol", 28, 14, &hbm->ncol)) goto fail;
      if (scan_int(dsa, "nnzero", 42, 14, &hbm->nnzero)) goto fail;
      if (scan_int(dsa, "neltvl", 56, 14, &hbm->neltvl)) goto fail;
      print("mxtype = %s; nrow = %d; ncol = %d; nnzero = %d; neltvl = %"
         "d", hbm->mxtype, hbm->nrow, hbm->ncol, hbm->nnzero,
         hbm->neltvl);
      /* read the fourth heading card */
      if (read_card(dsa)) goto fail;
      memcpy(hbm->ptrfmt, dsa->card, 16), hbm->ptrfmt[16] = '\0';
      strspx(hbm->ptrfmt);
      memcpy(hbm->indfmt, dsa->card+16, 16), hbm->indfmt[16] = '\0';
      strspx(hbm->indfmt);
      memcpy(hbm->valfmt, dsa->card+32, 20), hbm->valfmt[20] = '\0';
      strspx(hbm->valfmt);
      memcpy(hbm->rhsfmt, dsa->card+52, 20), hbm->rhsfmt[20] = '\0';
      strspx(hbm->rhsfmt);
      print("ptrfmt = %s; indfmt = %s; valfmt = %s; rhsfmt = %s",
         hbm->ptrfmt, hbm->indfmt, hbm->valfmt, hbm->rhsfmt);
      /* read the fifth heading card (optional) */
      if (hbm->rhscrd <= 0)
      {  strcpy(hbm->rhstyp, "???");
         hbm->nrhs = 0;
         hbm->nrhsix = 0;
      }
      else
      {  if (read_card(dsa)) goto fail;
         memcpy(hbm->rhstyp, dsa->card, 3), hbm->rhstyp[3] = '\0';
         if (scan_int(dsa, "nrhs", 14, 14, &hbm->nrhs)) goto fail;
         if (scan_int(dsa, "nrhsix", 28, 14, &hbm->nrhsix)) goto fail;
         print("rhstyp = `%s'; nrhs = %d; nrhsix = %d", hbm->rhstyp,
            hbm->nrhs, hbm->nrhsix);
      }
      /* read matrix structure */
      hbm->colptr = ucalloc(1+hbm->ncol+1, sizeof(int));
      if (read_int_array(dsa, "colptr", hbm->ptrfmt, hbm->ncol+1,
         hbm->colptr)) goto fail;
      hbm->rowind = ucalloc(1+hbm->nnzero, sizeof(int));
      if (read_int_array(dsa, "rowind", hbm->indfmt, hbm->nnzero,
         hbm->rowind)) goto fail;
      /* read matrix values */
      if (hbm->valcrd <= 0) goto done;
      if (hbm->mxtype[2] == 'A')
      {  /* assembled matrix */
         hbm->values = ucalloc(1+hbm->nnzero, sizeof(double));
         if (read_real_array(dsa, "values", hbm->valfmt, hbm->nnzero,
            hbm->values)) goto fail;
      }
      else
      {  /* elemental (unassembled) matrix */
         hbm->values = ucalloc(1+hbm->neltvl, sizeof(double));
         if (read_real_array(dsa, "values", hbm->valfmt, hbm->neltvl,
            hbm->values)) goto fail;
      }
      /* read right-hand sides */
      if (hbm->nrhs <= 0) goto done;
      if (hbm->rhstyp[0] == 'F')
      {  /* dense format */
         hbm->nrhsvl = hbm->nrow * hbm->nrhs;
         hbm->rhsval = ucalloc(1+hbm->nrhsvl, sizeof(double));
         if (read_real_array(dsa, "rhsval", hbm->rhsfmt, hbm->nrhsvl,
            hbm->rhsval)) goto fail;
      }
      else if (hbm->rhstyp[0] == 'M' && hbm->mxtype[2] == 'A')
      {  /* sparse format */
         /* read pointers */
         hbm->rhsptr = ucalloc(1+hbm->nrhs+1, sizeof(int));
         if (read_int_array(dsa, "rhsptr", hbm->ptrfmt, hbm->nrhs+1,
            hbm->rhsptr)) goto fail;
         /* read sparsity pattern */
         hbm->rhsind = ucalloc(1+hbm->nrhsix, sizeof(int));
         if (read_int_array(dsa, "rhsind", hbm->indfmt, hbm->nrhsix,
            hbm->rhsind)) goto fail;
         /* read values */
         hbm->rhsval = ucalloc(1+hbm->nrhsix, sizeof(double));
         if (read_real_array(dsa, "rhsval", hbm->rhsfmt, hbm->nrhsix,
            hbm->rhsval)) goto fail;
      }
      else if (hbm->rhstyp[0] == 'M' && hbm->mxtype[2] == 'E')
      {  /* elemental format */
         hbm->rhsval = ucalloc(1+hbm->nrhsvl, sizeof(double));
         if (read_real_array(dsa, "rhsval", hbm->rhsfmt, hbm->nrhsvl,
            hbm->rhsval)) goto fail;
      }
      else
      {  print("%s:%d: right-hand side type `%c' not recognised",
            dsa->fname, dsa->seqn, hbm->rhstyp[0]);
         goto fail;
      }
      /* read starting guesses */
      if (hbm->rhstyp[1] == 'G')
      {  hbm->nguess = hbm->nrow * hbm->nrhs;
         hbm->sguess = ucalloc(1+hbm->nguess, sizeof(double));
         if (read_real_array(dsa, "sguess", hbm->rhsfmt, hbm->nguess,
            hbm->sguess)) goto fail;
      }
      /* read solution vectors */
      if (hbm->rhstyp[2] == 'X')
      {  hbm->nexact = hbm->nrow * hbm->nrhs;
         hbm->xexact = ucalloc(1+hbm->nexact, sizeof(double));
         if (read_real_array(dsa, "xexact", hbm->rhsfmt, hbm->nexact,
            hbm->xexact)) goto fail;
      }
done: /* reading has been completed */
      print("hbm_read_mat: %d cards were read", dsa->seqn);
      ufclose(dsa->fp);
      return hbm;
fail: /* something wrong in Danish kingdom */
      if (hbm != NULL)
      {  if (hbm->colptr != NULL) ufree(hbm->colptr);
         if (hbm->rowind != NULL) ufree(hbm->rowind);
         if (hbm->rhsptr != NULL) ufree(hbm->rhsptr);
         if (hbm->rhsind != NULL) ufree(hbm->rhsind);
         if (hbm->values != NULL) ufree(hbm->values);
         if (hbm->rhsval != NULL) ufree(hbm->rhsval);
         if (hbm->sguess != NULL) ufree(hbm->sguess);
         if (hbm->xexact != NULL) ufree(hbm->xexact);
         ufree(hbm);
      }
      if (dsa->fp != NULL) ufclose(dsa->fp);
      return NULL;
}

/*----------------------------------------------------------------------
-- hbm_free_mat - free sparse matrix in Harwell-Boeing format.
--
-- *Synopsis*
--
-- #include "glphbm.h"
-- void hbm_free_mat(HBM *hbm);
--
-- *Description*
--
-- The hbm_free_mat routine frees all the memory allocated to the data
-- structure containing a sparse matrix in the Harwell-Boeing format. */

void hbm_free_mat(HBM *hbm)
{     if (hbm->colptr != NULL) ufree(hbm->colptr);
      if (hbm->rowind != NULL) ufree(hbm->rowind);
      if (hbm->rhsptr != NULL) ufree(hbm->rhsptr);
      if (hbm->rhsind != NULL) ufree(hbm->rhsind);
      if (hbm->values != NULL) ufree(hbm->values);
      if (hbm->rhsval != NULL) ufree(hbm->rhsval);
      if (hbm->sguess != NULL) ufree(hbm->sguess);
      if (hbm->xexact != NULL) ufree(hbm->xexact);
      ufree(hbm);
      return;
}

/* eof */
