Planet
navi homePPSaboutscreenshotsdownloaddevelopmentforum

source: downloads/boost_1_34_1/libs/graph/example/iohb.c @ 29

Last change on this file since 29 was 29, checked in by landauf, 16 years ago

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