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

Annotation of /trunk/mmio/mmio.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 993 - (hide annotations) (download) (as text)
Fri Dec 22 11:03:51 2006 UTC (15 years, 9 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 johnpye 993 /*
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