4
* -- SuperLU routine (version 2.0) --
5
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
6
* and Lawrence Berkeley National Lab.
15
/* Eat up the rest of the current line */
16
int cDumpLine(FILE *fp)
19
while ((c = fgetc(fp)) != '\n') ;
23
int cParseIntFormat(char *buf, int *num, int *size)
28
while (*tmp++ != '(') ;
29
sscanf(tmp, "%d", num);
30
while (*tmp != 'I' && *tmp != 'i') ++tmp;
32
sscanf(tmp, "%d", size);
36
int cParseFloatFormat(char *buf, int *num, int *size)
41
while (*tmp++ != '(') ;
42
*num = atoi(tmp); /*sscanf(tmp, "%d", num);*/
43
while (*tmp != 'E' && *tmp != 'e' && *tmp != 'D' && *tmp != 'd'
44
&& *tmp != 'F' && *tmp != 'f') {
45
/* May find kP before nE/nD/nF, like (1P6F13.6). In this case the
46
num picked up refers to P, which should be skipped. */
47
if (*tmp=='p' || *tmp=='P') {
49
*num = atoi(tmp); /*sscanf(tmp, "%d", num);*/
56
while (*period != '.' && *period != ')') ++period ;
58
*size = atoi(tmp); /*sscanf(tmp, "%2d", size);*/
63
int cReadVector(FILE *fp, int n, int *where, int perline, int persize)
65
register int i, j, item;
70
fgets(buf, 100, fp); /* read a line at a time */
71
for (j=0; j<perline && i<n; j++) {
72
tmp = buf[(j+1)*persize]; /* save the char at that place */
73
buf[(j+1)*persize] = 0; /* null terminate */
74
item = atoi(&buf[j*persize]);
75
buf[(j+1)*persize] = tmp; /* recover the char at that place */
76
where[i++] = item - 1;
83
/* Read complex numbers as pairs of (real, imaginary) */
84
int cReadValues(FILE *fp, int n, complex *destination, int perline, int persize)
86
register int i, j, k, s, pair;
87
register float realpart;
92
fgets(buf, 100, fp); /* read a line at a time */
93
for (j=0; j<perline && i<n; j++) {
94
tmp = buf[(j+1)*persize]; /* save the char at that place */
95
buf[(j+1)*persize] = 0; /* null terminate */
97
for (k = 0; k < persize; ++k) /* No D_ format in C */
98
if ( buf[s+k] == 'D' || buf[s+k] == 'd' ) buf[s+k] = 'E';
100
/* The value is real part */
101
realpart = atof(&buf[s]);
104
/* The value is imaginary part */
105
destination[i].r = realpart;
106
destination[i++].i = atof(&buf[s]);
109
buf[(j+1)*persize] = tmp; /* recover the char at that place */
118
creadhb(int *nrow, int *ncol, int *nonz,
119
complex **nzval, int **rowind, int **colptr)
125
* Read a COMPLEX PRECISION matrix stored in Harwell-Boeing format
126
* as described below.
129
* Col. 1 - 72 Title (TITLE)
130
* Col. 73 - 80 Key (KEY)
133
* Col. 1 - 14 Total number of lines excluding header (TOTCRD)
134
* Col. 15 - 28 Number of lines for pointers (PTRCRD)
135
* Col. 29 - 42 Number of lines for row (or variable) indices (INDCRD)
136
* Col. 43 - 56 Number of lines for numerical values (VALCRD)
137
* Col. 57 - 70 Number of lines for right-hand sides (RHSCRD)
138
* (including starting guesses and solution vectors
140
* (zero indicates no right-hand side data is present)
142
* Line 3 (A3, 11X, 4I14)
143
* Col. 1 - 3 Matrix type (see below) (MXTYPE)
144
* Col. 15 - 28 Number of rows (or variables) (NROW)
145
* Col. 29 - 42 Number of columns (or elements) (NCOL)
146
* Col. 43 - 56 Number of row (or variable) indices (NNZERO)
147
* (equal to number of entries for assembled matrices)
148
* Col. 57 - 70 Number of elemental matrix entries (NELTVL)
149
* (zero in the case of assembled matrices)
150
* Line 4 (2A16, 2A20)
151
* Col. 1 - 16 Format for pointers (PTRFMT)
152
* Col. 17 - 32 Format for row (or variable) indices (INDFMT)
153
* Col. 33 - 52 Format for numerical values of coefficient matrix (VALFMT)
154
* Col. 53 - 72 Format for numerical values of right-hand sides (RHSFMT)
156
* Line 5 (A3, 11X, 2I14) Only present if there are right-hand sides present
157
* Col. 1 Right-hand side type:
158
* F for full storage or M for same format as matrix
159
* Col. 2 G if a starting vector(s) (Guess) is supplied. (RHSTYP)
160
* Col. 3 X if an exact solution vector(s) is supplied.
161
* Col. 15 - 28 Number of right-hand sides (NRHS)
162
* Col. 29 - 42 Number of row indices (NRHSIX)
163
* (ignored in case of unassembled matrices)
165
* The three character type field on line 3 describes the matrix type.
166
* The following table lists the permitted values for each of the three
167
* characters. As an example of the type field, RSA denotes that the matrix
168
* is real, symmetric, and assembled.
173
* P Pattern only (no numerical values supplied)
184
* E Elemental matrices (unassembled)
188
register int i, numer_lines = 0, rhscrd = 0;
189
int tmp, colnum, colsize, rownum, rowsize, valnum, valsize;
190
char buf[100], type[4], key[10];
199
fscanf(fp, "%72c", buf); buf[72] = 0;
200
printf("Title: %s", buf);
201
fscanf(fp, "%8c", key); key[8] = 0;
202
printf("Key: %s\n", key);
207
for (i=0; i<5; i++) {
208
fscanf(fp, "%14c", buf); buf[14] = 0;
209
sscanf(buf, "%d", &tmp);
210
if (i == 3) numer_lines = tmp;
211
if (i == 4 && tmp) rhscrd = tmp;
216
fscanf(fp, "%3c", type);
217
fscanf(fp, "%11c", buf); /* pad */
220
printf("Matrix type %s\n", type);
223
fscanf(fp, "%14c", buf); sscanf(buf, "%d", nrow);
224
fscanf(fp, "%14c", buf); sscanf(buf, "%d", ncol);
225
fscanf(fp, "%14c", buf); sscanf(buf, "%d", nonz);
226
fscanf(fp, "%14c", buf); sscanf(buf, "%d", &tmp);
229
printf("This is not an assembled matrix!\n");
231
printf("Matrix is not square.\n");
234
/* Allocate storage for the three arrays ( nzval, rowind, colptr ) */
235
callocateA(*ncol, *nonz, nzval, rowind, colptr);
237
/* Line 4: format statement */
238
fscanf(fp, "%16c", buf);
239
cParseIntFormat(buf, &colnum, &colsize);
240
fscanf(fp, "%16c", buf);
241
cParseIntFormat(buf, &rownum, &rowsize);
242
fscanf(fp, "%20c", buf);
243
cParseFloatFormat(buf, &valnum, &valsize);
244
fscanf(fp, "%20c", buf);
247
/* Line 5: right-hand side */
248
if ( rhscrd ) cDumpLine(fp); /* skip RHSFMT */
251
printf("%d rows, %d nonzeros\n", *nrow, *nonz);
252
printf("colnum %d, colsize %d\n", colnum, colsize);
253
printf("rownum %d, rowsize %d\n", rownum, rowsize);
254
printf("valnum %d, valsize %d\n", valnum, valsize);
257
cReadVector(fp, *ncol+1, *colptr, colnum, colsize);
258
cReadVector(fp, *nonz, *rowind, rownum, rowsize);
260
cReadValues(fp, *nonz, *nzval, valnum, valsize);