1
2 #include <stdio.h>
3 #include <string.h>
4 #include "mmio.h"
5 #include "hbio.h"
6 #include "lusolio.h"
7
8 /* Utility routines to read matrix files in the Coordinate Text File format*/
9
ctf_read_A(char * filename,int maxm,int maxn,int maxnz,int * m,int * n,int * nnzero,int * iA,int * jA,REAL * Aij)10 MYBOOL ctf_read_A(char *filename, int maxm, int maxn, int maxnz,
11 int *m, int *n, int *nnzero, int *iA, int *jA, REAL *Aij)
12 {
13 FILE *iofile;
14 int eof;
15 char buffer[100];
16 int k, i, j;
17 REAL Ak;
18 MYBOOL filldata;
19
20 *nnzero = 0;
21 *m = 0;
22 *n = 0;
23
24 iofile = fopen(filename, "r" );
25 if(iofile == NULL) {
26 printf("A file %s does not exist\n", filename);
27 return( FALSE );
28 }
29
30 filldata = (MYBOOL) !((iA == NULL) && (jA == NULL) && (Aij == NULL));
31 eof = TRUE;
32 for (k = 1; k <= maxnz; k++) {
33 eof = feof(iofile);
34 if(eof)
35 break;
36 eof = fscanf(iofile, "%d %d %s", &i, &j, buffer);
37 if(eof == 0 || eof == EOF || i <= 0 || j <= 0 || strlen(buffer) == 0)
38 break;
39 Ak = atof(buffer);
40 (*nnzero)++;
41 if (filldata) {
42 iA[k] = i;
43 jA[k] = j;
44 Aij[k] = Ak;
45 }
46 if (i > *m) *m = i;
47 if (j > *n) *n = j;
48 }
49 fclose( iofile );
50 if(!eof) {
51 printf("Too much data in A file. Increase maxnz\n");
52 printf("Current maxnz = %d\n", maxnz);
53 return( FALSE );
54 }
55 printf("A read successfully\n");
56 printf("m = %d n = %d nnzero = %d\n",
57 *m, *n, *nnzero);
58 if (*m > maxm || *n > maxn) {
59 printf("However, matrix dimensions exceed maxm or maxn\n");
60 return( FALSE );
61 }
62 return( TRUE );
63 }
64
ctf_size_A(char * filename,int * m,int * n,int * nnzero)65 MYBOOL ctf_size_A(char *filename, int *m, int *n, int *nnzero)
66 {
67 int maxint = 200000000;
68
69 return( ctf_read_A(filename, maxint, maxint, maxint,
70 m, n, nnzero, NULL, NULL, NULL) );
71 }
72
ctf_read_b(char * filename,int m,REAL * b)73 MYBOOL ctf_read_b(char *filename, int m, REAL *b)
74 {
75 FILE *iofile;
76 int eof;
77 char buffer[100];
78 int i;
79
80 iofile = fopen(filename, "r");
81 if(iofile == NULL) {
82 printf("b file %s does not exist\n", filename);
83 return( FALSE );
84 }
85
86 for (i = 1; i <= m; i++) {
87 if(feof(iofile))
88 goto x350;
89 eof = fscanf(iofile, "%s", buffer);
90 if(eof == 0 || eof == EOF)
91 goto x350;
92 b[i] = atof(buffer);
93 }
94
95 fclose( iofile );
96 printf("b read successfully\n");
97 return( TRUE );
98
99 x350:
100 fclose( iofile );
101 printf("Not enough data in b file.\n");
102 return( FALSE );
103 }
104
105
106 /* Utility routines to read matrix files in the MatrixMarket format*/
107 #define mmf_recsize 255
mmf_read_A(char * filename,int maxM,int maxN,int maxnz,int * M,int * N,int * nz,int * iA,int * jA,REAL * Aij)108 MYBOOL mmf_read_A(char *filename, int maxM, int maxN, int maxnz,
109 int *M, int *N, int *nz, int *iA, int *jA, REAL *Aij)
110 {
111 MM_typecode matcode;
112 FILE *f;
113 int i, k, ret_code, ival, jval;
114 REAL Aval;
115 MYBOOL status = FALSE, ispat, filldata;
116 char buf[mmf_recsize];
117
118 f = fopen(filename, "r");
119 if(f == NULL)
120 return( status );
121
122 if(mm_read_banner(f, &matcode) != 0) {
123 printf("Could not process Matrix Market banner.\n");
124 goto x900;
125 }
126
127 /* Screen matrix types since LUSOL only supports a
128 subset of the Matrix Market data types. */
129 if(mm_is_complex(matcode) || mm_is_pattern(matcode)) {
130 printf("Sorry, this application does not support ");
131 printf("Market Market type: [%s]\n", mm_typecode_to_str(matcode));
132 goto x900;
133 }
134
135 /* Verify that we have sufficient array storage */
136 filldata = (MYBOOL) !((iA == NULL) && (jA == NULL) && (Aij == NULL));
137 if(filldata && maxN > 1 && jA == NULL) {
138 printf("Market Market insufficient array storage specified\n");
139 goto x900;
140 }
141
142 /* Find out size of sparse matrix .... */
143 ret_code = mm_read_mtx_crd_size(f, M, N, nz);
144 if(ret_code != 0 || !filldata || (*M > maxM) || (*N > maxN) || (*nz > maxnz)) {
145 status = !filldata;
146 goto x900;
147 }
148
149
150 /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
151 /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
152 /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
153
154 /* Read dense matrix in column order */
155 ispat = (MYBOOL) mm_is_pattern(matcode);
156 k = 1;
157 if(mm_is_dense(matcode)) {
158 maxN = MIN(maxN, *N);
159 for (jval = 1; jval <= maxN; jval++) {
160 for (i = 1; i <= *M; i++) {
161 if(fgets(buf, mmf_recsize-1, f) == NULL)
162 break;
163 if(sscanf(buf, "%lg\n", &Aval) == 0)
164 break;
165 if(Aval != 0) {
166 if(iA != NULL)
167 iA[k] = i;
168 if(jA != NULL)
169 jA[k] = jval;
170
171 /* Make sure we handle dense vector reading properly */
172 if(iA == NULL && jA == NULL)
173 Aij[i] = Aval;
174 else
175 Aij[k] = Aval;
176 k++;
177 }
178 }
179 }
180 }
181 /* Read sparse matrix by coordinate */
182 else {
183 for (i = 1; i <= *nz; i++) {
184 if(fgets(buf, mmf_recsize-1, f) == NULL)
185 break;
186 if(buf[0] == '%')
187 continue;
188 if(ispat) {
189 if(sscanf(buf, "%d %d\n", &ival, &jval) == 0)
190 continue;
191 Aij[k] = 1.0;
192 }
193 else
194 if(sscanf(buf, "%d %d %lg\n", &ival, &jval, &Aval) == 0)
195 continue;
196
197 /* Check if it is a nonzero and we are within column dimension */
198 if(Aval != 0 && jval <= maxN) {
199 Aij[k] = Aval;
200 if(iA != NULL)
201 iA[k] = ival;
202 if(jA != NULL)
203 jA[k] = jval;
204 k++;
205 }
206 }
207 }
208 *nz = k - 1;
209
210 /* Handle case where only the lower triangular parts are given */
211 if(!mm_is_general(matcode)) {
212 if((M != N) || (maxN != maxM) || (2*(*nz) > maxnz)) {
213 printf("Market Market cannot fill in symmetry data\n");
214 goto x900;
215 }
216 ispat = mm_is_skew(matcode);
217 for(i = 1; i <= *nz; i++) {
218 iA[k] = jA[i];
219 jA[k] = iA[i];
220 if(ispat)
221 Aij[k] = -Aij[i];
222 else
223 Aij[k] = Aij[i];
224 k++;
225 }
226 *nz = k - 1;
227 }
228 status = TRUE;
229
230 /* Finish up */
231 x900:
232 fclose( f );
233 return( status );
234 }
235
mmf_size_A(char * filename,int * M,int * N,int * nz)236 MYBOOL mmf_size_A(char *filename, int *M, int *N, int *nz)
237 {
238 int maxint = 200000000;
239
240 return( mmf_read_A(filename, maxint, maxint, maxint,
241 M, N, nz, NULL, NULL, NULL) );
242 }
243
mmf_read_b(char * filename,int m,REAL * b)244 MYBOOL mmf_read_b(char *filename, int m, REAL *b)
245 {
246 int im, jn, nnzero;
247
248 return( mmf_read_A(filename, m, 1, m,
249 &im, &jn, &nnzero, NULL, NULL, b));
250 }
251
252
253 /* Utility routines to read matrix files in Harwell-Boeing format*/
254
hbf_read_A(char * filename,int maxM,int maxN,int maxnz,int * M,int * N,int * nz,int * iA,int * jA,REAL * Aij)255 MYBOOL hbf_read_A(char *filename, int maxM, int maxN, int maxnz,
256 int *M, int *N, int *nz, int *iA, int *jA, REAL *Aij)
257 {
258 MYBOOL success;
259
260 success = hbf_size_A(filename, M, N, nz);
261 if(!success)
262 return( success );
263
264 Aij[1] = 0;
265 success = (MYBOOL) readHB_mat_double(filename, jA, iA-1, Aij-1);
266
267 /* Test if we have a pattern matrix and fill it with all zeros */
268 if(Aij[1] == 0) {
269 int i;
270 for(i = 1; i <= *nz; i++)
271 Aij[i] = 1;
272 }
273
274 /* Expand the column nz counts to triplet format */
275 if(success) {
276 int i, j, ii, k;
277 k = *nz;
278 for(j = *N; j > 0; j--) {
279 ii = jA[j];
280 for(i = jA[j-1]; i < ii; i++, k--)
281 jA[k] = j;
282 }
283 }
284 return( success );
285 }
286
hbf_size_A(char * filename,int * M,int * N,int * nz)287 MYBOOL hbf_size_A(char *filename, int *M, int *N, int *nz)
288 {
289 int Nrhs;
290 char *Type;
291
292 return( (MYBOOL) readHB_info(filename, M, N, nz, &Type, &Nrhs) );
293 }
294
hbf_read_b(char * filename,int m,REAL * b)295 MYBOOL hbf_read_b(char *filename, int m, REAL *b)
296 {
297 return( (MYBOOL) readHB_aux_double(filename, 'F', b) ); /* Same format as matrix */
298 }
299