1 /*
2 Matrix Market I/O library for ANSI C
3
4 See https://math.nist.gov/MatrixMarket/ for details.
5 */
6
7 #include <stdlib.h>
8 #include <stdio.h>
9 #include <string.h>
10 #include <ctype.h>
11
12 #include "mmio.h"
13
14 static char mm_buffer[MM_MAX_LINE_LENGTH];
15
mm_read_unsymmetric_sparse(const char * fname,int * M_,int * N_,int * nz_,double ** val_,int ** I_,int ** J_)16 int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_, double **val_, int **I_, int **J_)
17 {
18 FILE *f;
19 MM_typecode matcode;
20 int M, N, nz;
21 int i;
22 double *val;
23 int *ia, *ja;
24
25 if ((f = fopen(fname, "r")) == NULL) return -1;
26
27 if (mm_read_banner(f, &matcode) != 0) {
28 printf("mm_read_unsymmetric: Could not process Matrix Market banner ");
29 printf(" in file [%s]\n", fname);
30 return -1;
31 }
32
33 if (!(mm_is_real(matcode) && mm_is_matrix(matcode) && mm_is_sparse(matcode))) {
34 fprintf(stderr, "This application does not support ");
35 fprintf(stderr, "Market Market type: [%s]\n", mm_typecode_to_str(matcode));
36 return -1;
37 }
38
39 /* find out size of sparse matrix: M, N, nz .... */
40
41 if (mm_read_mtx_crd_size(f, &M, &N, &nz) != 0) {
42 fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n");
43 return -1;
44 }
45
46 *M_ = M;
47 *N_ = N;
48 *nz_ = nz;
49
50 /* reserve memory for matrices */
51
52 ia = (int *)malloc(nz * sizeof(int));
53 ja = (int *)malloc(nz * sizeof(int));
54 val = (double *)malloc(nz * sizeof(double));
55
56 *val_ = val;
57 *I_ = ia;
58 *J_ = ja;
59
60 /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
61 /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
62 /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
63
64 for (i = 0; i < nz; i++) {
65 if (fscanf(f, "%d %d %lg\n", &ia[i], &ja[i], &val[i]) != 3) {
66 fprintf(stderr, "read_unsymmetric_sparse(): could not parse i, j and nonzero.\n");
67 return -1;
68 }
69 ia[i]--; /* adjust from 1-based to 0-based */
70 ja[i]--;
71 }
72 fclose(f);
73
74 return 0;
75 }
76
mm_is_valid(MM_typecode matcode)77 int mm_is_valid(MM_typecode matcode)
78 {
79 if (!mm_is_matrix(matcode)) return 0;
80 if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
81 if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
82 if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) || mm_is_skew(matcode))) return 0;
83 return 1;
84 }
85
mm_read_banner(FILE * f,MM_typecode * matcode)86 int mm_read_banner(FILE *f, MM_typecode *matcode)
87 {
88 char line[MM_MAX_LINE_LENGTH];
89 char banner[MM_MAX_TOKEN_LENGTH];
90 char mtx[MM_MAX_TOKEN_LENGTH];
91 char crd[MM_MAX_TOKEN_LENGTH];
92 char data_type[MM_MAX_TOKEN_LENGTH];
93 char storage_scheme[MM_MAX_TOKEN_LENGTH];
94 char *p;
95
96 mm_clear_typecode(matcode);
97
98 if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF;
99
100 if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type, storage_scheme) != 5) return MM_PREMATURE_EOF;
101
102 for (p = mtx; *p != '\0'; *p = (char)tolower(*p), p++); /* convert to lower case */
103 for (p = crd; *p != '\0'; *p = (char)tolower(*p), p++);
104 for (p = data_type; *p != '\0'; *p = (char)tolower(*p), p++);
105 for (p = storage_scheme; *p != '\0'; *p = (char)tolower(*p), p++);
106
107 /* check for banner */
108 if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0) return MM_NO_HEADER;
109
110 /* first field should be "mtx" */
111 if (strcmp(mtx, MM_MTX_STR) != 0) return MM_UNSUPPORTED_TYPE;
112 mm_set_matrix(matcode);
113
114 /* second field describes whether this is a sparse matrix (in coordinate
115 storgae) or a dense array */
116
117 if (strcmp(crd, MM_SPARSE_STR) == 0) mm_set_sparse(matcode);
118 else if (strcmp(crd, MM_DENSE_STR) == 0) mm_set_dense(matcode);
119 else return MM_UNSUPPORTED_TYPE;
120
121 /* third field */
122
123 if (strcmp(data_type, MM_REAL_STR) == 0) mm_set_real(matcode);
124 else if (strcmp(data_type, MM_COMPLEX_STR) == 0) mm_set_complex(matcode);
125 else if (strcmp(data_type, MM_PATTERN_STR) == 0) mm_set_pattern(matcode);
126 else if (strcmp(data_type, MM_INT_STR) == 0) mm_set_integer(matcode);
127 else return MM_UNSUPPORTED_TYPE;
128
129 /* fourth field */
130
131 if (strcmp(storage_scheme, MM_GENERAL_STR) == 0) mm_set_general(matcode);
132 else if (strcmp(storage_scheme, MM_SYMM_STR) == 0) mm_set_symmetric(matcode);
133 else if (strcmp(storage_scheme, MM_HERM_STR) == 0) mm_set_hermitian(matcode);
134 else if (strcmp(storage_scheme, MM_SKEW_STR) == 0) mm_set_skew(matcode);
135 else return MM_UNSUPPORTED_TYPE;
136
137 return 0;
138 }
139
mm_write_mtx_crd_size(FILE * f,int M,int N,int nz)140 int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
141 {
142 if (fprintf(f, "%d %d %d\n", M, N, nz) < 0) return MM_COULD_NOT_WRITE_FILE;
143 else return 0;
144 }
145
mm_read_mtx_crd_size(FILE * f,int * M,int * N,int * nz)146 int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz)
147 {
148 char line[MM_MAX_LINE_LENGTH];
149 int num_items_read;
150
151 /* set return null parameter values, in case we exit with errors */
152 *M = *N = *nz = 0;
153
154 /* now continue scanning until you reach the end-of-comments */
155 do {
156 if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF;
157 } while (line[0] == '%');
158
159 /* line[] is either blank or has M,N, nz */
160 if (sscanf(line, "%d %d %d", M, N, nz) == 3) return 0;
161
162 else do {
163 num_items_read = fscanf(f, "%d %d %d", M, N, nz);
164 if (num_items_read == EOF) return MM_PREMATURE_EOF;
165 } while (num_items_read != 3);
166
167 return 0;
168 }
169
mm_read_mtx_array_size(FILE * f,int * M,int * N)170 int mm_read_mtx_array_size(FILE *f, int *M, int *N)
171 {
172 char line[MM_MAX_LINE_LENGTH];
173 int num_items_read;
174 /* set return null parameter values, in case we exit with errors */
175 *M = *N = 0;
176
177 /* now continue scanning until you reach the end-of-comments */
178 do {
179 if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF;
180 } while (line[0] == '%');
181
182 /* line[] is either blank or has M,N, nz */
183 if (sscanf(line, "%d %d", M, N) == 2) return 0;
184
185 else /* we have a blank line */ do {
186 num_items_read = fscanf(f, "%d %d", M, N);
187 if (num_items_read == EOF) return MM_PREMATURE_EOF;
188 } while (num_items_read != 2);
189
190 return 0;
191 }
192
mm_write_mtx_array_size(FILE * f,int M,int N)193 int mm_write_mtx_array_size(FILE *f, int M, int N)
194 {
195 if (fprintf(f, "%d %d\n", M, N) < 0) return MM_COULD_NOT_WRITE_FILE;
196 else return 0;
197 }
198
199 /*-------------------------------------------------------------------------*/
200
201 /******************************************************************/
202 /* use when ia[], ja[], and val[] are already allocated */
203 /******************************************************************/
204
mm_read_mtx_crd_data(FILE * f,int M,int N,int nz,int ia[],int ja[],double val[],MM_typecode matcode)205 int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int ia[], int ja[], double val[], MM_typecode matcode)
206 {
207 int i;
208 if (mm_is_complex(matcode)) {
209 for (i = 0; i < nz; i++)
210 if (fscanf(f, "%d %d %lg %lg", &ia[i], &ja[i], &val[2 * i], &val[2 * i + 1]) != 4) return MM_PREMATURE_EOF;
211 } else if (mm_is_real(matcode)) {
212 for (i = 0; i < nz; i++) {
213 if (fscanf(f, "%d %d %lg\n", &ia[i], &ja[i], &val[i]) != 3) return MM_PREMATURE_EOF;
214 }
215 }
216
217 else if (mm_is_pattern(matcode)) {
218 for (i = 0; i < nz; i++)
219 if (fscanf(f, "%d %d", &ia[i], &ja[i]) != 2) return MM_PREMATURE_EOF;
220 } else return MM_UNSUPPORTED_TYPE;
221
222 return 0;
223 }
224
mm_read_mtx_crd_entry(FILE * f,int * ia,int * ja,double * real,double * imag,MM_typecode matcode)225 int mm_read_mtx_crd_entry(FILE *f, int *ia, int *ja, double *real, double *imag, MM_typecode matcode)
226 {
227 if (mm_is_complex(matcode)) {
228 if (fscanf(f, "%d %d %lg %lg", ia, ja, real, imag) != 4) return MM_PREMATURE_EOF;
229 } else if (mm_is_real(matcode)) {
230 if (fscanf(f, "%d %d %lg\n", ia, ja, real) != 3) return MM_PREMATURE_EOF;
231
232 }
233
234 else if (mm_is_pattern(matcode)) {
235 if (fscanf(f, "%d %d", ia, ja) != 2) return MM_PREMATURE_EOF;
236 } else return MM_UNSUPPORTED_TYPE;
237
238 return 0;
239 }
240
241 /************************************************************************
242 mm_read_mtx_crd() fills M, N, nz, array of values, and return
243 type code, e.g. 'MCRS'
244
245 if matrix is complex, values[] is of size 2*nz,
246 (nz pairs of real/imaginary values)
247 ************************************************************************/
248
mm_read_mtx_crd(char * fname,int * M,int * N,int * nz,int ** ia,int ** ja,double ** val,MM_typecode * matcode)249 int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **ia, int **ja, double **val, MM_typecode *matcode)
250 {
251 int ret_code;
252 FILE *f;
253
254 if (strcmp(fname, "stdin") == 0) f = stdin;
255 else if ((f = fopen(fname, "r")) == NULL) return MM_COULD_NOT_READ_FILE;
256
257 if ((ret_code = mm_read_banner(f, matcode)) != 0) return ret_code;
258
259 if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) && mm_is_matrix(*matcode))) return MM_UNSUPPORTED_TYPE;
260
261 if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0) return ret_code;
262
263 *ia = (int *)malloc(*nz * sizeof(int));
264 *ja = (int *)malloc(*nz * sizeof(int));
265 *val = NULL;
266
267 if (mm_is_complex(*matcode)) {
268 *val = (double *)malloc(*nz * 2 * sizeof(double));
269 ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *ia, *ja, *val, *matcode);
270 if (ret_code != 0) return ret_code;
271 } else if (mm_is_real(*matcode)) {
272 *val = (double *)malloc(*nz * sizeof(double));
273 ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *ia, *ja, *val, *matcode);
274 if (ret_code != 0) return ret_code;
275 }
276
277 else if (mm_is_pattern(*matcode)) {
278 ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *ia, *ja, *val, *matcode);
279 if (ret_code != 0) return ret_code;
280 }
281
282 if (f != stdin) fclose(f);
283 return 0;
284 }
285
mm_write_banner(FILE * f,MM_typecode matcode)286 int mm_write_banner(FILE *f, MM_typecode matcode)
287 {
288 char *str = mm_typecode_to_str(matcode);
289 int ret_code;
290
291 ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str);
292 if (ret_code < 0) return MM_COULD_NOT_WRITE_FILE;
293 else return 0;
294 }
295
mm_write_mtx_crd(char fname[],int M,int N,int nz,int ia[],int ja[],double val[],MM_typecode matcode)296 int mm_write_mtx_crd(char fname[], int M, int N, int nz, int ia[], int ja[], double val[], MM_typecode matcode)
297 {
298 FILE *f;
299 int i;
300
301 if (strcmp(fname, "stdout") == 0) f = stdout;
302 else if ((f = fopen(fname, "w")) == NULL) return MM_COULD_NOT_WRITE_FILE;
303
304 /* print banner followed by typecode */
305 fprintf(f, "%s ", MatrixMarketBanner);
306 fprintf(f, "%s\n", mm_typecode_to_str(matcode));
307
308 /* print matrix sizes and nonzeros */
309 fprintf(f, "%d %d %d\n", M, N, nz);
310
311 /* print values */
312 if (mm_is_pattern(matcode))
313 for (i = 0; i < nz; i++) fprintf(f, "%d %d\n", ia[i], ja[i]);
314 else if (mm_is_real(matcode))
315 for (i = 0; i < nz; i++) fprintf(f, "%d %d %20.16g\n", ia[i], ja[i], val[i]);
316 else if (mm_is_complex(matcode))
317 for (i = 0; i < nz; i++) fprintf(f, "%d %d %20.16g %20.16g\n", ia[i], ja[i], val[2 * i], val[2 * i + 1]);
318 else {
319 if (f != stdout) fclose(f);
320 return MM_UNSUPPORTED_TYPE;
321 }
322
323 if (f != stdout) fclose(f);
324
325 return 0;
326 }
327
mm_typecode_to_str(MM_typecode matcode)328 char *mm_typecode_to_str(MM_typecode matcode)
329 {
330 const char *types[4];
331
332 /* check for MTX type */
333 if (mm_is_matrix(matcode)) types[0] = MM_MTX_STR;
334 else return NULL;
335
336 /* check for CRD or ARR matrix */
337 if (mm_is_sparse(matcode)) types[1] = MM_SPARSE_STR;
338 else if (mm_is_dense(matcode)) types[1] = MM_DENSE_STR;
339 else return NULL;
340
341 /* check for element data type */
342 if (mm_is_real(matcode)) types[2] = MM_REAL_STR;
343 else if (mm_is_complex(matcode)) types[2] = MM_COMPLEX_STR;
344 else if (mm_is_pattern(matcode)) types[2] = MM_PATTERN_STR;
345 else if (mm_is_integer(matcode)) types[2] = MM_INT_STR;
346 else return NULL;
347
348 /* check for symmetry type */
349 if (mm_is_general(matcode)) types[3] = MM_GENERAL_STR;
350 else if (mm_is_symmetric(matcode)) types[3] = MM_SYMM_STR;
351 else if (mm_is_hermitian(matcode)) types[3] = MM_HERM_STR;
352 else if (mm_is_skew(matcode)) types[3] = MM_SKEW_STR;
353 else return NULL;
354
355 mm_buffer[0] = 0;
356 snprintf(mm_buffer, sizeof(mm_buffer) / sizeof(mm_buffer[0]), "%s %s %s %s", types[0], types[1], types[2], types[3]);
357 return mm_buffer;
358 }
359