Actual source code: mmio.c

  1: /*
  2:    Matrix Market I/O library for ANSI C

  4:    See https://math.nist.gov/MatrixMarket/ for details.
  5: */

  7: #include <stdlib.h>
  8: #include <stdio.h>
  9: #include <string.h>
 10: #include <ctype.h>

 12: #include "mmio.h"

 14: static char mm_buffer[MM_MAX_LINE_LENGTH];

 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;

 25:   if ((f = fopen(fname, "r")) == NULL) return -1;

 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:   }

 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:   }

 39:   /* find out size of sparse matrix: M, N, nz .... */

 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:   }

 46:   *M_  = M;
 47:   *N_  = N;
 48:   *nz_ = nz;

 50:   /* reserve memory for matrices */

 52:   ia  = (int *)malloc(nz * sizeof(int));
 53:   ja  = (int *)malloc(nz * sizeof(int));
 54:   val = (double *)malloc(nz * sizeof(double));

 56:   *val_ = val;
 57:   *I_   = ia;
 58:   *J_   = ja;

 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)            */

 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);

 74:   return 0;
 75: }

 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: }

 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;

 96:   mm_clear_typecode(matcode);

 98:   if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF;

100:   if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type, storage_scheme) != 5) return MM_PREMATURE_EOF;

102:   for (p = mtx; *p != '\0'; *p = tolower(*p), p++)
103:     ; /* convert to lower case */
104:   for (p = crd; *p != '\0'; *p = tolower(*p), p++)
105:     ;
106:   for (p = data_type; *p != '\0'; *p = tolower(*p), p++)
107:     ;
108:   for (p = storage_scheme; *p != '\0'; *p = tolower(*p), p++)
109:     ;

111:   /* check for banner */
112:   if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0) return MM_NO_HEADER;

114:   /* first field should be "mtx" */
115:   if (strcmp(mtx, MM_MTX_STR) != 0) return MM_UNSUPPORTED_TYPE;
116:   mm_set_matrix(matcode);

118:   /* second field describes whether this is a sparse matrix (in coordinate
119:             storgae) or a dense array */

121:   if (strcmp(crd, MM_SPARSE_STR) == 0) mm_set_sparse(matcode);
122:   else if (strcmp(crd, MM_DENSE_STR) == 0) mm_set_dense(matcode);
123:   else return MM_UNSUPPORTED_TYPE;

125:   /* third field */

127:   if (strcmp(data_type, MM_REAL_STR) == 0) mm_set_real(matcode);
128:   else if (strcmp(data_type, MM_COMPLEX_STR) == 0) mm_set_complex(matcode);
129:   else if (strcmp(data_type, MM_PATTERN_STR) == 0) mm_set_pattern(matcode);
130:   else if (strcmp(data_type, MM_INT_STR) == 0) mm_set_integer(matcode);
131:   else return MM_UNSUPPORTED_TYPE;

133:   /* fourth field */

135:   if (strcmp(storage_scheme, MM_GENERAL_STR) == 0) mm_set_general(matcode);
136:   else if (strcmp(storage_scheme, MM_SYMM_STR) == 0) mm_set_symmetric(matcode);
137:   else if (strcmp(storage_scheme, MM_HERM_STR) == 0) mm_set_hermitian(matcode);
138:   else if (strcmp(storage_scheme, MM_SKEW_STR) == 0) mm_set_skew(matcode);
139:   else return MM_UNSUPPORTED_TYPE;

141:   return 0;
142: }

144: int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
145: {
146:   if (fprintf(f, "%d %d %d\n", M, N, nz) < 0) return MM_COULD_NOT_WRITE_FILE;
147:   else return 0;
148: }

150: int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz)
151: {
152:   char line[MM_MAX_LINE_LENGTH];
153:   int  num_items_read;

155:   /* set return null parameter values, in case we exit with errors */
156:   *M = *N = *nz = 0;

158:   /* now continue scanning until you reach the end-of-comments */
159:   do {
160:     if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF;
161:   } while (line[0] == '%');

163:   /* line[] is either blank or has M,N, nz */
164:   if (sscanf(line, "%d %d %d", M, N, nz) == 3) return 0;

166:   else do {
167:       num_items_read = fscanf(f, "%d %d %d", M, N, nz);
168:       if (num_items_read == EOF) return MM_PREMATURE_EOF;
169:     } while (num_items_read != 3);

171:   return 0;
172: }

