EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_mmio.cpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ***********************************************************************
4 //
5 // EpetraExt: Epetra Extended - Linear Algebra Services Package
6 // Copyright (2011) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //@HEADER
42 */
43 
44 /*
45 * Matrix Market I/O library for ANSI C
46 *
47 * See http://math.nist.gov/MatrixMarket for details.
48 *
49 *
50 */
51 
52 /* JW - This is essentially a C file that was "converted" to C++, but still
53  contains C function calls. Since it is independent of the rest of
54  the package, we will not include the central ConfigDefs file, but
55  will rather #include the C headers that we need. This should fix the
56  build problem on sass2889.
57 #include <Epetra_ConfigDefs.h> */
58 #include <string.h>
59 #include <stdio.h>
60 #include <ctype.h>
61 #include "Teuchos_Assert.hpp"
62 #include "EpetraExt_mmio.h"
63 
64 using namespace EpetraExt;
65 namespace EpetraExt {
66 int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
67  double **val_, int **I_, int **J_)
68 {
69  FILE *f;
70  MM_typecode matcode;
71  int M, N, nz;
72  int i;
73  double *val;
74  int *I, *J;
75 
76  if ((f = fopen(fname, "r")) == NULL)
77  return -1;
78 
79 
80  if (mm_read_banner(f, &matcode) != 0)
81  {
82  printf("mm_read_unsymetric: Could not process Matrix Market banner ");
83  printf(" in file [%s]\n", fname);
84  return -1;
85  }
86 
87 
88 
89  if ( !(mm_is_real(matcode) && mm_is_matrix(matcode) &&
90  mm_is_sparse(matcode)))
91  {
92  char buffer[MM_MAX_LINE_LENGTH];
93  mm_typecode_to_str(matcode, buffer);
94  fprintf(stderr, "Sorry, this application does not support ");
95  fprintf(stderr, "Market Market type: [%s]\n",buffer);
96  return -1;
97  }
98 
99  /* find out size of sparse matrix: M, N, nz .... */
100 
101  if (mm_read_mtx_crd_size(f, &M, &N, &nz) !=0)
102  {
103  fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n");
104  return -1;
105  }
106 
107  *M_ = M;
108  *N_ = N;
109  *nz_ = nz;
110 
111  /* reseve memory for matrices */
112 
113  //I = (int *) malloc(nz * sizeof(int));
114  //J = (int *) malloc(nz * sizeof(int));
115  //val = (double *) malloc(nz * sizeof(double));
116 
117  I = new int[nz];
118  J = new int[nz];
119  val = new double[nz];
120 
121  *val_ = val;
122  *I_ = I;
123  *J_ = J;
124 
125  /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
126  /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
127  /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
128 
129  for (i=0; i<nz; i++)
130  {
131  TEUCHOS_ASSERT(fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]) != EOF);
132  I[i]--; /* adjust from 1-based to 0-based */
133  J[i]--;
134  }
135  fclose(f);
136 
137  return 0;
138 }
139 
141 {
142  if (!mm_is_matrix(matcode)) return 0;
143  if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
144  if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
145  if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) ||
146  mm_is_skew(matcode))) return 0;
147  return 1;
148 }
149 
150 int mm_read_banner(FILE *f, MM_typecode *matcode)
151 {
152  char line[MM_MAX_LINE_LENGTH];
153  char banner[MM_MAX_TOKEN_LENGTH];
154  char mtx[MM_MAX_TOKEN_LENGTH];
155  char crd[MM_MAX_TOKEN_LENGTH];
156  char data_type[MM_MAX_TOKEN_LENGTH];
157  char storage_scheme[MM_MAX_TOKEN_LENGTH];
158  char *p;
159 
160 
161  mm_clear_typecode(matcode);
162 
163  if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL)
164  return MM_PREMATURE_EOF;
165 
166  if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type,
167  storage_scheme) != 5)
168  return MM_PREMATURE_EOF;
169 
170  for (p=mtx; *p!='\0'; *p=tolower(*p),p++); /* convert to lower case */
171  for (p=crd; *p!='\0'; *p=tolower(*p),p++);
172  for (p=data_type; *p!='\0'; *p=tolower(*p),p++);
173  for (p=storage_scheme; *p!='\0'; *p=tolower(*p),p++);
174 
175  /* check for banner */
176  if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
177  return MM_NO_HEADER;
178 
179  /* first field should be "mtx" */
180  if (strcmp(mtx, MM_MTX_STR) != 0)
181  return MM_UNSUPPORTED_TYPE;
182  mm_set_matrix(matcode);
183 
184 
185  /* second field describes whether this is a sparse matrix (in coordinate
186  storgae) or a dense array */
187 
188 
189  if (strcmp(crd, MM_SPARSE_STR) == 0)
190  mm_set_sparse(matcode);
191  else
192  if (strcmp(crd, MM_DENSE_STR) == 0)
193  mm_set_dense(matcode);
194  else
195  return MM_UNSUPPORTED_TYPE;
196 
197 
198  /* third field */
199 
200  if (strcmp(data_type, MM_REAL_STR) == 0)
201  mm_set_real(matcode);
202  else
203  if (strcmp(data_type, MM_COMPLEX_STR) == 0)
204  mm_set_complex(matcode);
205  else
206  if (strcmp(data_type, MM_PATTERN_STR) == 0)
207  mm_set_pattern(matcode);
208  else
209  if (strcmp(data_type, MM_INT_STR) == 0)
210  mm_set_integer(matcode);
211  else
212  return MM_UNSUPPORTED_TYPE;
213 
214 
215  /* fourth field */
216 
217  if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
218  mm_set_general(matcode);
219  else
220  if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
221  mm_set_symmetric(matcode);
222  else
223  if (strcmp(storage_scheme, MM_HERM_STR) == 0)
224  mm_set_hermitian(matcode);
225  else
226  if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
227  mm_set_skew(matcode);
228  else
229  return MM_UNSUPPORTED_TYPE;
230 
231 
232  return 0;
233 }
234 
235 int mm_write_mtx_crd_size(FILE *f, long long M, long long N, long long nz)
236 {
237  fprintf(f, "%lld %lld %lld\n", M, N, nz);
238  return 0;
239 }
240 
241 int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
242 {
243  char line[MM_MAX_LINE_LENGTH];
244  int num_items_read;
245 
246  /* set return null parameter values, in case we exit with errors */
247  *M = *N = *nz = 0;
248 
249  /* now continue scanning until you reach the end-of-comments */
250  do
251  {
252  if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
253  return MM_PREMATURE_EOF;
254  }while (line[0] == '%');
255 
256  /* line[] is either blank or has M,N, nz */
257  if (sscanf(line, "%d %d %d", M, N, nz) == 3)
258  return 0;
259 
260  else
261  do
262  {
263  num_items_read = fscanf(f, "%d %d %d", M, N, nz);
264  if (num_items_read == EOF) return MM_PREMATURE_EOF;
265  }
266  while (num_items_read != 3);
267 
268  return 0;
269 }
270 
271 int mm_read_mtx_crd_size(FILE *f, long long *M, long long *N, long long *nz )
272 {
273  char line[MM_MAX_LINE_LENGTH];
274  int num_items_read;
275 
276  /* set return null parameter values, in case we exit with errors */
277  *M = *N = *nz = 0;
278 
279  /* now continue scanning until you reach the end-of-comments */
280  do
281  {
282  if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
283  return MM_PREMATURE_EOF;
284  }while (line[0] == '%');
285 
286  /* line[] is either blank or has M,N, nz */
287  if (sscanf(line, "%lld %lld %lld", M, N, nz) == 3)
288  return 0;
289 
290  else
291  do
292  {
293  num_items_read = fscanf(f, "%lld %lld %lld", M, N, nz);
294  if (num_items_read == EOF) return MM_PREMATURE_EOF;
295  }
296  while (num_items_read != 3);
297 
298  return 0;
299 }
300 
301 int mm_read_mtx_array_size(FILE *f, int *M, int *N)
302 {
303  char line[MM_MAX_LINE_LENGTH];
304  int num_items_read;
305  /* set return null parameter values, in case we exit with errors */
306  *M = *N = 0;
307 
308  /* now continue scanning until you reach the end-of-comments */
309  do
310  {
311  if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
312  return MM_PREMATURE_EOF;
313  }while (line[0] == '%');
314 
315  /* line[] is either blank or has M,N, nz */
316  if (sscanf(line, "%d %d", M, N) == 2)
317  return 0;
318 
319  else /* we have a blank line */
320  do
321  {
322  num_items_read = fscanf(f, "%d %d", M, N);
323  if (num_items_read == EOF) return MM_PREMATURE_EOF;
324  }
325  while (num_items_read != 2);
326 
327  return 0;
328 }
329 
330 int mm_write_mtx_array_size(FILE *f, long long M, long long N)
331 {
332  fprintf(f, "%lld %lld\n", M, N);
333  return 0;
334 }
335 
336 
337 
338 /*-------------------------------------------------------------------------*/
339 
340 /******************************************************************/
341 /* use when I[], J[], and val[]J, and val[] are already allocated */
342 /******************************************************************/
343 
344 int mm_read_mtx_crd_data(FILE *f, int /* M */, int /* N */, int nz, int I[], int J[],
345  double val[], MM_typecode matcode)
346 {
347  int i;
348  if (mm_is_complex(matcode))
349  {
350  for (i=0; i<nz; i++)
351  if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2*i], &val[2*i+1])
352  != 4) return MM_PREMATURE_EOF;
353  }
354  else if (mm_is_real(matcode))
355  {
356  for (i=0; i<nz; i++)
357  {
358  if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i])
359  != 3) return MM_PREMATURE_EOF;
360 
361  }
362  }
363 
364  else if (mm_is_pattern(matcode))
365  {
366  for (i=0; i<nz; i++)
367  if (fscanf(f, "%d %d", &I[i], &J[i])
368  != 2) return MM_PREMATURE_EOF;
369  }
370  else
371  return MM_UNSUPPORTED_TYPE;
372 
373  return 0;
374 
375 }
376 
377 int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
378  double *real, double *imag, MM_typecode matcode)
379 {
380  if (mm_is_complex(matcode))
381  {
382  if (fscanf(f, "%d %d %lg %lg", I, J, real, imag)
383  != 4) return MM_PREMATURE_EOF;
384  }
385  else if (mm_is_real(matcode))
386  {
387  if (fscanf(f, "%d %d %lg\n", I, J, real)
388  != 3) return MM_PREMATURE_EOF;
389 
390  }
391 
392  else if (mm_is_pattern(matcode))
393  {
394  if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF;
395  }
396  else
397  return MM_UNSUPPORTED_TYPE;
398 
399  return 0;
400 
401 }
402 
403 int mm_read_mtx_crd_entry(FILE *f, long long *I, long long *J,
404  double *real, double *imag, MM_typecode matcode)
405 {
406  if (mm_is_complex(matcode))
407  {
408  if (fscanf(f, "%lld %lld %lg %lg", I, J, real, imag)
409  != 4) return MM_PREMATURE_EOF;
410  }
411  else if (mm_is_real(matcode))
412  {
413  if (fscanf(f, "%lld %lld %lg\n", I, J, real)
414  != 3) return MM_PREMATURE_EOF;
415 
416  }
417 
418  else if (mm_is_pattern(matcode))
419  {
420  if (fscanf(f, "%lld %lld", I, J) != 2) return MM_PREMATURE_EOF;
421  }
422  else
423  return MM_UNSUPPORTED_TYPE;
424 
425  return 0;
426 
427 }
428 
429 /************************************************************************
430  mm_read_mtx_crd() fills M, N, nz, array of values, and return
431  type code, e.g. 'MCRS'
432 
433  if matrix is complex, values[] is of size 2*nz,
434  (nz pairs of real/imaginary values)
435 ************************************************************************/
436 
437 int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J,
438  double **val, MM_typecode *matcode)
439 {
440  int ret_code;
441  FILE *f;
442 
443  if (strcmp(fname, "stdin") == 0) f=stdin;
444  else
445  if ((f = fopen(fname, "r")) == NULL)
446  return MM_COULD_NOT_READ_FILE;
447 
448 
449  if ((ret_code = mm_read_banner(f, matcode)) != 0)
450  return ret_code;
451 
452  if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) &&
453  mm_is_matrix(*matcode)))
454  return MM_UNSUPPORTED_TYPE;
455 
456  if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
457  return ret_code;
458 
459 
460  //*I = (int *) malloc(*nz * sizeof(int));
461  //*J = (int *) malloc(*nz * sizeof(int));
462  //*val = NULL;
463 
464  *I = new int[*nz];
465  *J = new int[*nz];
466  *val = 0;
467 
468  if (mm_is_complex(*matcode))
469  {
470  //*val = (double *) malloc(*nz * 2 * sizeof(double));
471  *val = new double[2*(*nz)];
472  ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
473  *matcode);
474  if (ret_code != 0) return ret_code;
475  }
476  else if (mm_is_real(*matcode))
477  {
478  //*val = (double *) malloc(*nz * sizeof(double));
479  *val = new double[*nz];
480  ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
481  *matcode);
482  if (ret_code != 0) return ret_code;
483  }
484 
485  else if (mm_is_pattern(*matcode))
486  {
487  ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
488  *matcode);
489  if (ret_code != 0) return ret_code;
490  }
491 
492  if (f != stdin) fclose(f);
493  return 0;
494 }
495 
496 int mm_write_banner(FILE *f, MM_typecode matcode)
497 {
498  char buffer[MM_MAX_LINE_LENGTH];
499 
500  mm_typecode_to_str(matcode, buffer);
501  const int ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, buffer);
502 
503  if (ret_code < 0) {
504  return -1;
505  } else {
506  return 0;
507  }
508 }
509 
510 int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
511  double val[], MM_typecode matcode)
512 {
513  FILE *f;
514  int i;
515 
516  if (strcmp(fname, "stdout") == 0)
517  f = stdout;
518  else
519  if ((f = fopen(fname, "w")) == NULL)
521 
522  /* print banner followed by typecode */
523  char buffer[MM_MAX_LINE_LENGTH];
524  mm_typecode_to_str(matcode, buffer);
525  fprintf(f, "%s ", MatrixMarketBanner);
526  fprintf(f, "%s\n", buffer);
527 
528  /* print matrix sizes and nonzeros */
529  fprintf(f, "%d %d %d\n", M, N, nz);
530 
531  /* print values */
532  if (mm_is_pattern(matcode))
533  for (i=0; i<nz; i++)
534  fprintf(f, "%d %d\n", I[i], J[i]);
535  else
536  if (mm_is_real(matcode))
537  for (i=0; i<nz; i++)
538  fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]);
539  else
540  if (mm_is_complex(matcode))
541  for (i=0; i<nz; i++)
542  fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2*i],
543  val[2*i+1]);
544  else
545  {
546  if (f != stdout) fclose(f);
547  return MM_UNSUPPORTED_TYPE;
548  }
549 
550  if (f !=stdout) fclose(f);
551 
552  return 0;
553 }
554 
555 
556 void mm_typecode_to_str(MM_typecode matcode, char * buffer)
557 {
558  char type0[20];
559  char type1[20];
560  char type2[20];
561  char type3[20];
562  //int error =0; // unused
563 
564  /* check for MTX type */
565  if (mm_is_matrix(matcode))
566  strcpy(type0, MM_MTX_STR);
567  else {
568  // FIXME (mfh 27 Mar 2015) The code I found here set error=1,
569  // but kept going. It would make more sense to return
570  // immediately on failure, as the code below does. It would
571  // make even _more_ sense to return an error code, like some
572  // other functions in this file, but I don't want to change
573  // the public interface.
574  return;
575  //error=1; // unused
576  }
577 
578  /* check for CRD or ARR matrix */
579  if (mm_is_sparse(matcode))
580  strcpy(type1, MM_SPARSE_STR);
581  else
582  if (mm_is_dense(matcode))
583  strcpy(type1, MM_DENSE_STR);
584  else
585  return;
586 
587  /* check for element data type */
588  if (mm_is_real(matcode))
589  strcpy(type2, MM_REAL_STR);
590  else
591  if (mm_is_complex(matcode))
592  strcpy(type2, MM_COMPLEX_STR);
593  else
594  if (mm_is_pattern(matcode))
595  strcpy(type2, MM_PATTERN_STR);
596  else
597  if (mm_is_integer(matcode))
598  strcpy(type2, MM_INT_STR);
599  else
600  return;
601 
602  /* check for symmetry type */
603  if (mm_is_general(matcode))
604  strcpy(type3, MM_GENERAL_STR);
605  else
606  if (mm_is_symmetric(matcode))
607  strcpy(type3, MM_SYMM_STR);
608  else
609  if (mm_is_hermitian(matcode))
610  strcpy(type3, MM_HERM_STR);
611  else
612  if (mm_is_skew(matcode))
613  strcpy(type3, MM_SKEW_STR);
614  else
615  return;
616 
617  sprintf(buffer,"%s %s %s %s", type0, type1, type2, type3);
618  return;
619 }
620 } // namespace EpetraExt
#define MM_PATTERN_STR
#define mm_is_complex(typecode)
int mm_read_mtx_crd_data(FILE *f, int, int, int nz, int I[], int J[], double val[], MM_typecode matcode)
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz)
int mm_read_mtx_array_size(FILE *f, int *M, int *N)
#define mm_set_real(typecode)
int mm_is_valid(MM_typecode matcode)
#define MatrixMarketBanner
#define mm_set_hermitian(typecode)
#define MM_COULD_NOT_WRITE_FILE
#define mm_set_dense(typecode)
#define MM_HERM_STR
#define mm_set_matrix(typecode)
#define MM_SPARSE_STR
#define mm_is_sparse(typecode)
int mm_write_banner(FILE *f, MM_typecode matcode)
#define mm_set_integer(typecode)
#define MM_COMPLEX_STR
#define MM_INT_STR
#define MM_MAX_LINE_LENGTH
#define mm_is_pattern(typecode)
#define mm_is_dense(typecode)
#define mm_set_general(typecode)
#define MM_PREMATURE_EOF
#define mm_is_skew(typecode)
#define mm_is_symmetric(typecode)
void mm_typecode_to_str(MM_typecode matcode, char *buffer)
int mm_read_banner(FILE *f, MM_typecode *matcode)
int mm_write_mtx_crd_size(FILE *f, long long M, long long N, long long nz)
#define MM_UNSUPPORTED_TYPE
#define MM_MTX_STR
#define MM_NO_HEADER
int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[], double val[], MM_typecode matcode)
char MM_typecode[4]
#define mm_is_hermitian(typecode)
#define mm_is_real(typecode)
#define MM_SYMM_STR
#define mm_set_skew(typecode)
int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J, double **val, MM_typecode *matcode)
#define mm_set_complex(typecode)
#define MM_DENSE_STR
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J, double *real, double *imag, MM_typecode matcode)
#define MM_GENERAL_STR
#define mm_is_matrix(typecode)
#define mm_is_general(typecode)
#define mm_set_symmetric(typecode)
#define mm_is_integer(typecode)
int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_, double **val_, int **I_, int **J_)
#define mm_set_pattern(typecode)
#define mm_clear_typecode(typecode)
#define MM_REAL_STR
#define MM_COULD_NOT_READ_FILE
#define mm_set_sparse(typecode)
#define MM_SKEW_STR
int mm_write_mtx_array_size(FILE *f, long long M, long long N)
#define MM_MAX_TOKEN_LENGTH