/[ascend]/trunk/mmio/mmio.c
ViewVC logotype

Contents of /trunk/mmio/mmio.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 993 - (show annotations) (download) (as text)
Fri Dec 22 11:03:51 2006 UTC (13 years, 5 months ago) by johnpye
File MIME type: text/x-csrc
File size: 12870 byte(s)
Added Matrix Market export routines (from math.nist.gov).
Added necessary build commands for above.
Work on Jacobi preconditioner for IDA (ongoing)
Set integrator_analye to make a call to slv_block_partition. Not sure if that's a good idea or not.
Tidied up comments in linsol.h
Moved unnecessary #defines from model_reorder.h into model_reorder.c.
Cleaned up codedocs in mtx*.h
Added WITH_MMIO config.h flag.

1 /*
2 * Matrix Market I/O library for ANSI C
3 *
4 * See http://math.nist.gov/MatrixMarket for details.
5 *
6 *
7 */
8
9
10 #include <stdio.h>
11 #include <string.h>
12 #include <malloc.h>
13 #include <ctype.h>
14
15 #include "mmio.h"
16
17 int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
18 double **val_, int **I_, int **J_)
19 {
20 FILE *f;
21 MM_typecode matcode;
22 int M, N, nz;
23 int i;
24 double *val;
25 int *I, *J;
26
27 if ((f = fopen(fname, "r")) == NULL)
28 return -1;
29
30
31 if (mm_read_banner(f, &matcode) != 0)
32 {
33 printf("mm_read_unsymetric: Could not process Matrix Market banner ");
34 printf(" in file [%s]\n", fname);
35 return -1;
36 }
37
38
39
40 if ( !(mm_is_real(matcode) && mm_is_matrix(matcode) &&
41 mm_is_sparse(matcode)))
42 {
43 fprintf(stderr, "Sorry, this application does not support ");
44 fprintf(stderr, "Market Market type: [%s]\n",
45 mm_typecode_to_str(matcode));
46 return -1;
47 }
48
49 /* find out size of sparse matrix: M, N, nz .... */
50
51 if (mm_read_mtx_crd_size(f, &M, &N, &nz) !=0)
52 {
53 fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n");
54 return -1;
55 }
56
57 *M_ = M;
58 *N_ = N;
59 *nz_ = nz;
60
61 /* reseve memory for matrices */
62
63 I = (int *) malloc(nz * sizeof(int));
64 J = (int *) malloc(nz * sizeof(int));
65 val = (double *) malloc(nz * sizeof(double));
66
67 *val_ = val;
68 *I_ = I;
69 *J_ = J;
70
71 /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
72 /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
73 /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
74
75 for (i=0; i<nz; i++)
76 {
77 fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]);
78 I[i]--; /* adjust from 1-based to 0-based */
79 J[i]--;
80 }
81 fclose(f);
82
83 return 0;
84 }
85
86 int mm_is_valid(MM_typecode matcode)
87 {
88 if (!mm_is_matrix(matcode)) return 0;
89 if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
90 if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
91 if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) ||
92 mm_is_skew(matcode))) return 0;
93 return 1;
94 }
95
96 int mm_read_banner(FILE *f, MM_typecode *matcode)
97 {
98 char line[MM_MAX_LINE_LENGTH];
99 char banner[MM_MAX_TOKEN_LENGTH];
100 char mtx[MM_MAX_TOKEN_LENGTH];
101 char crd[MM_MAX_TOKEN_LENGTH];
102 char data_type[MM_MAX_TOKEN_LENGTH];
103 char storage_scheme[MM_MAX_TOKEN_LENGTH];
104 char *p;
105
106
107 mm_clear_typecode(matcode);
108
109 if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL)
110 return MM_PREMATURE_EOF;
111
112 if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type,
113 storage_scheme) != 5)
114 return MM_PREMATURE_EOF;
115
116 for (p=mtx; *p!='\0'; *p=tolower(*p),p++); /* convert to lower case */
117 for (p=crd; *p!='\0'; *p=tolower(*p),p++);
118 for (p=data_type; *p!='\0'; *p=tolower(*p),p++);
119 for (p=storage_scheme; *p!='\0'; *p=tolower(*p),p++);
120
121 /* check for banner */
122 if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
123 return MM_NO_HEADER;
124
125 /* first field should be "mtx" */
126 if (strcmp(mtx, MM_MTX_STR) != 0)
127 return MM_UNSUPPORTED_TYPE;
128 mm_set_matrix(matcode);
129
130
131 /* second field describes whether this is a sparse matrix (in coordinate
132 storgae) or a dense array */
133
134
135 if (strcmp(crd, MM_SPARSE_STR) == 0)
136 mm_set_sparse(matcode);
137 else
138 if (strcmp(crd, MM_DENSE_STR) == 0)
139 mm_set_dense(matcode);
140 else
141 return MM_UNSUPPORTED_TYPE;
142
143
144 /* third field */
145
146 if (strcmp(data_type, MM_REAL_STR) == 0)
147 mm_set_real(matcode);
148 else
149 if (strcmp(data_type, MM_COMPLEX_STR) == 0)
150 mm_set_complex(matcode);
151 else
152 if (strcmp(data_type, MM_PATTERN_STR) == 0)
153 mm_set_pattern(matcode);
154 else
155 if (strcmp(data_type, MM_INT_STR) == 0)
156 mm_set_integer(matcode);
157 else
158 return MM_UNSUPPORTED_TYPE;
159
160
161 /* fourth field */
162
163 if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
164 mm_set_general(matcode);
165 else
166 if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
167 mm_set_symmetric(matcode);
168 else
169 if (strcmp(storage_scheme, MM_HERM_STR) == 0)
170 mm_set_hermitian(matcode);
171 else
172 if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
173 mm_set_skew(matcode);
174 else
175 return MM_UNSUPPORTED_TYPE;
176
177
178 return 0;
179 }
180
181 int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
182 {
183 if (fprintf(f, "%d %d %d\n", M, N, nz) != 3)
184 return MM_COULD_NOT_WRITE_FILE;
185 else
186 return 0;
187 }
188
189 int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
190 {
191 char line[MM_MAX_LINE_LENGTH];
192 int num_items_read;
193
194 /* set return null parameter values, in case we exit with errors */
195 *M = *N = *nz = 0;
196
197 /* now continue scanning until you reach the end-of-comments */
198 do
199 {
200 if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
201 return MM_PREMATURE_EOF;
202 }while (line[0] == '%');
203
204 /* line[] is either blank or has M,N, nz */
205 if (sscanf(line, "%d %d %d", M, N, nz) == 3)
206 return 0;
207
208 else
209 do
210 {
211 num_items_read = fscanf(f, "%d %d %d", M, N, nz);
212 if (num_items_read == EOF) return MM_PREMATURE_EOF;
213 }
214 while (num_items_read != 3);
215
216 return 0;
217 }
218
219
220 int mm_read_mtx_array_size(FILE *f, int *M, int *N)
221 {
222 char line[MM_MAX_LINE_LENGTH];
223 int num_items_read;
224 /* set return null parameter values, in case we exit with errors */
225 *M = *N = 0;
226
227 /* now continue scanning until you reach the end-of-comments */
228 do
229 {
230 if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
231 return MM_PREMATURE_EOF;
232 }while (line[0] == '%');
233
234 /* line[] is either blank or has M,N, nz */
235 if (sscanf(line, "%d %d", M, N) == 2)
236 return 0;
237
238 else /* we have a blank line */
239 do
240 {
241 num_items_read = fscanf(f, "%d %d", M, N);
242 if (num_items_read == EOF) return MM_PREMATURE_EOF;
243 }
244 while (num_items_read != 2);
245
246 return 0;
247 }
248
249 int mm_write_mtx_array_size(FILE *f, int M, int N)
250 {
251 if (fprintf(f, "%d %d\n", M, N) != 2)
252 return MM_COULD_NOT_WRITE_FILE;
253 else
254 return 0;
255 }
256
257
258
259 /*-------------------------------------------------------------------------*/
260
261 /******************************************************************/
262 /* use when I[], J[], and val[]J, and val[] are already allocated */
263 /******************************************************************/
264
265 int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
266 double val[], MM_typecode matcode)
267 {
268 int i;
269 if (mm_is_complex(matcode))
270 {
271 for (i=0; i<nz; i++)
272 if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2*i], &val[2*i+1])
273 != 4) return MM_PREMATURE_EOF;
274 }
275 else if (mm_is_real(matcode))
276 {
277 for (i=0; i<nz; i++)
278 {
279 if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i])
280 != 3) return MM_PREMATURE_EOF;
281
282 }
283 }
284
285 else if (mm_is_pattern(matcode))
286 {
287 for (i=0; i<nz; i++)
288 if (fscanf(f, "%d %d", &I[i], &J[i])
289 != 2) return MM_PREMATURE_EOF;
290 }
291 else
292 return MM_UNSUPPORTED_TYPE;
293
294 return 0;
295
296 }
297
298 int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
299 double *real, double *imag, MM_typecode matcode)
300 {
301 if (mm_is_complex(matcode))
302 {
303 if (fscanf(f, "%d %d %lg %lg", I, J, real, imag)
304 != 4) return MM_PREMATURE_EOF;
305 }
306 else if (mm_is_real(matcode))
307 {
308 if (fscanf(f, "%d %d %lg\n", I, J, real)
309 != 3) return MM_PREMATURE_EOF;
310
311 }
312
313 else if (mm_is_pattern(matcode))
314 {
315 if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF;
316 }
317 else
318 return MM_UNSUPPORTED_TYPE;
319
320 return 0;
321
322 }
323
324
325 /************************************************************************
326 mm_read_mtx_crd() fills M, N, nz, array of values, and return
327 type code, e.g. 'MCRS'
328
329 if matrix is complex, values[] is of size 2*nz,
330 (nz pairs of real/imaginary values)
331 ************************************************************************/
332
333 int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J,
334 double **val, MM_typecode *matcode)
335 {
336 int ret_code;
337 FILE *f;
338
339 if (strcmp(fname, "stdin") == 0) f=stdin;
340 else
341 if ((f = fopen(fname, "r")) == NULL)
342 return MM_COULD_NOT_READ_FILE;
343
344
345 if ((ret_code = mm_read_banner(f, matcode)) != 0)
346 return ret_code;
347
348 if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) &&
349 mm_is_matrix(*matcode)))
350 return MM_UNSUPPORTED_TYPE;
351
352 if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
353 return ret_code;
354
355
356 *I = (int *) malloc(*nz * sizeof(int));
357 *J = (int *) malloc(*nz * sizeof(int));
358 *val = NULL;
359
360 if (mm_is_complex(*matcode))
361 {
362 *val = (double *) malloc(*nz * 2 * sizeof(double));
363 ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
364 *matcode);
365 if (ret_code != 0) return ret_code;
366 }
367 else if (mm_is_real(*matcode))
368 {
369 *val = (double *) malloc(*nz * sizeof(double));
370 ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
371 *matcode);
372 if (ret_code != 0) return ret_code;
373 }
374
375 else if (mm_is_pattern(*matcode))
376 {
377 ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
378 *matcode);
379 if (ret_code != 0) return ret_code;
380 }
381
382 if (f != stdin) fclose(f);
383 return 0;
384 }
385
386 int mm_write_banner(FILE *f, MM_typecode matcode)
387 {
388 char *str = mm_typecode_to_str(matcode);
389 int ret_code;
390
391 ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str);
392 free(str);
393 if (ret_code !=2 )
394 return MM_COULD_NOT_WRITE_FILE;
395 else
396 return 0;
397 }
398
399 int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
400 double val[], MM_typecode matcode)
401 {
402 FILE *f;
403 int i;
404
405 if (strcmp(fname, "stdout") == 0)
406 f = stdout;
407 else
408 if ((f = fopen(fname, "w")) == NULL)
409 return MM_COULD_NOT_WRITE_FILE;
410
411 /* print banner followed by typecode */
412 fprintf(f, "%s ", MatrixMarketBanner);
413 fprintf(f, "%s\n", mm_typecode_to_str(matcode));
414
415 /* print matrix sizes and nonzeros */
416 fprintf(f, "%d %d %d\n", M, N, nz);
417
418 /* print values */
419 if (mm_is_pattern(matcode))
420 for (i=0; i<nz; i++)
421 fprintf(f, "%d %d\n", I[i], J[i]);
422 else
423 if (mm_is_real(matcode))
424 for (i=0; i<nz; i++)
425 fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]);
426 else
427 if (mm_is_complex(matcode))
428 for (i=0; i<nz; i++)
429 fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2*i],
430 val[2*i+1]);
431 else
432 {
433 if (f != stdout) fclose(f);
434 return MM_UNSUPPORTED_TYPE;
435 }
436
437 if (f !=stdout) fclose(f);
438
439 return 0;
440 }
441
442
443 /**
444 * Create a new copy of a string s. strdup() is a common routine, but
445 * not part of ANSI C, so it is included here. Used by mm_typecode_to_str().
446 *
447 */
448 char *strdup(const char *s)
449 {
450 int len = strlen(s);
451 char *s2 = (char *) malloc((len+1)*sizeof(char));
452 return strcpy(s2, s);
453 }
454
455 char *mm_typecode_to_str(MM_typecode matcode)
456 {
457 char buffer[MM_MAX_LINE_LENGTH];
458 char *types[4];
459 char *strdup(const char *);
460 int error =0;
461
462 /* check for MTX type */
463 if (mm_is_matrix(matcode))
464 types[0] = MM_MTX_STR;
465 else
466 error=1;
467
468 /* check for CRD or ARR matrix */
469 if (mm_is_sparse(matcode))
470 types[1] = MM_SPARSE_STR;
471 else
472 if (mm_is_dense(matcode))
473 types[1] = MM_DENSE_STR;
474 else
475 return NULL;
476
477 /* check for element data type */
478 if (mm_is_real(matcode))
479 types[2] = MM_REAL_STR;
480 else
481 if (mm_is_complex(matcode))
482 types[2] = MM_COMPLEX_STR;
483 else
484 if (mm_is_pattern(matcode))
485 types[2] = MM_PATTERN_STR;
486 else
487 if (mm_is_integer(matcode))
488 types[2] = MM_INT_STR;
489 else
490 return NULL;
491
492
493 /* check for symmetry type */
494 if (mm_is_general(matcode))
495 types[3] = MM_GENERAL_STR;
496 else
497 if (mm_is_symmetric(matcode))
498 types[3] = MM_SYMM_STR;
499 else
500 if (mm_is_hermitian(matcode))
501 types[3] = MM_HERM_STR;
502 else
503 if (mm_is_skew(matcode))
504 types[3] = MM_SKEW_STR;
505 else
506 return NULL;
507
508 sprintf(buffer,"%s %s %s %s", types[0], types[1], types[2], types[3]);
509 return strdup(buffer);
510
511 }

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22