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