/[escript]/trunk/paso/src/mmio.c
ViewVC logotype

Annotation of /trunk/paso/src/mmio.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 4 months ago) by ksteube
File MIME type: text/plain
File size: 12931 byte(s)
Copyright updated in all files

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26