1 /*! \file
2 Copyright (c) 2003, The Regents of the University of California, through
3 Lawrence Berkeley National Laboratory (subject to receipt of any required
4 approvals from U.S. Dept. of Energy)
5
6 All rights reserved.
7
8 The source code is distributed under BSD license, see the file License.txt
9 at the top-level directory.
10 */
11
12 /*! @file creadhb.c
13 * \brief Read a matrix stored in Harwell-Boeing format
14 *
15 * <pre>
16 * -- SuperLU routine (version 2.0) --
17 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
18 * and Lawrence Berkeley National Lab.
19 * November 15, 1997
20 *
21 * Purpose
22 * =======
23 *
24 * Read a COMPLEX PRECISION matrix stored in Harwell-Boeing format
25 * as described below.
26 *
27 * Line 1 (A72,A8)
28 * Col. 1 - 72 Title (TITLE)
29 * Col. 73 - 80 Key (KEY)
30 *
31 * Line 2 (5I14)
32 * Col. 1 - 14 Total number of lines excluding header (TOTCRD)
33 * Col. 15 - 28 Number of lines for pointers (PTRCRD)
34 * Col. 29 - 42 Number of lines for row (or variable) indices (INDCRD)
35 * Col. 43 - 56 Number of lines for numerical values (VALCRD)
36 * Col. 57 - 70 Number of lines for right-hand sides (RHSCRD)
37 * (including starting guesses and solution vectors
38 * if present)
39 * (zero indicates no right-hand side data is present)
40 *
41 * Line 3 (A3, 11X, 4I14)
42 * Col. 1 - 3 Matrix type (see below) (MXTYPE)
43 * Col. 15 - 28 Number of rows (or variables) (NROW)
44 * Col. 29 - 42 Number of columns (or elements) (NCOL)
45 * Col. 43 - 56 Number of row (or variable) indices (NNZERO)
46 * (equal to number of entries for assembled matrices)
47 * Col. 57 - 70 Number of elemental matrix entries (NELTVL)
48 * (zero in the case of assembled matrices)
49 * Line 4 (2A16, 2A20)
50 * Col. 1 - 16 Format for pointers (PTRFMT)
51 * Col. 17 - 32 Format for row (or variable) indices (INDFMT)
52 * Col. 33 - 52 Format for numerical values of coefficient matrix (VALFMT)
53 * Col. 53 - 72 Format for numerical values of right-hand sides (RHSFMT)
54 *
55 * Line 5 (A3, 11X, 2I14) Only present if there are right-hand sides present
56 * Col. 1 Right-hand side type:
57 * F for full storage or M for same format as matrix
58 * Col. 2 G if a starting vector(s) (Guess) is supplied. (RHSTYP)
59 * Col. 3 X if an exact solution vector(s) is supplied.
60 * Col. 15 - 28 Number of right-hand sides (NRHS)
61 * Col. 29 - 42 Number of row indices (NRHSIX)
62 * (ignored in case of unassembled matrices)
63 *
64 * The three character type field on line 3 describes the matrix type.
65 * The following table lists the permitted values for each of the three
66 * characters. As an example of the type field, RSA denotes that the matrix
67 * is real, symmetric, and assembled.
68 *
69 * First Character:
70 * R Real matrix
71 * C Complex matrix
72 * P Pattern only (no numerical values supplied)
73 *
74 * Second Character:
75 * S Symmetric
76 * U Unsymmetric
77 * H Hermitian
78 * Z Skew symmetric
79 * R Rectangular
80 *
81 * Third Character:
82 * A Assembled
83 * E Elemental matrices (unassembled)
84 *
85 * </pre>
86 */
87 #include <stdio.h>
88 #include <stdlib.h>
89 #include "slu_cdefs.h"
90
91
92 /*! \brief Eat up the rest of the current line */
cDumpLine(FILE * fp)93 int cDumpLine(FILE *fp)
94 {
95 register int c;
96 while ((c = fgetc(fp)) != '\n') ;
97 return 0;
98 }
99
cParseIntFormat(char * buf,int * num,int * size)100 int cParseIntFormat(char *buf, int *num, int *size)
101 {
102 char *tmp;
103
104 tmp = buf;
105 while (*tmp++ != '(') ;
106 sscanf(tmp, "%d", num);
107 while (*tmp != 'I' && *tmp != 'i') ++tmp;
108 ++tmp;
109 sscanf(tmp, "%d", size);
110 return 0;
111 }
112
cParseFloatFormat(char * buf,int * num,int * size)113 int cParseFloatFormat(char *buf, int *num, int *size)
114 {
115 char *tmp, *period;
116
117 tmp = buf;
118 while (*tmp++ != '(') ;
119 *num = atoi(tmp); /*sscanf(tmp, "%d", num);*/
120 while (*tmp != 'E' && *tmp != 'e' && *tmp != 'D' && *tmp != 'd'
121 && *tmp != 'F' && *tmp != 'f') {
122 /* May find kP before nE/nD/nF, like (1P6F13.6). In this case the
123 num picked up refers to P, which should be skipped. */
124 if (*tmp=='p' || *tmp=='P') {
125 ++tmp;
126 *num = atoi(tmp); /*sscanf(tmp, "%d", num);*/
127 } else {
128 ++tmp;
129 }
130 }
131 ++tmp;
132 period = tmp;
133 while (*period != '.' && *period != ')') ++period ;
134 *period = '\0';
135 *size = atoi(tmp); /*sscanf(tmp, "%2d", size);*/
136
137 return 0;
138 }
139
ReadVector(FILE * fp,int n,int * where,int perline,int persize)140 static int ReadVector(FILE *fp, int n, int *where, int perline, int persize)
141 {
142 register int i, j, item;
143 char tmp, buf[100];
144
145 i = 0;
146 while (i < n) {
147 fgets(buf, 100, fp); /* read a line at a time */
148 for (j=0; j<perline && i<n; j++) {
149 tmp = buf[(j+1)*persize]; /* save the char at that place */
150 buf[(j+1)*persize] = 0; /* null terminate */
151 item = atoi(&buf[j*persize]);
152 buf[(j+1)*persize] = tmp; /* recover the char at that place */
153 where[i++] = item - 1;
154 }
155 }
156
157 return 0;
158 }
159
160 /*! \brief Read complex numbers as pairs of (real, imaginary) */
cReadValues(FILE * fp,int n,complex * destination,int perline,int persize)161 int cReadValues(FILE *fp, int n, complex *destination, int perline, int persize)
162 {
163 register int i, j, k, s, pair;
164 register float realpart;
165 char tmp, buf[100];
166
167 i = pair = 0;
168 while (i < n) {
169 fgets(buf, 100, fp); /* read a line at a time */
170 for (j=0; j<perline && i<n; j++) {
171 tmp = buf[(j+1)*persize]; /* save the char at that place */
172 buf[(j+1)*persize] = 0; /* null terminate */
173 s = j*persize;
174 for (k = 0; k < persize; ++k) /* No D_ format in C */
175 if ( buf[s+k] == 'D' || buf[s+k] == 'd' ) buf[s+k] = 'E';
176 if ( pair == 0 ) {
177 /* The value is real part */
178 realpart = atof(&buf[s]);
179 pair = 1;
180 } else {
181 /* The value is imaginary part */
182 destination[i].r = realpart;
183 destination[i++].i = atof(&buf[s]);
184 pair = 0;
185 }
186 buf[(j+1)*persize] = tmp; /* recover the char at that place */
187 }
188 }
189
190 return 0;
191 }
192
193 /*! \brief
194 *
195 * <pre>
196 * On input, nonz/nzval/rowind/colptr represents lower part of a symmetric
197 * matrix. On exit, it represents the full matrix with lower and upper parts.
198 * </pre>
199 */
200 static void
FormFullA(int n,int * nonz,complex ** nzval,int ** rowind,int ** colptr)201 FormFullA(int n, int *nonz, complex **nzval, int **rowind, int **colptr)
202 {
203 register int i, j, k, col, new_nnz;
204 int *t_rowind, *t_colptr, *al_rowind, *al_colptr, *a_rowind, *a_colptr;
205 int *marker;
206 complex *t_val, *al_val, *a_val;
207
208 al_rowind = *rowind;
209 al_colptr = *colptr;
210 al_val = *nzval;
211
212 if ( !(marker =(int *) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
213 ABORT("SUPERLU_MALLOC fails for marker[]");
214 if ( !(t_colptr = (int *) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
215 ABORT("SUPERLU_MALLOC t_colptr[]");
216 if ( !(t_rowind = (int *) SUPERLU_MALLOC( *nonz * sizeof(int)) ) )
217 ABORT("SUPERLU_MALLOC fails for t_rowind[]");
218 if ( !(t_val = (complex*) SUPERLU_MALLOC( *nonz * sizeof(complex)) ) )
219 ABORT("SUPERLU_MALLOC fails for t_val[]");
220
221 /* Get counts of each column of T, and set up column pointers */
222 for (i = 0; i < n; ++i) marker[i] = 0;
223 for (j = 0; j < n; ++j) {
224 for (i = al_colptr[j]; i < al_colptr[j+1]; ++i)
225 ++marker[al_rowind[i]];
226 }
227 t_colptr[0] = 0;
228 for (i = 0; i < n; ++i) {
229 t_colptr[i+1] = t_colptr[i] + marker[i];
230 marker[i] = t_colptr[i];
231 }
232
233 /* Transpose matrix A to T */
234 for (j = 0; j < n; ++j)
235 for (i = al_colptr[j]; i < al_colptr[j+1]; ++i) {
236 col = al_rowind[i];
237 t_rowind[marker[col]] = j;
238 t_val[marker[col]] = al_val[i];
239 ++marker[col];
240 }
241
242 new_nnz = *nonz * 2 - n;
243 if ( !(a_colptr = (int *) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
244 ABORT("SUPERLU_MALLOC a_colptr[]");
245 if ( !(a_rowind = (int *) SUPERLU_MALLOC( new_nnz * sizeof(int)) ) )
246 ABORT("SUPERLU_MALLOC fails for a_rowind[]");
247 if ( !(a_val = (complex*) SUPERLU_MALLOC( new_nnz * sizeof(complex)) ) )
248 ABORT("SUPERLU_MALLOC fails for a_val[]");
249
250 a_colptr[0] = 0;
251 k = 0;
252 for (j = 0; j < n; ++j) {
253 for (i = t_colptr[j]; i < t_colptr[j+1]; ++i) {
254 if ( t_rowind[i] != j ) { /* not diagonal */
255 a_rowind[k] = t_rowind[i];
256 a_val[k] = t_val[i];
257 #ifdef DEBUG
258 if ( fabs(a_val[k]) < 4.047e-300 )
259 printf("%5d: %e\n", k, a_val[k]);
260 #endif
261 ++k;
262 }
263 }
264
265 for (i = al_colptr[j]; i < al_colptr[j+1]; ++i) {
266 a_rowind[k] = al_rowind[i];
267 a_val[k] = al_val[i];
268 #ifdef DEBUG
269 if ( fabs(a_val[k]) < 4.047e-300 )
270 printf("%5d: %e\n", k, a_val[k]);
271 #endif
272 ++k;
273 }
274
275 a_colptr[j+1] = k;
276 }
277
278 printf("FormFullA: new_nnz = %d, k = %d\n", new_nnz, k);
279
280 SUPERLU_FREE(al_val);
281 SUPERLU_FREE(al_rowind);
282 SUPERLU_FREE(al_colptr);
283 SUPERLU_FREE(marker);
284 SUPERLU_FREE(t_val);
285 SUPERLU_FREE(t_rowind);
286 SUPERLU_FREE(t_colptr);
287
288 *nzval = a_val;
289 *rowind = a_rowind;
290 *colptr = a_colptr;
291 *nonz = new_nnz;
292 }
293
294 void
creadhb(FILE * fp,int * nrow,int * ncol,int * nonz,complex ** nzval,int ** rowind,int ** colptr)295 creadhb(FILE *fp, int *nrow, int *ncol, int *nonz,
296 complex **nzval, int **rowind, int **colptr)
297 {
298
299 register int i, numer_lines = 0, rhscrd = 0;
300 int tmp, colnum, colsize, rownum, rowsize, valnum, valsize;
301 char buf[100], type[4], key[10];
302 int sym;
303
304 /* Line 1 */
305 fgets(buf, 100, fp);
306 fputs(buf, stdout);
307 #if 0
308 fscanf(fp, "%72c", buf); buf[72] = 0;
309 printf("Title: %s", buf);
310 fscanf(fp, "%8c", key); key[8] = 0;
311 printf("Key: %s\n", key);
312 cDumpLine(fp);
313 #endif
314
315 /* Line 2 */
316 for (i=0; i<5; i++) {
317 fscanf(fp, "%14c", buf); buf[14] = 0;
318 sscanf(buf, "%d", &tmp);
319 if (i == 3) numer_lines = tmp;
320 if (i == 4 && tmp) rhscrd = tmp;
321 }
322 cDumpLine(fp);
323
324 /* Line 3 */
325 fscanf(fp, "%3c", type);
326 fscanf(fp, "%11c", buf); /* pad */
327 type[3] = 0;
328 #ifdef DEBUG
329 printf("Matrix type %s\n", type);
330 #endif
331
332 fscanf(fp, "%14c", buf); sscanf(buf, "%d", nrow);
333 fscanf(fp, "%14c", buf); sscanf(buf, "%d", ncol);
334 fscanf(fp, "%14c", buf); sscanf(buf, "%d", nonz);
335 fscanf(fp, "%14c", buf); sscanf(buf, "%d", &tmp);
336
337 if (tmp != 0)
338 printf("This is not an assembled matrix!\n");
339 if (*nrow != *ncol)
340 printf("Matrix is not square.\n");
341 cDumpLine(fp);
342
343 /* Allocate storage for the three arrays ( nzval, rowind, colptr ) */
344 callocateA(*ncol, *nonz, nzval, rowind, colptr);
345
346 /* Line 4: format statement */
347 fscanf(fp, "%16c", buf);
348 cParseIntFormat(buf, &colnum, &colsize);
349 fscanf(fp, "%16c", buf);
350 cParseIntFormat(buf, &rownum, &rowsize);
351 fscanf(fp, "%20c", buf);
352 cParseFloatFormat(buf, &valnum, &valsize);
353 fscanf(fp, "%20c", buf);
354 cDumpLine(fp);
355
356 /* Line 5: right-hand side */
357 if ( rhscrd ) cDumpLine(fp); /* skip RHSFMT */
358
359 #ifdef DEBUG
360 printf("%d rows, %d nonzeros\n", *nrow, *nonz);
361 printf("colnum %d, colsize %d\n", colnum, colsize);
362 printf("rownum %d, rowsize %d\n", rownum, rowsize);
363 printf("valnum %d, valsize %d\n", valnum, valsize);
364 #endif
365
366 ReadVector(fp, *ncol+1, *colptr, colnum, colsize);
367 ReadVector(fp, *nonz, *rowind, rownum, rowsize);
368 if ( numer_lines ) {
369 cReadValues(fp, *nonz, *nzval, valnum, valsize);
370 }
371
372 sym = (type[1] == 'S' || type[1] == 's');
373 if ( sym ) {
374 FormFullA(*ncol, nonz, nzval, rowind, colptr);
375 }
376
377 fclose(fp);
378 }
379
380