PageRenderTime 56ms CodeModel.GetById 16ms app.highlight 35ms RepoModel.GetById 1ms app.codeStats 0ms

/farmR/src/glphbm.c

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