EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
64using namespace EpetraExt;
65namespace EpetraExt {
66int 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
150int 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
235int 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
241int 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
271int 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
301int 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
330int 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
344int 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
377int 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
403int 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
437int 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)
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
496int 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
510int 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
556void 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_is_dense(typecode)
#define mm_is_complex(typecode)
#define MM_MTX_STR
#define mm_set_complex(typecode)
#define mm_set_hermitian(typecode)
#define MM_COULD_NOT_READ_FILE
#define MM_PATTERN_STR
#define mm_set_general(typecode)
char MM_typecode[4]
#define MM_MAX_LINE_LENGTH
#define mm_set_skew(typecode)
#define MM_SKEW_STR
#define MM_PREMATURE_EOF
#define mm_set_sparse(typecode)
#define mm_is_matrix(typecode)
#define mm_is_hermitian(typecode)
#define mm_clear_typecode(typecode)
#define MM_GENERAL_STR
#define MM_MAX_TOKEN_LENGTH
#define MM_COMPLEX_STR
#define mm_set_symmetric(typecode)
#define mm_is_real(typecode)
#define mm_is_skew(typecode)
#define MM_SPARSE_STR
#define MM_UNSUPPORTED_TYPE
#define MM_COULD_NOT_WRITE_FILE
#define mm_set_integer(typecode)
#define MM_NO_HEADER
#define MM_HERM_STR
#define mm_is_pattern(typecode)
#define mm_is_general(typecode)
#define mm_set_matrix(typecode)
#define mm_set_dense(typecode)
#define MM_SYMM_STR
#define mm_set_pattern(typecode)
#define mm_is_integer(typecode)
#define MM_DENSE_STR
#define MM_REAL_STR
#define mm_is_symmetric(typecode)
#define MM_INT_STR
#define mm_is_sparse(typecode)
#define mm_set_real(typecode)
#define MatrixMarketBanner
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_, double **val_, int **I_, int **J_)
void mm_typecode_to_str(MM_typecode matcode, char *buffer)
int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J, double **val, MM_typecode *matcode)
int mm_write_mtx_array_size(FILE *f, long long M, long long N)
int mm_write_banner(FILE *f, MM_typecode matcode)
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz)
int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[], double val[], MM_typecode matcode)
int mm_write_mtx_crd_size(FILE *f, long long M, long long N, long long nz)
int mm_read_banner(FILE *f, MM_typecode *matcode)
int mm_is_valid(MM_typecode matcode)
const int N
int mm_read_mtx_array_size(FILE *f, int *M, int *N)
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_entry(FILE *f, int *I, int *J, double *real, double *imag, MM_typecode matcode)