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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26