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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1636 - (show annotations)
Mon Jul 14 03:40:22 2008 UTC (11 years, 2 months ago) by trankine
File MIME type: text/plain
File size: 12966 byte(s)
somehow, this was not correctly merged from the branch??
1
2 /* $Id$ */
3
4 /*******************************************************
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
16 /*
17 * 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 #include "Common.h"
25 #include "mmio.h"
26
27 #include <string.h>
28 #include <ctype.h>
29
30
31
32 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 int *Ip, *Jp;
41
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 Ip = MEMALLOC(nz, int);
79 Jp = MEMALLOC(nz, int);
80 val = MEMALLOC(nz, double);
81
82 *val_ = val;
83 *I_ = Ip;
84 *J_ = Jp;
85
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 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 }
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 /* use when Ip[], Jp[], and val[] are already allocated */
279 /******************************************************************/
280
281 int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int Ip[], int Jp[],
282 double val[], MM_typecode matcode)
283 {
284 int i;
285 if (mm_is_complex(matcode))
286 {
287 for (i=0; i<nz; i++)
288 if (fscanf(f, "%d %d %lg %lg", &Ip[i], &Jp[i], &val[2*i], &val[2*i+1])
289 != 4) return MM_PREMATURE_EOF;
290 }
291 else if (mm_is_real(matcode))
292 {
293 for (i=0; i<nz; i++)
294 {
295 if (fscanf(f, "%d %d %lg\n", &Ip[i], &Jp[i], &val[i])
296 != 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 if (fscanf(f, "%d %d", &Ip[i], &Jp[i])
305 != 2) return MM_PREMATURE_EOF;
306 }
307 else
308 return MM_UNSUPPORTED_TYPE;
309
310 return 0;
311
312 }
313
314 int mm_read_mtx_crd_entry(FILE *f, int *Ip, int *Jp,
315 double *real, double *imag, MM_typecode matcode)
316 {
317 if (mm_is_complex(matcode))
318 {
319 if (fscanf(f, "%d %d %lg %lg", Ip, Jp, real, imag)
320 != 4) return MM_PREMATURE_EOF;
321 }
322 else if (mm_is_real(matcode))
323 {
324 if (fscanf(f, "%d %d %lg\n", Ip, Jp, real)
325 != 3) return MM_PREMATURE_EOF;
326
327 }
328
329 else if (mm_is_pattern(matcode))
330 {
331 if (fscanf(f, "%d %d", Ip, Jp) != 2) return MM_PREMATURE_EOF;
332 }
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 int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **Ip, int **Jp,
350 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 *Ip = MEMALLOC(*nz, int);
373 *Jp = MEMALLOC(*nz, int);
374 *val = NULL;
375
376 if (mm_is_complex(*matcode))
377 {
378 *val = MEMALLOC(*nz * 2, double);
379 ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *Ip, *Jp, *val,
380 *matcode);
381 if (ret_code != 0) return ret_code;
382 }
383 else if (mm_is_real(*matcode))
384 {
385 *val = MEMALLOC(*nz, double);
386 ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *Ip, *Jp, *val,
387 *matcode);
388 if (ret_code != 0) return ret_code;
389 }
390
391 else if (mm_is_pattern(*matcode))
392 {
393 ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *Ip, *Jp, *val,
394 *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 int mm_write_mtx_crd(char fname[], int M, int N, int nz, int Ip[], int Jp[],
412 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 fprintf(f, "%d %d\n", Ip[i], Jp[i]);
434 else
435 if (mm_is_real(matcode))
436 for (i=0; i<nz; i++)
437 fprintf(f, "%d %d %20.16g\n", Ip[i], Jp[i], val[i]);
438 else
439 if (mm_is_complex(matcode))
440 for (i=0; i<nz; i++)
441 fprintf(f, "%d %d %20.16g %20.16g\n", Ip[i], Jp[i], val[2*i],
442 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