lemon-project-template-glpk

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