Actual source code: mmloader.c
1: #include "mmloader.h"
3: PetscErrorCode MatCreateFromMTX(Mat *A, const char *filein, PetscBool aijonly)
4: {
5: MM_typecode matcode;
6: FILE *file;
7: PetscInt M, N, ninput;
8: PetscInt *ia, *ja;
9: PetscInt i, j, nz, *rownz;
10: PetscScalar *val;
11: PetscBool sametype, symmetric = PETSC_FALSE, skew = PETSC_FALSE;
13: /* Read in matrix */
14: PetscFunctionBeginUser;
15: PetscCall(PetscFOpen(PETSC_COMM_SELF, filein, "r", &file));
16: PetscCheck(mm_read_banner(file, &matcode) == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not process Matrix Market banner.");
17: /* This is how one can screen matrix types if their application */
18: /* only supports a subset of the Matrix Market data types. */
19: PetscCheck(mm_is_matrix(matcode) && mm_is_sparse(matcode) && (mm_is_real(matcode) || mm_is_integer(matcode)), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input must be a sparse real or integer matrix. Market Market type: [%s]", mm_typecode_to_str(matcode));
21: if (mm_is_symmetric(matcode)) symmetric = PETSC_TRUE;
22: if (mm_is_skew(matcode)) skew = PETSC_TRUE;
24: /* Find out size of sparse matrix .... */
25: PetscCheck(mm_read_mtx_crd_size(file, &M, &N, &nz) == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Size of sparse matrix is wrong.");
27: /* Reserve memory for matrices */
28: PetscCall(PetscMalloc4(nz, &ia, nz, &ja, nz, &val, M, &rownz));
29: for (i = 0; i < M; i++) rownz[i] = 0;
31: /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
32: /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
33: /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
34: for (i = 0; i < nz; i++) {
35: ninput = fscanf(file, "%d %d %lg\n", &ia[i], &ja[i], &val[i]);
36: PetscCheck(ninput >= 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file");
37: ia[i]--;
38: ja[i]--; /* adjust from 1-based to 0-based */
39: if ((symmetric && aijonly) || skew) { /* transpose */
40: rownz[ia[i]]++;
41: if (ja[i] != ia[i]) rownz[ja[i]]++;
42: } else {
43: if (symmetric) rownz[ja[i]]++;
44: else rownz[ia[i]]++;
45: }
46: }
47: PetscCall(PetscFClose(PETSC_COMM_SELF, file));
49: /* Create, preallocate, and then assemble the matrix */
50: PetscCall(MatCreate(PETSC_COMM_SELF, A));
51: PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, M, N));
53: if (symmetric && !aijonly) {
54: PetscCall(MatSetType(*A, MATSEQSBAIJ));
55: PetscCall(MatSetFromOptions(*A));
56: PetscCall(MatSeqSBAIJSetPreallocation(*A, 1, 0, rownz));
57: PetscCall(PetscObjectTypeCompare((PetscObject)(*A), MATSEQSBAIJ, &sametype));
58: PetscCheck(sametype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Only AIJ and SBAIJ are supported. Your mattype is not supported");
59: } else {
60: PetscCall(MatSetType(*A, MATSEQAIJ));
61: PetscCall(MatSetFromOptions(*A));
62: PetscCall(MatSeqAIJSetPreallocation(*A, 0, rownz));
63: PetscCall(PetscObjectTypeCompare((PetscObject)(*A), MATSEQAIJ, &sametype));
64: PetscCheck(sametype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Only AIJ and SBAIJ are supported. Your mattype is not supported");
65: }
66: /* Add values to the matrix, these correspond to lower triangular part for symmetric or skew matrices */
67: if (!(symmetric && !aijonly))
68: for (j = 0; j < nz; j++) PetscCall(MatSetValues(*A, 1, &ia[j], 1, &ja[j], &val[j], INSERT_VALUES));
70: /* Add values to upper triangular part for some cases */
71: if (symmetric && !aijonly) {
72: /* MatrixMarket matrix stores symm matrix in lower triangular part. Take its transpose */
73: for (j = 0; j < nz; j++) PetscCall(MatSetValues(*A, 1, &ja[j], 1, &ia[j], &val[j], INSERT_VALUES));
74: }
75: if (skew) {
76: for (j = 0; j < nz; j++) {
77: val[j] = -val[j];
78: PetscCall(MatSetValues(*A, 1, &ja[j], 1, &ia[j], &val[j], INSERT_VALUES));
79: }
80: }
81: PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
82: PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
83: PetscCall(PetscFree4(ia, ja, val, rownz));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }