lemon-project-template-glpk

comparison deps/glpk/src/glphbm.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:3d6a71cbdc5e
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, 2011 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 */