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