Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
iohb.c
Go to the documentation of this file.
1 /*
2 Fri 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 
224 char* substr(const char* S, const int pos, const int len);
225 void upcase(char* S);
226 void IOHBTerminate(char* message);
227 
228 int 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 
293 int 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 
366 int 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 
512 int 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 
543 int 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 
708 int 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 
732 int 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 
914 int 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 
1060 int 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 
1107 int 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 
1275 int 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 
1315 int 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 
1489 int 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 
1508 int 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 
1575 char* 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>
1591 void 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 
1600 void IOHBTerminate(char* message)
1601 {
1602  fprintf(stderr,message);
1603  exit(1);
1604 }
1605 
substr
char * substr(const char *S, const int pos, const int len)
Definition: iohb.c:1575
ParseRfmt
int ParseRfmt(char *fmt, int *perline, int *width, int *prec, int *flag)
Definition: iohb.c:1508
writeHB_mat_char
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
IOHBTerminate
void IOHBTerminate(char *message)
Definition: iohb.c:1600
readHB_header
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
readHB_mat_char
int readHB_mat_char(const char *filename, int colptr[], int rowind[], char val[], char *Valfmt)
Definition: iohb.c:914
_SP_base
#define _SP_base
Definition: iohb.c:213
readHB_mat_double
int readHB_mat_double(const char *filename, int colptr[], int rowind[], double val[])
Definition: iohb.c:366
readHB_aux_char
int readHB_aux_char(const char *filename, const char AuxType, char b[])
Definition: iohb.c:1107
readHB_newmat_double
int readHB_newmat_double(const char *filename, int *M, int *N, int *nonzeros, int **colptr, int **rowind, double **val)
Definition: iohb.c:512
readHB_newaux_char
int readHB_newaux_char(const char *filename, const char AuxType, char **b, char **Rhsfmt)
Definition: iohb.c:1275
N
const int N
Definition: LL/hypre_UnitTest.cpp:73
readHB_newmat_char
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
readHB_aux_double
int readHB_aux_double(const char *filename, const char AuxType, double b[])
Definition: iohb.c:543
upcase
void upcase(char *S)
Definition: iohb.c:1591
readHB_info
int readHB_info(const char *filename, int *M, int *N, int *nz, char **Type, int *Nrhs)
Definition: iohb.c:228
writeHB_mat_double
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
ParseIfmt
int ParseIfmt(char *fmt, int *perline, int *width)
Definition: iohb.c:1489
n
int n
readHB_newaux_double
int readHB_newaux_double(const char *filename, const char AuxType, double **b)
Definition: iohb.c:708
Generated on Thu Nov 21 2024 09:19:19 for Ifpack Package Browser (Single Doxygen Collection) by   doxygen 1.8.5