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 dreadrb.c
13 * \brief Read a matrix stored in Rutherford-Boeing format
14 *
15 * <pre>
16 * -- SuperLU routine (version 4.0) --
17 * Lawrence Berkeley National Laboratory.
18 * June 30, 2009
19 * </pre>
20 *
21 * Purpose
22 * =======
23 *
24 * Read a DOUBLE PRECISION matrix stored in Rutherford-Boeing format
25 * as described below.
26 *
27 * Line 1 (A72, A8)
28 * Col. 1 - 72 Title (TITLE)
29 * Col. 73 - 80 Matrix name / identifier (MTRXID)
30 *
31 * Line 2 (I14, 3(1X, I13))
32 * Col. 1 - 14 Total number of lines excluding header (TOTCRD)
33 * Col. 16 - 28 Number of lines for pointers (PTRCRD)
34 * Col. 30 - 42 Number of lines for row (or variable) indices (INDCRD)
35 * Col. 44 - 56 Number of lines for numerical values (VALCRD)
36 *
37 * Line 3 (A3, 11X, 4(1X, I13))
38 * Col. 1 - 3 Matrix type (see below) (MXTYPE)
39 * Col. 15 - 28 Compressed Column: Number of rows (NROW)
40 * Elemental: Largest integer used to index variable (MVAR)
41 * Col. 30 - 42 Compressed Column: Number of columns (NCOL)
42 * Elemental: Number of element matrices (NELT)
43 * Col. 44 - 56 Compressed Column: Number of entries (NNZERO)
44 * Elemental: Number of variable indeces (NVARIX)
45 * Col. 58 - 70 Compressed Column: Unused, explicitly zero
46 * Elemental: Number of elemental matrix entries (NELTVL)
47 *
48 * Line 4 (2A16, A20)
49 * Col. 1 - 16 Fortran format for pointers (PTRFMT)
50 * Col. 17 - 32 Fortran format for row (or variable) indices (INDFMT)
51 * Col. 33 - 52 Fortran format for numerical values of coefficient matrix
52 * (VALFMT)
53 * (blank in the case of matrix patterns)
54 *
55 * The three character type field on line 3 describes the matrix type.
56 * The following table lists the permitted values for each of the three
57 * characters. As an example of the type field, RSA denotes that the matrix
58 * is real, symmetric, and assembled.
59 *
60 * First Character:
61 * R Real matrix
62 * C Complex matrix
63 * I integer matrix
64 * P Pattern only (no numerical values supplied)
65 * Q Pattern only (numerical values supplied in associated auxiliary value
66 * file)
67 *
68 * Second Character:
69 * S Symmetric
70 * U Unsymmetric
71 * H Hermitian
72 * Z Skew symmetric
73 * R Rectangular
74 *
75 * Third Character:
76 * A Compressed column form
77 * E Elemental form
78 *
79 * </pre>
80 */
81
82 #include <stdio.h>
83 #include <stdlib.h>
84 #include "slu_ddefs.h"
85
86
87 /*! \brief Eat up the rest of the current line */
dDumpLine(FILE * fp)88 static int dDumpLine(FILE *fp)
89 {
90 register int c;
91 while ((c = fgetc(fp)) != '\n') ;
92 return 0;
93 }
94
dParseIntFormat(char * buf,int * num,int * size)95 static int dParseIntFormat(char *buf, int *num, int *size)
96 {
97 char *tmp;
98
99 tmp = buf;
100 while (*tmp++ != '(') ;
101 sscanf(tmp, "%d", num);
102 while (*tmp != 'I' && *tmp != 'i') ++tmp;
103 ++tmp;
104 sscanf(tmp, "%d", size);
105 return 0;
106 }
107
dParseFloatFormat(char * buf,int * num,int * size)108 static int dParseFloatFormat(char *buf, int *num, int *size)
109 {
110 char *tmp, *period;
111
112 tmp = buf;
113 while (*tmp++ != '(') ;
114 *num = atoi(tmp); /*sscanf(tmp, "%d", num);*/
115 while (*tmp != 'E' && *tmp != 'e' && *tmp != 'D' && *tmp != 'd'
116 && *tmp != 'F' && *tmp != 'f') {
117 /* May find kP before nE/nD/nF, like (1P6F13.6). In this case the
118 num picked up refers to P, which should be skipped. */
119 if (*tmp=='p' || *tmp=='P') {
120 ++tmp;
121 *num = atoi(tmp); /*sscanf(tmp, "%d", num);*/
122 } else {
123 ++tmp;
124 }
125 }
126 ++tmp;
127 period = tmp;
128 while (*period != '.' && *period != ')') ++period ;
129 *period = '\0';
130 *size = atoi(tmp); /*sscanf(tmp, "%2d", size);*/
131
132 return 0;
133 }
134
ReadVector(FILE * fp,int n,int * where,int perline,int persize)135 static int ReadVector(FILE *fp, int n, int *where, int perline, int persize)
136 {
137 register int i, j, item;
138 char tmp, buf[100];
139
140 i = 0;
141 while (i < n) {
142 fgets(buf, 100, fp); /* read a line at a time */
143 for (j=0; j<perline && i<n; j++) {
144 tmp = buf[(j+1)*persize]; /* save the char at that place */
145 buf[(j+1)*persize] = 0; /* null terminate */
146 item = atoi(&buf[j*persize]);
147 buf[(j+1)*persize] = tmp; /* recover the char at that place */
148 where[i++] = item - 1;
149 }
150 }
151
152 return 0;
153 }
154
dReadValues(FILE * fp,int n,double * destination,int perline,int persize)155 static int dReadValues(FILE *fp, int n, double *destination, int perline,
156 int persize)
157 {
158 register int i, j, k, s;
159 char tmp, buf[100];
160
161 i = 0;
162 while (i < n) {
163 fgets(buf, 100, fp); /* read a line at a time */
164 for (j=0; j<perline && i<n; j++) {
165 tmp = buf[(j+1)*persize]; /* save the char at that place */
166 buf[(j+1)*persize] = 0; /* null terminate */
167 s = j*persize;
168 for (k = 0; k < persize; ++k) /* No D_ format in C */
169 if ( buf[s+k] == 'D' || buf[s+k] == 'd' ) buf[s+k] = 'E';
170 destination[i++] = atof(&buf[s]);
171 buf[(j+1)*persize] = tmp; /* recover the char at that place */
172 }
173 }
174
175 return 0;
176 }
177
178
179
180 /*! \brief
181 *
182 * <pre>
183 * On input, nonz/nzval/rowind/colptr represents lower part of a symmetric
184 * matrix. On exit, it represents the full matrix with lower and upper parts.
185 * </pre>
186 */
187 static void
FormFullA(int n,int * nonz,double ** nzval,int ** rowind,int ** colptr)188 FormFullA(int n, int *nonz, double **nzval, int **rowind, int **colptr)
189 {
190 register int i, j, k, col, new_nnz;
191 int *t_rowind, *t_colptr, *al_rowind, *al_colptr, *a_rowind, *a_colptr;
192 int *marker;
193 double *t_val, *al_val, *a_val;
194
195 al_rowind = *rowind;
196 al_colptr = *colptr;
197 al_val = *nzval;
198
199 if ( !(marker =(int *) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
200 ABORT("SUPERLU_MALLOC fails for marker[]");
201 if ( !(t_colptr = (int *) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
202 ABORT("SUPERLU_MALLOC t_colptr[]");
203 if ( !(t_rowind = (int *) SUPERLU_MALLOC( *nonz * sizeof(int)) ) )
204 ABORT("SUPERLU_MALLOC fails for t_rowind[]");
205 if ( !(t_val = (double*) SUPERLU_MALLOC( *nonz * sizeof(double)) ) )
206 ABORT("SUPERLU_MALLOC fails for t_val[]");
207
208 /* Get counts of each column of T, and set up column pointers */
209 for (i = 0; i < n; ++i) marker[i] = 0;
210 for (j = 0; j < n; ++j) {
211 for (i = al_colptr[j]; i < al_colptr[j+1]; ++i)
212 ++marker[al_rowind[i]];
213 }
214 t_colptr[0] = 0;
215 for (i = 0; i < n; ++i) {
216 t_colptr[i+1] = t_colptr[i] + marker[i];
217 marker[i] = t_colptr[i];
218 }
219
220 /* Transpose matrix A to T */
221 for (j = 0; j < n; ++j)
222 for (i = al_colptr[j]; i < al_colptr[j+1]; ++i) {
223 col = al_rowind[i];
224 t_rowind[marker[col]] = j;
225 t_val[marker[col]] = al_val[i];
226 ++marker[col];
227 }
228
229 new_nnz = *nonz * 2 - n;
230 if ( !(a_colptr = (int *) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
231 ABORT("SUPERLU_MALLOC a_colptr[]");
232 if ( !(a_rowind = (int *) SUPERLU_MALLOC( new_nnz * sizeof(int)) ) )
233 ABORT("SUPERLU_MALLOC fails for a_rowind[]");
234 if ( !(a_val = (double*) SUPERLU_MALLOC( new_nnz * sizeof(double)) ) )
235 ABORT("SUPERLU_MALLOC fails for a_val[]");
236
237 a_colptr[0] = 0;
238 k = 0;
239 for (j = 0; j < n; ++j) {
240 for (i = t_colptr[j]; i < t_colptr[j+1]; ++i) {
241 if ( t_rowind[i] != j ) { /* not diagonal */
242 a_rowind[k] = t_rowind[i];
243 a_val[k] = t_val[i];
244 #ifdef DEBUG
245 if ( fabs(a_val[k]) < 4.047e-300 )
246 printf("%5d: %e\n", k, a_val[k]);
247 #endif
248 ++k;
249 }
250 }
251
252 for (i = al_colptr[j]; i < al_colptr[j+1]; ++i) {
253 a_rowind[k] = al_rowind[i];
254 a_val[k] = al_val[i];
255 #ifdef DEBUG
256 if ( fabs(a_val[k]) < 4.047e-300 )
257 printf("%5d: %e\n", k, a_val[k]);
258 #endif
259 ++k;
260 }
261
262 a_colptr[j+1] = k;
263 }
264
265 printf("FormFullA: new_nnz = %d, k = %d\n", new_nnz, k);
266
267 SUPERLU_FREE(al_val);
268 SUPERLU_FREE(al_rowind);
269 SUPERLU_FREE(al_colptr);
270 SUPERLU_FREE(marker);
271 SUPERLU_FREE(t_val);
272 SUPERLU_FREE(t_rowind);
273 SUPERLU_FREE(t_colptr);
274
275 *nzval = a_val;
276 *rowind = a_rowind;
277 *colptr = a_colptr;
278 *nonz = new_nnz;
279 }
280
281 void
dreadrb(int * nrow,int * ncol,int * nonz,double ** nzval,int ** rowind,int ** colptr)282 dreadrb(int *nrow, int *ncol, int *nonz,
283 double **nzval, int **rowind, int **colptr)
284 {
285
286 register int i, numer_lines = 0;
287 int tmp, colnum, colsize, rownum, rowsize, valnum, valsize;
288 char buf[100], type[4];
289 int sym;
290 FILE *fp;
291
292 fp = stdin;
293
294 /* Line 1 */
295 fgets(buf, 100, fp);
296 fputs(buf, stdout);
297
298 /* Line 2 */
299 for (i=0; i<4; i++) {
300 fscanf(fp, "%14c", buf); buf[14] = 0;
301 sscanf(buf, "%d", &tmp);
302 if (i == 3) numer_lines = tmp;
303 }
304 dDumpLine(fp);
305
306 /* Line 3 */
307 fscanf(fp, "%3c", type);
308 fscanf(fp, "%11c", buf); /* pad */
309 type[3] = 0;
310 #ifdef DEBUG
311 printf("Matrix type %s\n", type);
312 #endif
313
314 fscanf(fp, "%14c", buf); sscanf(buf, "%d", nrow);
315 fscanf(fp, "%14c", buf); sscanf(buf, "%d", ncol);
316 fscanf(fp, "%14c", buf); sscanf(buf, "%d", nonz);
317 fscanf(fp, "%14c", buf); sscanf(buf, "%d", &tmp);
318
319 if (tmp != 0)
320 printf("This is not an assembled matrix!\n");
321 if (*nrow != *ncol)
322 printf("Matrix is not square.\n");
323 dDumpLine(fp);
324
325 /* Allocate storage for the three arrays ( nzval, rowind, colptr ) */
326 dallocateA(*ncol, *nonz, nzval, rowind, colptr);
327
328 /* Line 4: format statement */
329 fscanf(fp, "%16c", buf);
330 dParseIntFormat(buf, &colnum, &colsize);
331 fscanf(fp, "%16c", buf);
332 dParseIntFormat(buf, &rownum, &rowsize);
333 fscanf(fp, "%20c", buf);
334 dParseFloatFormat(buf, &valnum, &valsize);
335 dDumpLine(fp);
336
337 #ifdef DEBUG
338 printf("%d rows, %d nonzeros\n", *nrow, *nonz);
339 printf("colnum %d, colsize %d\n", colnum, colsize);
340 printf("rownum %d, rowsize %d\n", rownum, rowsize);
341 printf("valnum %d, valsize %d\n", valnum, valsize);
342 #endif
343
344 ReadVector(fp, *ncol+1, *colptr, colnum, colsize);
345 ReadVector(fp, *nonz, *rowind, rownum, rowsize);
346 if ( numer_lines ) {
347 dReadValues(fp, *nonz, *nzval, valnum, valsize);
348 }
349
350 sym = (type[1] == 'S' || type[1] == 's');
351 if ( sym ) {
352 FormFullA(*ncol, nonz, nzval, rowind, colptr);
353 }
354
355 fclose(fp);
356 }
357