174: int mm_read_mtx_array_size(FILE *f, int *M, int *N)
175: {
176:   char line[MM_MAX_LINE_LENGTH];
177:   int  num_items_read;
178:   /* set return null parameter values, in case we exit with errors */
179:   *M = *N = 0;

181:   /* now continue scanning until you reach the end-of-comments */
182:   do {
183:     if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) return MM_PREMATURE_EOF;
184:   } while (line[0] == '%');

186:   /* line[] is either blank or has M,N, nz */
187:   if (sscanf(line, "%d %d", M, N) == 2) return 0;

189:   else /* we have a blank line */ do {
190:       num_items_read = fscanf(f, "%d %d", M, N);
191:       if (num_items_read == EOF) return MM_PREMATURE_EOF;
192:     } while (num_items_read != 2);

194:   return 0;
195: }

197: int mm_write_mtx_array_size(FILE *f, int M, int N)
198: {
199:   if (fprintf(f, "%d %d\n", M, N) < 0) return MM_COULD_NOT_WRITE_FILE;
200:   else return 0;
201: }

203: /*-------------------------------------------------------------------------*/

205: /******************************************************************/
206: /* use when ia[], ja[], and val[] are already allocated */
207: /******************************************************************/

209: int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int ia[], int ja[], double val[], MM_typecode matcode)
210: {
211:   int i;
212:   if (mm_is_complex(matcode)) {
213:     for (i = 0; i < nz; i++)
214:       if (fscanf(f, "%d %d %lg %lg", &ia[i], &ja[i], &val[2 * i], &val[2 * i + 1]) != 4) return MM_PREMATURE_EOF;
215:   } else if (mm_is_real(matcode)) {
216:     for (i = 0; i < nz; i++) {
217:       if (fscanf(f, "%d %d %lg\n", &ia[i], &ja[i], &val[i]) != 3) return MM_PREMATURE_EOF;
218:     }
219:   }

221:   else if (mm_is_pattern(matcode)) {
222:     for (i = 0; i < nz; i++)
223:       if (fscanf(f, "%d %d", &ia[i], &ja[i]) != 2) return MM_PREMATURE_EOF;
224:   } else return MM_UNSUPPORTED_TYPE;

226:   return 0;
227: }

229: int mm_read_mtx_crd_entry(FILE *f, int *ia, int *ja, double *real, double *imag, MM_typecode matcode)
230: {
231:   if (mm_is_complex(matcode)) {
232:     if (fscanf(f, "%d %d %lg %lg", ia, ja, real, imag) != 4) return MM_PREMATURE_EOF;
233:   } else if (mm_is_real(matcode)) {
234:     if (fscanf(f, "%d %d %lg\n", ia, ja, real) != 3) return MM_PREMATURE_EOF;

236:   }

238:   else if (mm_is_pattern(matcode)) {
239:     if (fscanf(f, "%d %d", ia, ja) != 2) return MM_PREMATURE_EOF;
240:   } else return MM_UNSUPPORTED_TYPE;

242:   return 0;
243: }

245: /************************************************************************
246:     mm_read_mtx_crd()  fills M, N, nz, array of values, and return
247:                         type code, e.g. 'MCRS'

249:                         if matrix is complex, values[] is of size 2*nz,
250:                             (nz pairs of real/imaginary values)
251: ************************************************************************/

253: int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **ia, int **ja, double **val, MM_typecode *matcode)
254: {
255:   int   ret_code;
256:   FILE *f;

258:   if (strcmp(fname, "stdin") == 0) f = stdin;
259:   else if ((f = fopen(fname, "r")) == NULL) return MM_COULD_NOT_READ_FILE;

261:   if ((ret_code = mm_read_banner(f, matcode)) != 0) return ret_code;

263:   if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) && mm_is_matrix(*matcode))) return MM_UNSUPPORTED_TYPE;

