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