/[escript]/branches/doubleplusgood/paso/src/mmio.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/mmio.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26