265:   if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0) return ret_code;

267:   *ia  = (int *)malloc(*nz * sizeof(int));
268:   *ja  = (int *)malloc(*nz * sizeof(int));
269:   *val = NULL;

271:   if (mm_is_complex(*matcode)) {
272:     *val     = (double *)malloc(*nz * 2 * 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:   } else if (mm_is_real(*matcode)) {
276:     *val     = (double *)malloc(*nz * sizeof(double));
277:     ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *ia, *ja, *val, *matcode);
278:     if (ret_code != 0) return ret_code;
279:   }

281:   else if (mm_is_pattern(*matcode)) {
282:     ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *ia, *ja, *val, *matcode);
283:     if (ret_code != 0) return ret_code;
284:   }

286:   if (f != stdin) fclose(f);
287:   return 0;
288: }

290: int mm_write_banner(FILE *f, MM_typecode matcode)
291: {
292:   char *str = mm_typecode_to_str(matcode);
293:   int   ret_code;

295:   ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str);
296:   if (ret_code < 0) return MM_COULD_NOT_WRITE_FILE;
297:   else return 0;
298: }

300: int mm_write_mtx_crd(char fname[], int M, int N, int nz, int ia[], int ja[], double val[], MM_typecode matcode)
301: {
302:   FILE *f;
303:   int   i;

305:   if (strcmp(fname, "stdout") == 0) f = stdout;
306:   else if ((f = fopen(fname, "w")) == NULL) return MM_COULD_NOT_WRITE_FILE;

308:   /* print banner followed by typecode */
309:   fprintf(f, "%s ", MatrixMarketBanner);
310:   fprintf(f, "%s\n", mm_typecode_to_str(matcode));

312:   /* print matrix sizes and nonzeros */
313:   fprintf(f, "%d %d %d\n", M, N, nz);

315:   /* print values */
316:   if (mm_is_pattern(matcode))
317:     for (i = 0; i < nz; i++) fprintf(f, "%d %d\n", ia[i], ja[i]);
318:   else if (mm_is_real(matcode))
319:     for (i = 0; i < nz; i++) fprintf(f, "%d %d %20.16g\n", ia[i], ja[i], val[i]);
320:   else if (mm_is_complex(matcode))
321:     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]);
322:   else {
323:     if (f != stdout) fclose(f);
324:     return MM_UNSUPPORTED_TYPE;
325:   }

327:   if (f != stdout) fclose(f);

329:   return 0;
330: }

332: char *mm_typecode_to_str(MM_typecode matcode)
333: {
334:   const char *types[4];

336:   /* check for MTX type */
337:   if (mm_is_matrix(matcode)) types[0] = MM_MTX_STR;
338:   else return NULL;

340:   /* check for CRD or ARR matrix */
341:   if (mm_is_sparse(matcode)) types[1] = MM_SPARSE_STR;
342:   else if (mm_is_dense(matcode)) types[1] = MM_DENSE_STR;
343:   else return NULL;

345:   /* check for element data type */
346:   if (mm_is_real(matcode)) types[2] = MM_REAL_STR;
347:   else if (mm_is_complex(matcode)) types[2] = MM_COMPLEX_STR;
348:   else if (mm_is_pattern(matcode)) types[2] = MM_PATTERN_STR;
349:   else if (mm_is_integer(matcode)) types[2] = MM_INT_STR;
350:   else return NULL;

352:   /* check for symmetry type */
353:   if (mm_is_general(matcode)) types[3] = MM_GENERAL_STR;
354:   else if (mm_is_symmetric(matcode)) types[3] = MM_SYMM_STR;
355:   else if (mm_is_hermitian(matcode)) types[3] = MM_HERM_STR;
356:   else if (mm_is_skew(matcode)) types[3] = MM_SKEW_STR;
357:   else return NULL;

359:   mm_buffer[0] = 0;
360:   snprintf(mm_buffer, sizeof(mm_buffer) / sizeof(mm_buffer[0]), "%s %s %s %s", types[0], types[1], types[2], types[3]);
361:   return mm_buffer;
362: }