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