COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glphbm.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 18.9 KB
Line 
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
57struct 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
86static 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
133static 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
155static 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] != '(')
160fail: {  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;
190iter:    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
236static 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
275static 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
325HBM *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      }
471done: /* reading has been completed */
472      xprintf("hbm_read_mat: %d cards were read\n", dsa->seqn);
473      fclose(dsa->fp);
474      return hbm;
475fail: /* 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
506void 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 */
Note: See TracBrowser for help on using the repository browser.