1 // (C) Copyright Jeremy Siek 2004
2 // Distributed under the Boost Software License, Version 1.0. (See
3 // accompanying file LICENSE_1_0.txt or copy at
4 // http://www.boost.org/LICENSE_1_0.txt)
5
6 /*
7 Fri Aug 15 16:29:47 EDT 1997
8
9 Harwell-Boeing File I/O in C
10 V. 1.0
11
12 National Institute of Standards and Technology, MD.
13 K.A. Remington
14
15 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
16 NOTICE
17
18 Permission to use, copy, modify, and distribute this software and
19 its documentation for any purpose and without fee is hereby granted
20 provided that the above copyright notice appear in all copies and
21 that both the copyright notice and this permission notice appear in
22 supporting documentation.
23
24 Neither the Author nor the Institution (National Institute of Standards
25 and Technology) make any representations about the suitability of this
26 software for any purpose. This software is provided "as is" without
27 expressed or implied warranty.
28 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
29
30 ---------------------
31 INTERFACE DESCRIPTION
32 ---------------------
33 ---------------
34 QUERY FUNCTIONS
35 ---------------
36
37 FUNCTION:
38
39 int readHB_info(const char *filename, int *M, int *N, int *nz,
40 char **Type, int *Nrhs)
41
42 DESCRIPTION:
43
44 The readHB_info function opens and reads the header information from
45 the specified Harwell-Boeing file, and reports back the number of rows
46 and columns in the stored matrix (M and N), the number of nonzeros in
47 the matrix (nz), the 3-character matrix type(Type), and the number of
48 right-hand-sides stored along with the matrix (Nrhs). This function
49 is designed to retrieve basic size information which can be used to
50 allocate arrays.
51
52 FUNCTION:
53
54 int readHB_header(FILE* in_file, char* Title, char* Key, char* Type,
55 int* Nrow, int* Ncol, int* Nnzero, int* Nrhs,
56 char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
57 int* Ptrcrd, int* Indcrd, int* Valcrd, int* Rhscrd,
58 char *Rhstype)
59
60 DESCRIPTION:
61
62 More detailed than the readHB_info function, readHB_header() reads from
63 the specified Harwell-Boeing file all of the header information.
64
65
66 ------------------------------
67 DOUBLE PRECISION I/O FUNCTIONS
68 ------------------------------
69 FUNCTION:
70
71 int readHB_newmat_double(const char *filename, int *M, int *N, *int nz,
72 int **colptr, int **rowind, double**val)
73
74 int readHB_mat_double(const char *filename, int *colptr, int *rowind,
75 double*val)
76
77
78 DESCRIPTION:
79
80 This function opens and reads the specified file, interpreting its
81 contents as a sparse matrix stored in the Harwell/Boeing standard
82 format. (See readHB_aux_double to read auxillary vectors.)
83 -- Values are interpreted as double precision numbers. --
84
85 The "mat" function uses _pre-allocated_ vectors to hold the index and
86 nonzero value information.
87
88 The "newmat" function allocates vectors to hold the index and nonzero
89 value information, and returns pointers to these vectors along with
90 matrix dimension and number of nonzeros.
91
92 FUNCTION:
93
94 int readHB_aux_double(const char* filename, const char AuxType, double b[])
95
96 int readHB_newaux_double(const char* filename, const char AuxType, double** b)
97
98 DESCRIPTION:
99
100 This function opens and reads from the specified file auxillary vector(s).
101 The char argument Auxtype determines which type of auxillary vector(s)
102 will be read (if present in the file).
103
104 AuxType = 'F' right-hand-side
105 AuxType = 'G' initial estimate (Guess)
106 AuxType = 'X' eXact solution
107
108 If Nrhs > 1, all of the Nrhs vectors of the given type are read and
109 stored in column-major order in the vector b.
110
111 The "newaux" function allocates a vector to hold the values retrieved.
112 The "mat" function uses a _pre-allocated_ vector to hold the values.
113
114 FUNCTION:
115
116 int writeHB_mat_double(const char* filename, int M, int N,
117 int nz, const int colptr[], const int rowind[],
118 const double val[], int Nrhs, const double rhs[],
119 const double guess[], const double exact[],
120 const char* Title, const char* Key, const char* Type,
121 char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
122 const char* Rhstype)
123
124 DESCRIPTION:
125
126 The writeHB_mat_double function opens the named file and writes the specified
127 matrix and optional auxillary vector(s) to that file in Harwell-Boeing
128 format. The format arguments (Ptrfmt,Indfmt,Valfmt, and Rhsfmt) are
129 character strings specifying "Fortran-style" output formats -- as they
130 would appear in a Harwell-Boeing file. They are used to produce output
131 which is as close as possible to what would be produced by Fortran code,
132 but note that "D" and "P" edit descriptors are not supported.
133 If NULL, the following defaults will be used:
134 Ptrfmt = Indfmt = "(8I10)"
135 Valfmt = Rhsfmt = "(4E20.13)"
136
137 -----------------------
138 CHARACTER I/O FUNCTIONS
139 -----------------------
140 FUNCTION:
141
142 int readHB_mat_char(const char* filename, int colptr[], int rowind[],
143 char val[], char* Valfmt)
144 int readHB_newmat_char(const char* filename, int* M, int* N, int* nonzeros,
145 int** colptr, int** rowind, char** val, char** Valfmt)
146
147 DESCRIPTION:
148
149 This function opens and reads the specified file, interpreting its
150 contents as a sparse matrix stored in the Harwell/Boeing standard
151 format. (See readHB_aux_char to read auxillary vectors.)
152 -- Values are interpreted as char strings. --
153 (Used to translate exact values from the file into a new storage format.)
154
155 The "mat" function uses _pre-allocated_ arrays to hold the index and
156 nonzero value information.
157
158 The "newmat" function allocates char arrays to hold the index
159 and nonzero value information, and returns pointers to these arrays
160 along with matrix dimension and number of nonzeros.
161
162 FUNCTION:
163
164 int readHB_aux_char(const char* filename, const char AuxType, char b[])
165 int readHB_newaux_char(const char* filename, const char AuxType, char** b,
166 char** Rhsfmt)
167
168 DESCRIPTION:
169
170 This function opens and reads from the specified file auxillary vector(s).
171 The char argument Auxtype determines which type of auxillary vector(s)
172 will be read (if present in the file).
173
174 AuxType = 'F' right-hand-side
175 AuxType = 'G' initial estimate (Guess)
176 AuxType = 'X' eXact solution
177
178 If Nrhs > 1, all of the Nrhs vectors of the given type are read and
179 stored in column-major order in the vector b.
180
181 The "newaux" function allocates a character array to hold the values
182 retrieved.
183 The "mat" function uses a _pre-allocated_ array to hold the values.
184
185 FUNCTION:
186
187 int writeHB_mat_char(const char* filename, int M, int N,
188 int nz, const int colptr[], const int rowind[],
189 const char val[], int Nrhs, const char rhs[],
190 const char guess[], const char exact[],
191 const char* Title, const char* Key, const char* Type,
192 char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
193 const char* Rhstype)
194
195 DESCRIPTION:
196
197 The writeHB_mat_char function opens the named file and writes the specified
198 matrix and optional auxillary vector(s) to that file in Harwell-Boeing
199 format. The format arguments (Ptrfmt,Indfmt,Valfmt, and Rhsfmt) are
200 character strings specifying "Fortran-style" output formats -- as they
201 would appear in a Harwell-Boeing file. Valfmt and Rhsfmt must accurately
202 represent the character representation of the values stored in val[]
203 and rhs[].
204
205 If NULL, the following defaults will be used for the integer vectors:
206 Ptrfmt = Indfmt = "(8I10)"
207 Valfmt = Rhsfmt = "(4E20.13)"
208
209
210 */
211
212 /*---------------------------------------------------------------------*/
213 /* If zero-based indexing is desired, _SP_base should be set to 0 */
214 /* This will cause indices read from H-B files to be decremented by 1 */
215 /* and indices written to H-B files to be incremented by 1 */
216 /* <<< Standard usage is _SP_base = 1 >>> */
217 #ifndef _SP_base
218 #define _SP_base 1
219 #endif
220 /*---------------------------------------------------------------------*/
221
222 #include "iohb.h"
223 #include<stdio.h>
224 #include<stdlib.h>
225 #include<string.h>
226 #include<math.h>
227
228 char* substr(const char* S, const int pos, const int len);
229 void upcase(char* S);
230 void IOHBTerminate(char* message);
231
readHB_info(const char * filename,int * M,int * N,int * nz,char ** Type,int * Nrhs)232 int readHB_info(const char* filename, int* M, int* N, int* nz, char** Type,
233 int* Nrhs)
234 {
235 /****************************************************************************/
236 /* The readHB_info function opens and reads the header information from */
237 /* the specified Harwell-Boeing file, and reports back the number of rows */
238 /* and columns in the stored matrix (M and N), the number of nonzeros in */
239 /* the matrix (nz), and the number of right-hand-sides stored along with */
240 /* the matrix (Nrhs). */
241 /* */
242 /* For a description of the Harwell Boeing standard, see: */
243 /* Duff, et al., ACM TOMS Vol.15, No.1, March 1989 */
244 /* */
245 /* ---------- */
246 /* **CAVEAT** */
247 /* ---------- */
248 /* ** If the input file does not adhere to the H/B format, the ** */
249 /* ** results will be unpredictable. ** */
250 /* */
251 /****************************************************************************/
252 FILE *in_file;
253 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
254 int Nrow, Ncol, Nnzero;
255 char* mat_type;
256 char Title[73], Key[9], Rhstype[4];
257 char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
258
259 mat_type = *Type;
260 if ( mat_type == NULL ) IOHBTerminate("Insufficient memory for mat_typen");
261
262 if ( (in_file = fopen( filename, "r")) == NULL ) {
263 fprintf(stderr,"Error: Cannot open file: %s\n",filename);
264 return 0;
265 }
266
267 readHB_header(in_file, Title, Key, mat_type, &Nrow, &Ncol, &Nnzero, Nrhs,
268 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
269 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
270 fclose(in_file);
271 *Type = mat_type;
272 *(*Type+3) = (char) NULL;
273 *M = Nrow;
274 *N = Ncol;
275 *nz = Nnzero;
276 if (Rhscrd == 0) {*Nrhs = 0;}
277
278 /* In verbose mode, print some of the header information: */
279 /*
280 if (verbose == 1)
281 {
282 printf("Reading from Harwell-Boeing file %s (verbose on)...\n",filename);
283 printf(" Title: %s\n",Title);
284 printf(" Key: %s\n",Key);
285 printf(" The stored matrix is %i by %i with %i nonzeros.\n",
286 *M, *N, *nz );
287 printf(" %i right-hand--side(s) stored.\n",*Nrhs);
288 }
289 */
290
291 return 1;
292
293 }
294
295
296
readHB_header(FILE * in_file,char * Title,char * Key,char * Type,int * Nrow,int * Ncol,int * Nnzero,int * Nrhs,char * Ptrfmt,char * Indfmt,char * Valfmt,char * Rhsfmt,int * Ptrcrd,int * Indcrd,int * Valcrd,int * Rhscrd,char * Rhstype)297 int readHB_header(FILE* in_file, char* Title, char* Key, char* Type,
298 int* Nrow, int* Ncol, int* Nnzero, int* Nrhs,
299 char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
300 int* Ptrcrd, int* Indcrd, int* Valcrd, int* Rhscrd,
301 char *Rhstype)
302 {
303 /*************************************************************************/
304 /* Read header information from the named H/B file... */
305 /*************************************************************************/
306 int Totcrd,Neltvl,Nrhsix;
307 char line[BUFSIZ];
308
309 /* First line: */
310 fgets(line, BUFSIZ, in_file);
311 if ( sscanf(line,"%*s") < 0 )
312 IOHBTerminate("iohb.c: Null (or blank) first line of HB file.\n");
313 (void) sscanf(line, "%72c%8[^\n]", Title, Key);
314 *(Key+8) = (char) NULL;
315 *(Title+72) = (char) NULL;
316
317 /* Second line: */
318 fgets(line, BUFSIZ, in_file);
319 if ( sscanf(line,"%*s") < 0 )
320 IOHBTerminate("iohb.c: Null (or blank) second line of HB file.\n");
321 if ( sscanf(line,"%i",&Totcrd) != 1) Totcrd = 0;
322 if ( sscanf(line,"%*i%i",Ptrcrd) != 1) *Ptrcrd = 0;
323 if ( sscanf(line,"%*i%*i%i",Indcrd) != 1) *Indcrd = 0;
324 if ( sscanf(line,"%*i%*i%*i%i",Valcrd) != 1) *Valcrd = 0;
325 if ( sscanf(line,"%*i%*i%*i%*i%i",Rhscrd) != 1) *Rhscrd = 0;
326
327 /* Third line: */
328 fgets(line, BUFSIZ, in_file);
329 if ( sscanf(line,"%*s") < 0 )
330 IOHBTerminate("iohb.c: Null (or blank) third line of HB file.\n");
331 if ( sscanf(line, "%3c", Type) != 1)
332 IOHBTerminate("iohb.c: Invalid Type info, line 3 of Harwell-Boeing file.\n");
333 upcase(Type);
334 if ( sscanf(line,"%*3c%i",Nrow) != 1) *Nrow = 0 ;
335 if ( sscanf(line,"%*3c%*i%i",Ncol) != 1) *Ncol = 0 ;
336 if ( sscanf(line,"%*3c%*i%*i%i",Nnzero) != 1) *Nnzero = 0 ;
337 if ( sscanf(line,"%*3c%*i%*i%*i%i",&Neltvl) != 1) Neltvl = 0 ;
338
339 /* Fourth line: */
340 fgets(line, BUFSIZ, in_file);
341 if ( sscanf(line,"%*s") < 0 )
342 IOHBTerminate("iohb.c: Null (or blank) fourth line of HB file.\n");
343 if ( sscanf(line, "%16c",Ptrfmt) != 1)
344 IOHBTerminate("iohb.c: Invalid format info, line 4 of Harwell-Boeing file.\n");
345 if ( sscanf(line, "%*16c%16c",Indfmt) != 1)
346 IOHBTerminate("iohb.c: Invalid format info, line 4 of Harwell-Boeing file.\n");
347 if ( sscanf(line, "%*16c%*16c%20c",Valfmt) != 1)
348 IOHBTerminate("iohb.c: Invalid format info, line 4 of Harwell-Boeing file.\n");
349 sscanf(line, "%*16c%*16c%*20c%20c",Rhsfmt);
350 *(Ptrfmt+16) = (char) NULL;
351 *(Indfmt+16) = (char) NULL;
352 *(Valfmt+20) = (char) NULL;
353 *(Rhsfmt+20) = (char) NULL;
354
355 /* (Optional) Fifth line: */
356 if (*Rhscrd != 0 )
357 {
358 fgets(line, BUFSIZ, in_file);
359 if ( sscanf(line,"%*s") < 0 )
360 IOHBTerminate("iohb.c: Null (or blank) fifth line of HB file.\n");
361 if ( sscanf(line, "%3c", Rhstype) != 1)
362 IOHBTerminate("iohb.c: Invalid RHS type information, line 5 of Harwell-Boeing file.\n");
363 if ( sscanf(line, "%*3c%i", Nrhs) != 1) *Nrhs = 0;
364 if ( sscanf(line, "%*3c%*i%i", &Nrhsix) != 1) Nrhsix = 0;
365 }
366 return 1;
367 }
368
369
readHB_mat_double(const char * filename,int colptr[],int rowind[],double val[])370 int readHB_mat_double(const char* filename, int colptr[], int rowind[],
371 double val[])
372 {
373 /****************************************************************************/
374 /* This function opens and reads the specified file, interpreting its */
375 /* contents as a sparse matrix stored in the Harwell/Boeing standard */
376 /* format and creating compressed column storage scheme vectors to hold */
377 /* the index and nonzero value information. */
378 /* */
379 /* ---------- */
380 /* **CAVEAT** */
381 /* ---------- */
382 /* Parsing real formats from Fortran is tricky, and this file reader */
383 /* does not claim to be foolproof. It has been tested for cases when */
384 /* the real values are printed consistently and evenly spaced on each */
385 /* line, with Fixed (F), and Exponential (E or D) formats. */
386 /* */
387 /* ** If the input file does not adhere to the H/B format, the ** */
388 /* ** results will be unpredictable. ** */
389 /* */
390 /****************************************************************************/
391 FILE *in_file;
392 int i,j,ind,col,offset,count,last,Nrhs;
393 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
394 int Nrow, Ncol, Nnzero, Nentries;
395 int Ptrperline, Ptrwidth, Indperline, Indwidth;
396 int Valperline, Valwidth, Valprec;
397 int Valflag; /* Indicates 'E','D', or 'F' float format */
398 char* ThisElement;
399 char Title[73], Key[8], Type[4], Rhstype[4];
400 char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
401 char line[BUFSIZ];
402
403 if ( (in_file = fopen( filename, "r")) == NULL ) {
404 fprintf(stderr,"Error: Cannot open file: %s\n",filename);
405 return 0;
406 }
407
408 readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
409 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
410 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
411
412 /* Parse the array input formats from Line 3 of HB file */
413 ParseIfmt(Ptrfmt,&Ptrperline,&Ptrwidth);
414 ParseIfmt(Indfmt,&Indperline,&Indwidth);
415 if ( Type[0] != 'P' ) { /* Skip if pattern only */
416 ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
417 }
418
419 /* Read column pointer array: */
420
421 offset = 1-_SP_base; /* if base 0 storage is declared (via macro definition), */
422 /* then storage entries are offset by 1 */
423
424 ThisElement = (char *) malloc(Ptrwidth+1);
425 if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
426 *(ThisElement+Ptrwidth) = (char) NULL;
427 count=0;
428 for (i=0;i<Ptrcrd;i++)
429 {
430 fgets(line, BUFSIZ, in_file);
431 if ( sscanf(line,"%*s") < 0 )
432 IOHBTerminate("iohb.c: Null (or blank) line in pointer data region of HB file.\n");
433 col = 0;
434 for (ind = 0;ind<Ptrperline;ind++)
435 {
436 if (count > Ncol) break;
437 strncpy(ThisElement,line+col,Ptrwidth);
438 /* ThisElement = substr(line,col,Ptrwidth); */
439 colptr[count] = atoi(ThisElement)-offset;
440 count++; col += Ptrwidth;
441 }
442 }
443 free(ThisElement);
444
445 /* Read row index array: */
446
447 ThisElement = (char *) malloc(Indwidth+1);
448 if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
449 *(ThisElement+Indwidth) = (char) NULL;
450 count = 0;
451 for (i=0;i<Indcrd;i++)
452 {
453 fgets(line, BUFSIZ, in_file);
454 if ( sscanf(line,"%*s") < 0 )
455 IOHBTerminate("iohb.c: Null (or blank) line in index data region of HB file.\n");
456 col = 0;
457 for (ind = 0;ind<Indperline;ind++)
458 {
459 if (count == Nnzero) break;
460 strncpy(ThisElement,line+col,Indwidth);
461 /* ThisElement = substr(line,col,Indwidth); */
462 rowind[count] = atoi(ThisElement)-offset;
463 count++; col += Indwidth;
464 }
465 }
466 free(ThisElement);
467
468 /* Read array of values: */
469
470 if ( Type[0] != 'P' ) { /* Skip if pattern only */
471
472 if ( Type[0] == 'C' ) Nentries = 2*Nnzero;
473 else Nentries = Nnzero;
474
475 ThisElement = (char *) malloc(Valwidth+1);
476 if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
477 *(ThisElement+Valwidth) = (char) NULL;
478 count = 0;
479 for (i=0;i<Valcrd;i++)
480 {
481 fgets(line, BUFSIZ, in_file);
482 if ( sscanf(line,"%*s") < 0 )
483 IOHBTerminate("iohb.c: Null (or blank) line in value data region of HB file.\n");
484 if (Valflag == 'D') {
485 while( strchr(line,'D') ) *strchr(line,'D') = 'E';
486 /* *strchr(Valfmt,'D') = 'E'; */
487 }
488 col = 0;
489 for (ind = 0;ind<Valperline;ind++)
490 {
491 if (count == Nentries) break;
492 strncpy(ThisElement,line+col,Valwidth);
493 /*ThisElement = substr(line,col,Valwidth);*/
494 if ( Valflag != 'F' && strchr(ThisElement,'E') == NULL ) {
495 /* insert a char prefix for exp */
496 last = strlen(ThisElement);
497 for (j=last+1;j>=0;j--) {
498 ThisElement[j] = ThisElement[j-1];
499 if ( ThisElement[j] == '+' || ThisElement[j] == '-' ) {
500 ThisElement[j-1] = Valflag;
501 break;
502 }
503 }
504 }
505 val[count] = atof(ThisElement);
506 count++; col += Valwidth;
507 }
508 }
509 free(ThisElement);
510 }
511
512 fclose(in_file);
513 return 1;
514 }
515
readHB_newmat_double(const char * filename,int * M,int * N,int * nonzeros,int ** colptr,int ** rowind,double ** val)516 int readHB_newmat_double(const char* filename, int* M, int* N, int* nonzeros,
517 int** colptr, int** rowind, double** val)
518 {
519 int Nrhs;
520 char *Type;
521
522 readHB_info(filename, M, N, nonzeros, &Type, &Nrhs);
523
524 *colptr = (int *)malloc((*N+1)*sizeof(int));
525 if ( *colptr == NULL ) IOHBTerminate("Insufficient memory for colptr.\n");
526 *rowind = (int *)malloc(*nonzeros*sizeof(int));
527 if ( *rowind == NULL ) IOHBTerminate("Insufficient memory for rowind.\n");
528 if ( Type[0] == 'C' ) {
529 /*
530 fprintf(stderr, "Warning: Reading complex data from HB file %s.\n",filename);
531 fprintf(stderr, " Real and imaginary parts will be interlaced in val[].\n");
532 */
533 /* Malloc enough space for real AND imaginary parts of val[] */
534 *val = (double *)malloc(*nonzeros*sizeof(double)*2);
535 if ( *val == NULL ) IOHBTerminate("Insufficient memory for val.\n");
536 } else {
537 if ( Type[0] != 'P' ) {
538 /* Malloc enough space for real array val[] */
539 *val = (double *)malloc(*nonzeros*sizeof(double));
540 if ( *val == NULL ) IOHBTerminate("Insufficient memory for val.\n");
541 }
542 } /* No val[] space needed if pattern only */
543 return readHB_mat_double(filename, *colptr, *rowind, *val);
544
545 }
546
readHB_aux_double(const char * filename,const char AuxType,double b[])547 int readHB_aux_double(const char* filename, const char AuxType, double b[])
548 {
549 /****************************************************************************/
550 /* This function opens and reads the specified file, placing auxillary */
551 /* vector(s) of the given type (if available) in b. */
552 /* Return value is the number of vectors successfully read. */
553 /* */
554 /* AuxType = 'F' full right-hand-side vector(s) */
555 /* AuxType = 'G' initial Guess vector(s) */
556 /* AuxType = 'X' eXact solution vector(s) */
557 /* */
558 /* ---------- */
559 /* **CAVEAT** */
560 /* ---------- */
561 /* Parsing real formats from Fortran is tricky, and this file reader */
562 /* does not claim to be foolproof. It has been tested for cases when */
563 /* the real values are printed consistently and evenly spaced on each */
564 /* line, with Fixed (F), and Exponential (E or D) formats. */
565 /* */
566 /* ** If the input file does not adhere to the H/B format, the ** */
567 /* ** results will be unpredictable. ** */
568 /* */
569 /****************************************************************************/
570 FILE *in_file;
571 int i,j,n,maxcol,start,stride,col,last,linel;
572 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
573 int Nrow, Ncol, Nnzero, Nentries;
574 int Nrhs, nvecs, rhsi;
575 int Rhsperline, Rhswidth, Rhsprec;
576 int Rhsflag;
577 char *ThisElement;
578 char Title[73], Key[9], Type[4], Rhstype[4];
579 char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
580 char line[BUFSIZ];
581
582 if ((in_file = fopen( filename, "r")) == NULL) {
583 fprintf(stderr,"Error: Cannot open file: %s\n",filename);
584 return 0;
585 }
586
587 readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
588 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
589 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
590
591 if (Nrhs <= 0)
592 {
593 fprintf(stderr, "Warn: Attempt to read auxillary vector(s) when none are present.\n");
594 return 0;
595 }
596 if (Rhstype[0] != 'F' )
597 {
598 fprintf(stderr,"Warn: Attempt to read auxillary vector(s) which are not stored in Full form.\n");
599 fprintf(stderr," Rhs must be specified as full. \n");
600 return 0;
601 }
602
603 /* If reading complex data, allow for interleaved real and imaginary values. */
604 if ( Type[0] == 'C' ) {
605 Nentries = 2*Nrow;
606 } else {
607 Nentries = Nrow;
608 }
609
610 nvecs = 1;
611
612 if ( Rhstype[1] == 'G' ) nvecs++;
613 if ( Rhstype[2] == 'X' ) nvecs++;
614
615 if ( AuxType == 'G' && Rhstype[1] != 'G' ) {
616 fprintf(stderr, "Warn: Attempt to read auxillary Guess vector(s) when none are present.\n");
617 return 0;
618 }
619 if ( AuxType == 'X' && Rhstype[2] != 'X' ) {
620 fprintf(stderr, "Warn: Attempt to read auxillary eXact solution vector(s) when none are present.\n");
621 return 0;
622 }
623
624 ParseRfmt(Rhsfmt, &Rhsperline, &Rhswidth, &Rhsprec,&Rhsflag);
625 maxcol = Rhsperline*Rhswidth;
626
627 /* Lines to skip before starting to read RHS values... */
628 n = Ptrcrd + Indcrd + Valcrd;
629
630 for (i = 0; i < n; i++)
631 fgets(line, BUFSIZ, in_file);
632
633 /* start - number of initial aux vector entries to skip */
634 /* to reach first vector requested */
635 /* stride - number of aux vector entries to skip between */
636 /* requested vectors */
637 if ( AuxType == 'F' ) start = 0;
638 else if ( AuxType == 'G' ) start = Nentries;
639 else start = (nvecs-1)*Nentries;
640 stride = (nvecs-1)*Nentries;
641
642 fgets(line, BUFSIZ, in_file);
643 linel= strchr(line,'\n')-line;
644 col = 0;
645 /* Skip to initial offset */
646
647 for (i=0;i<start;i++) {
648 if ( col >= ( maxcol<linel?maxcol:linel ) ) {
649 fgets(line, BUFSIZ, in_file);
650 linel= strchr(line,'\n')-line;
651 col = 0;
652 }
653 col += Rhswidth;
654 }
655 if (Rhsflag == 'D') {
656 while( strchr(line,'D') ) *strchr(line,'D') = 'E';
657 }
658
659 /* Read a vector of desired type, then skip to next */
660 /* repeating to fill Nrhs vectors */
661
662 ThisElement = (char *) malloc(Rhswidth+1);
663 if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
664 *(ThisElement+Rhswidth) = (char) NULL;
665 for (rhsi=0;rhsi<Nrhs;rhsi++) {
666
667 for (i=0;i<Nentries;i++) {
668 if ( col >= ( maxcol<linel?maxcol:linel ) ) {
669 fgets(line, BUFSIZ, in_file);
670 linel= strchr(line,'\n')-line;
671 if (Rhsflag == 'D') {
672 while( strchr(line,'D') ) *strchr(line,'D') = 'E';
673 }
674 col = 0;
675 }
676 strncpy(ThisElement,line+col,Rhswidth);
677 /*ThisElement = substr(line, col, Rhswidth);*/
678 if ( Rhsflag != 'F' && strchr(ThisElement,'E') == NULL ) {
679 /* insert a char prefix for exp */
680 last = strlen(ThisElement);
681 for (j=last+1;j>=0;j--) {
682 ThisElement[j] = ThisElement[j-1];
683 if ( ThisElement[j] == '+' || ThisElement[j] == '-' ) {
684 ThisElement[j-1] = Rhsflag;
685 break;
686 }
687 }
688 }
689 b[i] = atof(ThisElement);
690 col += Rhswidth;
691 }
692
693 /* Skip any interleaved Guess/eXact vectors */
694
695 for (i=0;i<stride;i++) {
696 if ( col >= ( maxcol<linel?maxcol:linel ) ) {
697 fgets(line, BUFSIZ, in_file);
698 linel= strchr(line,'\n')-line;
699 col = 0;
700 }
701 col += Rhswidth;
702 }
703
704 }
705 free(ThisElement);
706
707
708 fclose(in_file);
709 return Nrhs;
710 }
711
readHB_newaux_double(const char * filename,const char AuxType,double ** b)712 int readHB_newaux_double(const char* filename, const char AuxType, double** b)
713 {
714 int Nrhs,M,N,nonzeros;
715 char *Type;
716
717 readHB_info(filename, &M, &N, &nonzeros, &Type, &Nrhs);
718 if ( Nrhs <= 0 ) {
719 fprintf(stderr,"Warn: Requested read of aux vector(s) when none are present.\n");
720 return 0;
721 } else {
722 if ( Type[0] == 'C' ) {
723 fprintf(stderr, "Warning: Reading complex aux vector(s) from HB file %s.",filename);
724 fprintf(stderr, " Real and imaginary parts will be interlaced in b[].");
725 *b = (double *)malloc(M*Nrhs*sizeof(double)*2);
726 if ( *b == NULL ) IOHBTerminate("Insufficient memory for rhs.\n");
727 return readHB_aux_double(filename, AuxType, *b);
728 } else {
729 *b = (double *)malloc(M*Nrhs*sizeof(double));
730 if ( *b == NULL ) IOHBTerminate("Insufficient memory for rhs.\n");
731 return readHB_aux_double(filename, AuxType, *b);
732 }
733 }
734 }
735
writeHB_mat_double(const char * filename,int M,int N,int nz,const int colptr[],const int rowind[],const double val[],int Nrhs,const double rhs[],const double guess[],const double exact[],const char * Title,const char * Key,const char * Type,char * Ptrfmt,char * Indfmt,char * Valfmt,char * Rhsfmt,const char * Rhstype)736 int writeHB_mat_double(const char* filename, int M, int N,
737 int nz, const int colptr[], const int rowind[],
738 const double val[], int Nrhs, const double rhs[],
739 const double guess[], const double exact[],
740 const char* Title, const char* Key, const char* Type,
741 char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
742 const char* Rhstype)
743 {
744 /****************************************************************************/
745 /* The writeHB function opens the named file and writes the specified */
746 /* matrix and optional right-hand-side(s) to that file in Harwell-Boeing */
747 /* format. */
748 /* */
749 /* For a description of the Harwell Boeing standard, see: */
750 /* Duff, et al., ACM TOMS Vol.15, No.1, March 1989 */
751 /* */
752 /****************************************************************************/
753 FILE *out_file;
754 int i,j,entry,offset,acount,linemod;
755 int totcrd, ptrcrd, indcrd, valcrd, rhscrd;
756 int nvalentries, nrhsentries;
757 int Ptrperline, Ptrwidth, Indperline, Indwidth;
758 int Rhsperline, Rhswidth, Rhsprec;
759 int Rhsflag;
760 int Valperline, Valwidth, Valprec;
761 int Valflag; /* Indicates 'E','D', or 'F' float format */
762 char pformat[16],iformat[16],vformat[19],rformat[19];
763
764 if ( Type[0] == 'C' ) {
765 nvalentries = 2*nz;
766 nrhsentries = 2*M;
767 } else {
768 nvalentries = nz;
769 nrhsentries = M;
770 }
771
772 if ( filename != NULL ) {
773 if ( (out_file = fopen( filename, "w")) == NULL ) {
774 fprintf(stderr,"Error: Cannot open file: %s\n",filename);
775 return 0;
776 }
777 } else out_file = stdout;
778
779 if ( Ptrfmt == NULL ) Ptrfmt = "(8I10)";
780 ParseIfmt(Ptrfmt,&Ptrperline,&Ptrwidth);
781 sprintf(pformat,"%%%dd",Ptrwidth);
782 ptrcrd = (N+1)/Ptrperline;
783 if ( (N+1)%Ptrperline != 0) ptrcrd++;
784
785 if ( Indfmt == NULL ) Indfmt = Ptrfmt;
786 ParseIfmt(Indfmt,&Indperline,&Indwidth);
787 sprintf(iformat,"%%%dd",Indwidth);
788 indcrd = nz/Indperline;
789 if ( nz%Indperline != 0) indcrd++;
790
791 if ( Type[0] != 'P' ) { /* Skip if pattern only */
792 if ( Valfmt == NULL ) Valfmt = "(4E20.13)";
793 ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
794 if (Valflag == 'D') *strchr(Valfmt,'D') = 'E';
795 if (Valflag == 'F')
796 sprintf(vformat,"%% %d.%df",Valwidth,Valprec);
797 else
798 sprintf(vformat,"%% %d.%dE",Valwidth,Valprec);
799 valcrd = nvalentries/Valperline;
800 if ( nvalentries%Valperline != 0) valcrd++;
801 } else valcrd = 0;
802
803 if ( Nrhs > 0 ) {
804 if ( Rhsfmt == NULL ) Rhsfmt = Valfmt;
805 ParseRfmt(Rhsfmt,&Rhsperline,&Rhswidth,&Rhsprec, &Rhsflag);
806 if (Rhsflag == 'F')
807 sprintf(rformat,"%% %d.%df",Rhswidth,Rhsprec);
808 else
809 sprintf(rformat,"%% %d.%dE",Rhswidth,Rhsprec);
810 if (Rhsflag == 'D') *strchr(Rhsfmt,'D') = 'E';
811 rhscrd = nrhsentries/Rhsperline;
812 if ( nrhsentries%Rhsperline != 0) rhscrd++;
813 if ( Rhstype[1] == 'G' ) rhscrd+=rhscrd;
814 if ( Rhstype[2] == 'X' ) rhscrd+=rhscrd;
815 rhscrd*=Nrhs;
816 } else rhscrd = 0;
817
818 totcrd = 4+ptrcrd+indcrd+valcrd+rhscrd;
819
820
821 /* Print header information: */
822
823 fprintf(out_file,"%-72s%-8s\n%14d%14d%14d%14d%14d\n",Title, Key, totcrd,
824 ptrcrd, indcrd, valcrd, rhscrd);
825 fprintf(out_file,"%3s%11s%14d%14d%14d\n",Type," ", M, N, nz);
826 fprintf(out_file,"%-16s%-16s%-20s", Ptrfmt, Indfmt, Valfmt);
827 if ( Nrhs != 0 ) {
828 /* Print Rhsfmt on fourth line and */
829 /* optional fifth header line for auxillary vector information: */
830 fprintf(out_file,"%-20s\n%-14s%d\n",Rhsfmt,Rhstype,Nrhs);
831 } else fprintf(out_file,"\n");
832
833 offset = 1-_SP_base; /* if base 0 storage is declared (via macro definition), */
834 /* then storage entries are offset by 1 */
835
836 /* Print column pointers: */
837 for (i=0;i<N+1;i++)
838 {
839 entry = colptr[i]+offset;
840 fprintf(out_file,pformat,entry);
841 if ( (i+1)%Ptrperline == 0 ) fprintf(out_file,"\n");
842 }
843
844 if ( (N+1) % Ptrperline != 0 ) fprintf(out_file,"\n");
845
846 /* Print row indices: */
847 for (i=0;i<nz;i++)
848 {
849 entry = rowind[i]+offset;
850 fprintf(out_file,iformat,entry);
851 if ( (i+1)%Indperline == 0 ) fprintf(out_file,"\n");
852 }
853
854 if ( nz % Indperline != 0 ) fprintf(out_file,"\n");
855
856 /* Print values: */
857
858 if ( Type[0] != 'P' ) { /* Skip if pattern only */
859
860 for (i=0;i<nvalentries;i++)
861 {
862 fprintf(out_file,vformat,val[i]);
863 if ( (i+1)%Valperline == 0 ) fprintf(out_file,"\n");
864 }
865
866 if ( nvalentries % Valperline != 0 ) fprintf(out_file,"\n");
867
868 /* If available, print right hand sides,
869 guess vectors and exact solution vectors: */
870 acount = 1;
871 linemod = 0;
872 if ( Nrhs > 0 ) {
873 for (i=0;i<Nrhs;i++)
874 {
875 for ( j=0;j<nrhsentries;j++ ) {
876 fprintf(out_file,rformat,rhs[j]);
877 if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
878 }
879 if ( acount%Rhsperline != linemod ) {
880 fprintf(out_file,"\n");
881 linemod = (acount-1)%Rhsperline;
882 }
883 rhs += nrhsentries;
884 if ( Rhstype[1] == 'G' ) {
885 for ( j=0;j<nrhsentries;j++ ) {
886 fprintf(out_file,rformat,guess[j]);
887 if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
888 }
889 if ( acount%Rhsperline != linemod ) {
890 fprintf(out_file,"\n");
891 linemod = (acount-1)%Rhsperline;
892 }
893 guess += nrhsentries;
894 }
895 if ( Rhstype[2] == 'X' ) {
896 for ( j=0;j<nrhsentries;j++ ) {
897 fprintf(out_file,rformat,exact[j]);
898 if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
899 }
900 if ( acount%Rhsperline != linemod ) {
901 fprintf(out_file,"\n");
902 linemod = (acount-1)%Rhsperline;
903 }
904 exact += nrhsentries;
905 }
906 }
907 }
908
909 }
910
911 if ( fclose(out_file) != 0){
912 fprintf(stderr,"Error closing file in writeHB_mat_double().\n");
913 return 0;
914 } else return 1;
915
916 }
917
readHB_mat_char(const char * filename,int colptr[],int rowind[],char val[],char * Valfmt)918 int readHB_mat_char(const char* filename, int colptr[], int rowind[],
919 char val[], char* Valfmt)
920 {
921 /****************************************************************************/
922 /* This function opens and reads the specified file, interpreting its */
923 /* contents as a sparse matrix stored in the Harwell/Boeing standard */
924 /* format and creating compressed column storage scheme vectors to hold */
925 /* the index and nonzero value information. */
926 /* */
927 /* ---------- */
928 /* **CAVEAT** */
929 /* ---------- */
930 /* Parsing real formats from Fortran is tricky, and this file reader */
931 /* does not claim to be foolproof. It has been tested for cases when */
932 /* the real values are printed consistently and evenly spaced on each */
933 /* line, with Fixed (F), and Exponential (E or D) formats. */
934 /* */
935 /* ** If the input file does not adhere to the H/B format, the ** */
936 /* ** results will be unpredictable. ** */
937 /* */
938 /****************************************************************************/
939 FILE *in_file;
940 int i,j,ind,col,offset,count,last;
941 int Nrow,Ncol,Nnzero,Nentries,Nrhs;
942 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
943 int Ptrperline, Ptrwidth, Indperline, Indwidth;
944 int Valperline, Valwidth, Valprec;
945 int Valflag; /* Indicates 'E','D', or 'F' float format */
946 char* ThisElement;
947 char line[BUFSIZ];
948 char Title[73], Key[8], Type[4], Rhstype[4];
949 char Ptrfmt[17], Indfmt[17], Rhsfmt[21];
950
951 if ( (in_file = fopen( filename, "r")) == NULL ) {
952 fprintf(stderr,"Error: Cannot open file: %s\n",filename);
953 return 0;
954 }
955
956 readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
957 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
958 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
959
960 /* Parse the array input formats from Line 3 of HB file */
961 ParseIfmt(Ptrfmt,&Ptrperline,&Ptrwidth);
962 ParseIfmt(Indfmt,&Indperline,&Indwidth);
963 if ( Type[0] != 'P' ) { /* Skip if pattern only */
964 ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
965 if (Valflag == 'D') {
966 *strchr(Valfmt,'D') = 'E';
967 }
968 }
969
970 /* Read column pointer array: */
971
972 offset = 1-_SP_base; /* if base 0 storage is declared (via macro definition), */
973 /* then storage entries are offset by 1 */
974
975 ThisElement = (char *) malloc(Ptrwidth+1);
976 if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
977 *(ThisElement+Ptrwidth) = (char) NULL;
978 count=0;
979 for (i=0;i<Ptrcrd;i++)
980 {
981 fgets(line, BUFSIZ, in_file);
982 if ( sscanf(line,"%*s") < 0 )
983 IOHBTerminate("iohb.c: Null (or blank) line in pointer data region of HB file.\n");
984 col = 0;
985 for (ind = 0;ind<Ptrperline;ind++)
986 {
987 if (count > Ncol) break;
988 strncpy(ThisElement,line+col,Ptrwidth);
989 /*ThisElement = substr(line,col,Ptrwidth);*/
990 colptr[count] = atoi(ThisElement)-offset;
991 count++; col += Ptrwidth;
992 }
993 }
994 free(ThisElement);
995
996 /* Read row index array: */
997
998 ThisElement = (char *) malloc(Indwidth+1);
999 if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
1000 *(ThisElement+Indwidth) = (char) NULL;
1001 count = 0;
1002 for (i=0;i<Indcrd;i++)
1003 {
1004 fgets(line, BUFSIZ, in_file);
1005 if ( sscanf(line,"%*s") < 0 )
1006 IOHBTerminate("iohb.c: Null (or blank) line in index data region of HB file.\n");
1007 col = 0;
1008 for (ind = 0;ind<Indperline;ind++)
1009 {
1010 if (count == Nnzero) break;
1011 strncpy(ThisElement,line+col,Indwidth);
1012 /*ThisElement = substr(line,col,Indwidth);*/
1013 rowind[count] = atoi(ThisElement)-offset;
1014 count++; col += Indwidth;
1015 }
1016 }
1017 free(ThisElement);
1018
1019 /* Read array of values: AS CHARACTERS*/
1020
1021 if ( Type[0] != 'P' ) { /* Skip if pattern only */
1022
1023 if ( Type[0] == 'C' ) Nentries = 2*Nnzero;
1024 else Nentries = Nnzero;
1025
1026 ThisElement = (char *) malloc(Valwidth+1);
1027 if ( ThisElement == NULL ) IOHBTerminate("Insufficient memory for ThisElement.");
1028 *(ThisElement+Valwidth) = (char) NULL;
1029 count = 0;
1030 for (i=0;i<Valcrd;i++)
1031 {
1032 fgets(line, BUFSIZ, in_file);
1033 if ( sscanf(line,"%*s") < 0 )
1034 IOHBTerminate("iohb.c: Null (or blank) line in value data region of HB file.\n");
1035 if (Valflag == 'D') {
1036 while( strchr(line,'D') ) *strchr(line,'D') = 'E';
1037 }
1038 col = 0;
1039 for (ind = 0;ind<Valperline;ind++)
1040 {
1041 if (count == Nentries) break;
1042 ThisElement = &val[count*Valwidth];
1043 strncpy(ThisElement,line+col,Valwidth);
1044 /*strncpy(ThisElement,substr(line,col,Valwidth),Valwidth);*/
1045 if ( Valflag != 'F' && strchr(ThisElement,'E') == NULL ) {
1046 /* insert a char prefix for exp */
1047 last = strlen(ThisElement);
1048 for (j=last+1;j>=0;j--) {
1049 ThisElement[j] = ThisElement[j-1];
1050 if ( ThisElement[j] == '+' || ThisElement[j] == '-' ) {
1051 ThisElement[j-1] = Valflag;
1052 break;
1053 }
1054 }
1055 }
1056 count++; col += Valwidth;
1057 }
1058 }
1059 }
1060
1061 return 1;
1062 }
1063
readHB_newmat_char(const char * filename,int * M,int * N,int * nonzeros,int ** colptr,int ** rowind,char ** val,char ** Valfmt)1064 int readHB_newmat_char(const char* filename, int* M, int* N, int* nonzeros, int** colptr,
1065 int** rowind, char** val, char** Valfmt)
1066 {
1067 FILE *in_file;
1068 int Nrhs;
1069 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
1070 int Valperline, Valwidth, Valprec;
1071 int Valflag; /* Indicates 'E','D', or 'F' float format */
1072 char Title[73], Key[9], Type[4], Rhstype[4];
1073 char Ptrfmt[17], Indfmt[17], Rhsfmt[21];
1074
1075 if ((in_file = fopen( filename, "r")) == NULL) {
1076 fprintf(stderr,"Error: Cannot open file: %s\n",filename);
1077 return 0;
1078 }
1079
1080 *Valfmt = (char *)malloc(21*sizeof(char));
1081 if ( *Valfmt == NULL ) IOHBTerminate("Insufficient memory for Valfmt.");
1082 readHB_header(in_file, Title, Key, Type, M, N, nonzeros, &Nrhs,
1083 Ptrfmt, Indfmt, (*Valfmt), Rhsfmt,
1084 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
1085 fclose(in_file);
1086 ParseRfmt(*Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
1087
1088 *colptr = (int *)malloc((*N+1)*sizeof(int));
1089 if ( *colptr == NULL ) IOHBTerminate("Insufficient memory for colptr.\n");
1090 *rowind = (int *)malloc(*nonzeros*sizeof(int));
1091 if ( *rowind == NULL ) IOHBTerminate("Insufficient memory for rowind.\n");
1092 if ( Type[0] == 'C' ) {
1093 /*
1094 fprintf(stderr, "Warning: Reading complex data from HB file %s.\n",filename);
1095 fprintf(stderr, " Real and imaginary parts will be interlaced in val[].\n");
1096 */
1097 /* Malloc enough space for real AND imaginary parts of val[] */
1098 *val = (char *)malloc(*nonzeros*Valwidth*sizeof(char)*2);
1099 if ( *val == NULL ) IOHBTerminate("Insufficient memory for val.\n");
1100 } else {
1101 if ( Type[0] != 'P' ) {
1102 /* Malloc enough space for real array val[] */
1103 *val = (char *)malloc(*nonzeros*Valwidth*sizeof(char));
1104 if ( *val == NULL ) IOHBTerminate("Insufficient memory for val.\n");
1105 }
1106 } /* No val[] space needed if pattern only */
1107 return readHB_mat_char(filename, *colptr, *rowind, *val, *Valfmt);
1108
1109 }
1110
readHB_aux_char(const char * filename,const char AuxType,char b[])1111 int readHB_aux_char(const char* filename, const char AuxType, char b[])
1112 {
1113 /****************************************************************************/
1114 /* This function opens and reads the specified file, placing auxilary */
1115 /* vector(s) of the given type (if available) in b : */
1116 /* Return value is the number of vectors successfully read. */
1117 /* */
1118 /* AuxType = 'F' full right-hand-side vector(s) */
1119 /* AuxType = 'G' initial Guess vector(s) */
1120 /* AuxType = 'X' eXact solution vector(s) */
1121 /* */
1122 /* ---------- */
1123 /* **CAVEAT** */
1124 /* ---------- */
1125 /* Parsing real formats from Fortran is tricky, and this file reader */
1126 /* does not claim to be foolproof. It has been tested for cases when */
1127 /* the real values are printed consistently and evenly spaced on each */
1128 /* line, with Fixed (F), and Exponential (E or D) formats. */
1129 /* */
1130 /* ** If the input file does not adhere to the H/B format, the ** */
1131 /* ** results will be unpredictable. ** */
1132 /* */
1133 /****************************************************************************/
1134 FILE *in_file;
1135 int i,j,n,maxcol,start,stride,col,last,linel,nvecs,rhsi;
1136 int Nrow, Ncol, Nnzero, Nentries,Nrhs;
1137 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
1138 int Rhsperline, Rhswidth, Rhsprec;
1139 int Rhsflag;
1140 char Title[73], Key[9], Type[4], Rhstype[4];
1141 char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
1142 char line[BUFSIZ];
1143 char *ThisElement;
1144
1145 if ((in_file = fopen( filename, "r")) == NULL) {
1146 fprintf(stderr,"Error: Cannot open file: %s\n",filename);
1147 return 0;
1148 }
1149
1150 readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
1151 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
1152 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
1153
1154 if (Nrhs <= 0)
1155 {
1156 fprintf(stderr, "Warn: Attempt to read auxillary vector(s) when none are present.\n");
1157 return 0;
1158 }
1159 if (Rhstype[0] != 'F' )
1160 {
1161 fprintf(stderr,"Warn: Attempt to read auxillary vector(s) which are not stored in Full form.\n");
1162 fprintf(stderr," Rhs must be specified as full. \n");
1163 return 0;
1164 }
1165
1166 /* If reading complex data, allow for interleaved real and imaginary values. */
1167 if ( Type[0] == 'C' ) {
1168 Nentries = 2*Nrow;
1169 } else {
1170 Nentries = Nrow;
1171 }
1172
1173 nvecs = 1;
1174
1175 if ( Rhstype[1] == 'G' ) nvecs++;
1176 if ( Rhstype[2] == 'X' ) nvecs++;
1177
1178 if ( AuxType == 'G' && Rhstype[1] != 'G' ) {
1179 fprintf(stderr, "Warn: Attempt to read auxillary Guess vector(s) when none are present.\n");
1180 return 0;
1181 }
1182 if ( AuxType == 'X' && Rhstype[2] != 'X' ) {
1183 fprintf(stderr, "Warn: Attempt to read auxillary eXact solution vector(s) when none are present.\n");
1184 return 0;
1185 }
1186
1187 ParseRfmt(Rhsfmt, &Rhsperline, &Rhswidth, &Rhsprec,&Rhsflag);
1188 maxcol = Rhsperline*Rhswidth;
1189
1190 /* Lines to skip before starting to read RHS values... */
1191 n = Ptrcrd + Indcrd + Valcrd;
1192
1193 for (i = 0; i < n; i++)
1194 fgets(line, BUFSIZ, in_file);
1195
1196 /* start - number of initial aux vector entries to skip */
1197 /* to reach first vector requested */
1198 /* stride - number of aux vector entries to skip between */
1199 /* requested vectors */
1200 if ( AuxType == 'F' ) start = 0;
1201 else if ( AuxType == 'G' ) start = Nentries;
1202 else start = (nvecs-1)*Nentries;
1203 stride = (nvecs-1)*Nentries;
1204
1205 fgets(line, BUFSIZ, in_file);
1206 linel= strchr(line,'\n')-line;
1207 if ( sscanf(line,"%*s") < 0 )
1208 IOHBTerminate("iohb.c: Null (or blank) line in auxillary vector data region of HB file.\n");
1209 col = 0;
1210 /* Skip to initial offset */
1211
1212 for (i=0;i<start;i++) {
1213 col += Rhswidth;
1214 if ( col >= ( maxcol<linel?maxcol:linel ) ) {
1215 fgets(line, BUFSIZ, in_file);
1216 linel= strchr(line,'\n')-line;
1217 if ( sscanf(line,"%*s") < 0 )
1218 IOHBTerminate("iohb.c: Null (or blank) line in auxillary vector data region of HB file.\n");
1219 col = 0;
1220 }
1221 }
1222
1223 if (Rhsflag == 'D') {
1224 while( strchr(line,'D') ) *strchr(line,'D') = 'E';
1225 }
1226 /* Read a vector of desired type, then skip to next */
1227 /* repeating to fill Nrhs vectors */
1228
1229 for (rhsi=0;rhsi<Nrhs;rhsi++) {
1230
1231 for (i=0;i<Nentries;i++) {
1232 if ( col >= ( maxcol<linel?maxcol:linel ) ) {
1233 fgets(line, BUFSIZ, in_file);
1234 linel= strchr(line,'\n')-line;
1235 if ( sscanf(line,"%*s") < 0 )
1236 IOHBTerminate("iohb.c: Null (or blank) line in auxillary vector data region of HB file.\n");
1237 if (Rhsflag == 'D') {
1238 while( strchr(line,'D') ) *strchr(line,'D') = 'E';
1239 }
1240 col = 0;
1241 }
1242 ThisElement = &b[i*Rhswidth];
1243 strncpy(ThisElement,line+col,Rhswidth);
1244 if ( Rhsflag != 'F' && strchr(ThisElement,'E') == NULL ) {
1245 /* insert a char prefix for exp */
1246 last = strlen(ThisElement);
1247 for (j=last+1;j>=0;j--) {
1248 ThisElement[j] = ThisElement[j-1];
1249 if ( ThisElement[j] == '+' || ThisElement[j] == '-' ) {
1250 ThisElement[j-1] = Rhsflag;
1251 break;
1252 }
1253 }
1254 }
1255 col += Rhswidth;
1256 }
1257 b+=Nentries*Rhswidth;
1258
1259 /* Skip any interleaved Guess/eXact vectors */
1260
1261 for (i=0;i<stride;i++) {
1262 col += Rhswidth;
1263 if ( col >= ( maxcol<linel?maxcol:linel ) ) {
1264 fgets(line, BUFSIZ, in_file);
1265 linel= strchr(line,'\n')-line;
1266 if ( sscanf(line,"%*s") < 0 )
1267 IOHBTerminate("iohb.c: Null (or blank) line in auxillary vector data region of HB file.\n");
1268 col = 0;
1269 }
1270 }
1271
1272 }
1273
1274
1275 fclose(in_file);
1276 return Nrhs;
1277 }
1278
readHB_newaux_char(const char * filename,const char AuxType,char ** b,char ** Rhsfmt)1279 int readHB_newaux_char(const char* filename, const char AuxType, char** b, char** Rhsfmt)
1280 {
1281 FILE *in_file;
1282 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
1283 int Nrow,Ncol,Nnzero,Nrhs;
1284 int Rhsperline, Rhswidth, Rhsprec;
1285 int Rhsflag;
1286 char Title[73], Key[9], Type[4], Rhstype[4];
1287 char Ptrfmt[17], Indfmt[17], Valfmt[21];
1288
1289 if ((in_file = fopen( filename, "r")) == NULL) {
1290 fprintf(stderr,"Error: Cannot open file: %s\n",filename);
1291 return 0;
1292 }
1293
1294 *Rhsfmt = (char *)malloc(21*sizeof(char));
1295 if ( *Rhsfmt == NULL ) IOHBTerminate("Insufficient memory for Rhsfmt.");
1296 readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
1297 Ptrfmt, Indfmt, Valfmt, (*Rhsfmt),
1298 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
1299 fclose(in_file);
1300 if ( Nrhs == 0 ) {
1301 fprintf(stderr,"Warn: Requested read of aux vector(s) when none are present.\n");
1302 return 0;
1303 } else {
1304 ParseRfmt(*Rhsfmt,&Rhsperline,&Rhswidth,&Rhsprec,&Rhsflag);
1305 if ( Type[0] == 'C' ) {
1306 fprintf(stderr, "Warning: Reading complex aux vector(s) from HB file %s.",filename);
1307 fprintf(stderr, " Real and imaginary parts will be interlaced in b[].");
1308 *b = (char *)malloc(Nrow*Nrhs*Rhswidth*sizeof(char)*2);
1309 if ( *b == NULL ) IOHBTerminate("Insufficient memory for rhs.\n");
1310 return readHB_aux_char(filename, AuxType, *b);
1311 } else {
1312 *b = (char *)malloc(Nrow*Nrhs*Rhswidth*sizeof(char));
1313 if ( *b == NULL ) IOHBTerminate("Insufficient memory for rhs.\n");
1314 return readHB_aux_char(filename, AuxType, *b);
1315 }
1316 }
1317 }
1318
writeHB_mat_char(const char * filename,int M,int N,int nz,const int colptr[],const int rowind[],const char val[],int Nrhs,const char rhs[],const char guess[],const char exact[],const char * Title,const char * Key,const char * Type,char * Ptrfmt,char * Indfmt,char * Valfmt,char * Rhsfmt,const char * Rhstype)1319 int writeHB_mat_char(const char* filename, int M, int N,
1320 int nz, const int colptr[], const int rowind[],
1321 const char val[], int Nrhs, const char rhs[],
1322 const char guess[], const char exact[],
1323 const char* Title, const char* Key, const char* Type,
1324 char* Ptrfmt, char* Indfmt, char* Valfmt, char* Rhsfmt,
1325 const char* Rhstype)
1326 {
1327 /****************************************************************************/
1328 /* The writeHB function opens the named file and writes the specified */
1329 /* matrix and optional right-hand-side(s) to that file in Harwell-Boeing */
1330 /* format. */
1331 /* */
1332 /* For a description of the Harwell Boeing standard, see: */
1333 /* Duff, et al., ACM TOMS Vol.15, No.1, March 1989 */
1334 /* */
1335 /****************************************************************************/
1336 FILE *out_file;
1337 int i,j,acount,linemod,entry,offset;
1338 int totcrd, ptrcrd, indcrd, valcrd, rhscrd;
1339 int nvalentries, nrhsentries;
1340 int Ptrperline, Ptrwidth, Indperline, Indwidth;
1341 int Rhsperline, Rhswidth, Rhsprec;
1342 int Rhsflag;
1343 int Valperline, Valwidth, Valprec;
1344 int Valflag; /* Indicates 'E','D', or 'F' float format */
1345 char pformat[16],iformat[16],vformat[19],rformat[19];
1346
1347 if ( Type[0] == 'C' ) {
1348 nvalentries = 2*nz;
1349 nrhsentries = 2*M;
1350 } else {
1351 nvalentries = nz;
1352 nrhsentries = M;
1353 }
1354
1355 if ( filename != NULL ) {
1356 if ( (out_file = fopen( filename, "w")) == NULL ) {
1357 fprintf(stderr,"Error: Cannot open file: %s\n",filename);
1358 return 0;
1359 }
1360 } else out_file = stdout;
1361
1362 if ( Ptrfmt == NULL ) Ptrfmt = "(8I10)";
1363 ParseIfmt(Ptrfmt,&Ptrperline,&Ptrwidth);
1364 sprintf(pformat,"%%%dd",Ptrwidth);
1365
1366 if ( Indfmt == NULL ) Indfmt = Ptrfmt;
1367 ParseIfmt(Indfmt,&Indperline,&Indwidth);
1368 sprintf(iformat,"%%%dd",Indwidth);
1369
1370 if ( Type[0] != 'P' ) { /* Skip if pattern only */
1371 if ( Valfmt == NULL ) Valfmt = "(4E20.13)";
1372 ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
1373 sprintf(vformat,"%%%ds",Valwidth);
1374 }
1375
1376 ptrcrd = (N+1)/Ptrperline;
1377 if ( (N+1)%Ptrperline != 0) ptrcrd++;
1378
1379 indcrd = nz/Indperline;
1380 if ( nz%Indperline != 0) indcrd++;
1381
1382 valcrd = nvalentries/Valperline;
1383 if ( nvalentries%Valperline != 0) valcrd++;
1384
1385 if ( Nrhs > 0 ) {
1386 if ( Rhsfmt == NULL ) Rhsfmt = Valfmt;
1387 ParseRfmt(Rhsfmt,&Rhsperline,&Rhswidth,&Rhsprec, &Rhsflag);
1388 sprintf(rformat,"%%%ds",Rhswidth);
1389 rhscrd = nrhsentries/Rhsperline;
1390 if ( nrhsentries%Rhsperline != 0) rhscrd++;
1391 if ( Rhstype[1] == 'G' ) rhscrd+=rhscrd;
1392 if ( Rhstype[2] == 'X' ) rhscrd+=rhscrd;
1393 rhscrd*=Nrhs;
1394 } else rhscrd = 0;
1395
1396 totcrd = 4+ptrcrd+indcrd+valcrd+rhscrd;
1397
1398
1399 /* Print header information: */
1400
1401 fprintf(out_file,"%-72s%-8s\n%14d%14d%14d%14d%14d\n",Title, Key, totcrd,
1402 ptrcrd, indcrd, valcrd, rhscrd);
1403 fprintf(out_file,"%3s%11s%14d%14d%14d\n",Type," ", M, N, nz);
1404 fprintf(out_file,"%-16s%-16s%-20s", Ptrfmt, Indfmt, Valfmt);
1405 if ( Nrhs != 0 ) {
1406 /* Print Rhsfmt on fourth line and */
1407 /* optional fifth header line for auxillary vector information: */
1408 fprintf(out_file,"%-20s\n%-14s%d\n",Rhsfmt,Rhstype,Nrhs);
1409 } else fprintf(out_file,"\n");
1410
1411 offset = 1-_SP_base; /* if base 0 storage is declared (via macro definition), */
1412 /* then storage entries are offset by 1 */
1413
1414 /* Print column pointers: */
1415 for (i=0;i<N+1;i++)
1416 {
1417 entry = colptr[i]+offset;
1418 fprintf(out_file,pformat,entry);
1419 if ( (i+1)%Ptrperline == 0 ) fprintf(out_file,"\n");
1420 }
1421
1422 if ( (N+1) % Ptrperline != 0 ) fprintf(out_file,"\n");
1423
1424 /* Print row indices: */
1425 for (i=0;i<nz;i++)
1426 {
1427 entry = rowind[i]+offset;
1428 fprintf(out_file,iformat,entry);
1429 if ( (i+1)%Indperline == 0 ) fprintf(out_file,"\n");
1430 }
1431
1432 if ( nz % Indperline != 0 ) fprintf(out_file,"\n");
1433
1434 /* Print values: */
1435
1436 if ( Type[0] != 'P' ) { /* Skip if pattern only */
1437 for (i=0;i<nvalentries;i++)
1438 {
1439 fprintf(out_file,vformat,val+i*Valwidth);
1440 if ( (i+1)%Valperline == 0 ) fprintf(out_file,"\n");
1441 }
1442
1443 if ( nvalentries % Valperline != 0 ) fprintf(out_file,"\n");
1444
1445 /* Print right hand sides: */
1446 acount = 1;
1447 linemod=0;
1448 if ( Nrhs > 0 ) {
1449 for (j=0;j<Nrhs;j++) {
1450 for (i=0;i<nrhsentries;i++)
1451 {
1452 fprintf(out_file,rformat,rhs+i*Rhswidth);
1453 if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
1454 }
1455 if ( acount%Rhsperline != linemod ) {
1456 fprintf(out_file,"\n");
1457 linemod = (acount-1)%Rhsperline;
1458 }
1459 if ( Rhstype[1] == 'G' ) {
1460 for (i=0;i<nrhsentries;i++)
1461 {
1462 fprintf(out_file,rformat,guess+i*Rhswidth);
1463 if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
1464 }
1465 if ( acount%Rhsperline != linemod ) {
1466 fprintf(out_file,"\n");
1467 linemod = (acount-1)%Rhsperline;
1468 }
1469 }
1470 if ( Rhstype[2] == 'X' ) {
1471 for (i=0;i<nrhsentries;i++)
1472 {
1473 fprintf(out_file,rformat,exact+i*Rhswidth);
1474 if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
1475 }
1476 if ( acount%Rhsperline != linemod ) {
1477 fprintf(out_file,"\n");
1478 linemod = (acount-1)%Rhsperline;
1479 }
1480 }
1481 }
1482 }
1483
1484 }
1485
1486 if ( fclose(out_file) != 0){
1487 fprintf(stderr,"Error closing file in writeHB_mat_char().\n");
1488 return 0;
1489 } else return 1;
1490
1491 }
1492
ParseIfmt(char * fmt,int * perline,int * width)1493 int ParseIfmt(char* fmt, int* perline, int* width)
1494 {
1495 /*************************************************/
1496 /* Parse an *integer* format field to determine */
1497 /* width and number of elements per line. */
1498 /*************************************************/
1499 char *tmp;
1500 if (fmt == NULL ) {
1501 *perline = 0; *width = 0; return 0;
1502 }
1503 upcase(fmt);
1504 tmp = strchr(fmt,'(');
1505 tmp = substr(fmt,tmp - fmt + 1, strchr(fmt,'I') - tmp - 1);
1506 *perline = atoi(tmp);
1507 tmp = strchr(fmt,'I');
1508 tmp = substr(fmt,tmp - fmt + 1, strchr(fmt,')') - tmp - 1);
1509 return *width = atoi(tmp);
1510 }
1511
ParseRfmt(char * fmt,int * perline,int * width,int * prec,int * flag)1512 int ParseRfmt(char* fmt, int* perline, int* width, int* prec, int* flag)
1513 {
1514 /*************************************************/
1515 /* Parse a *real* format field to determine */
1516 /* width and number of elements per line. */
1517 /* Also sets flag indicating 'E' 'F' 'P' or 'D' */
1518 /* format. */
1519 /*************************************************/
1520 char* tmp;
1521 char* tmp2;
1522 char* tmp3;
1523 int len;
1524
1525 if (fmt == NULL ) {
1526 *perline = 0;
1527 *width = 0;
1528 flag = NULL;
1529 return 0;
1530 }
1531
1532 upcase(fmt);
1533 if (strchr(fmt,'(') != NULL) fmt = strchr(fmt,'(');
1534 if (strchr(fmt,')') != NULL) {
1535 tmp2 = strchr(fmt,')');
1536 while ( strchr(tmp2+1,')') != NULL ) {
1537 tmp2 = strchr(tmp2+1,')');
1538 }
1539 *(tmp2+1) = (int) NULL;
1540 }
1541 if (strchr(fmt,'P') != NULL) /* Remove any scaling factor, which */
1542 { /* affects output only, not input */
1543 if (strchr(fmt,'(') != NULL) {
1544 tmp = strchr(fmt,'P');
1545 if ( *(++tmp) == ',' ) tmp++;
1546 tmp3 = strchr(fmt,'(')+1;
1547 len = tmp-tmp3;
1548 tmp2 = tmp3;
1549 while ( *(tmp2+len) != (int) NULL ) {
1550 *tmp2=*(tmp2+len);
1551 tmp2++;
1552 }
1553 *(strchr(fmt,')')+1) = (int) NULL;
1554 }
1555 }
1556 if (strchr(fmt,'E') != NULL) {
1557 *flag = 'E';
1558 } else if (strchr(fmt,'D') != NULL) {
1559 *flag = 'D';
1560 } else if (strchr(fmt,'F') != NULL) {
1561 *flag = 'F';
1562 } else {
1563 fprintf(stderr,"Real format %s in H/B file not supported.\n",fmt);
1564 return 0;
1565 }
1566 tmp = strchr(fmt,'(');
1567 tmp = substr(fmt,tmp - fmt + 1, strchr(fmt,*flag) - tmp - 1);
1568 *perline = atoi(tmp);
1569 tmp = strchr(fmt,*flag);
1570 if ( strchr(fmt,'.') ) {
1571 *prec = atoi( substr( fmt, strchr(fmt,'.') - fmt + 1, strchr(fmt,')') - strchr(fmt,'.')-1) );
1572 tmp = substr(fmt,tmp - fmt + 1, strchr(fmt,'.') - tmp - 1);
1573 } else {
1574 tmp = substr(fmt,tmp - fmt + 1, strchr(fmt,')') - tmp - 1);
1575 }
1576 return *width = atoi(tmp);
1577 }
1578
substr(const char * S,const int pos,const int len)1579 char* substr(const char* S, const int pos, const int len)
1580 {
1581 int i;
1582 char *SubS;
1583 if ( pos+len <= strlen(S)) {
1584 SubS = (char *)malloc(len+1);
1585 if ( SubS == NULL ) IOHBTerminate("Insufficient memory for SubS.");
1586 for (i=0;i<len;i++) SubS[i] = S[pos+i];
1587 SubS[len] = (char) NULL;
1588 } else {
1589 SubS = NULL;
1590 }
1591 return SubS;
1592 }
1593
1594 #include<ctype.h>
upcase(char * S)1595 void upcase(char* S)
1596 {
1597 /* Convert S to uppercase */
1598 int i,len;
1599 if ( S == NULL ) return;
1600 len = strlen(S);
1601 for (i=0;i< len;i++)
1602 S[i] = toupper(S[i]);
1603 }
1604
IOHBTerminate(char * message)1605 void IOHBTerminate(char* message)
1606 {
1607 fprintf(stderr,message);
1608 exit(1);
1609 }
1610
1611