| 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 | |
|---|
| 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 */ |
|---|
