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