Actual source code: ex72.c

  1: static char help[] = "Read a non-complex sparse matrix from a Matrix Market (v. 2.0) file\n\
  2: and write it to a file in petsc sparse binary format. If the matrix is symmetric, the binary file is in \n\
  3: PETSc MATSBAIJ format, otherwise it is in MATAIJ format \n\
  4: Usage:  ./ex72 -fin <infile> -fout <outfile> \n\
  5: (See https://math.nist.gov/MatrixMarket/ for details.)\n\
  6: The option -permute <natural,rcm,nd,...> permutes the matrix using the ordering type.\n\
  7: The option -aij_only allows to use MATAIJ for all cases.\n\\n";

  9: /*
 10:    NOTES:

 12:    1) Matrix Market files are always 1-based, i.e. the index of the first
 13:       element of a matrix is (1,1), not (0,0) as in C.  ADJUST THESE
 14:       OFFSETS ACCORDINGLY offsets accordingly when reading and writing
 15:       to files.

 17:    2) ANSI C requires one to use the "l" format modifier when reading
 18:       double precision floating point numbers in scanf() and
 19:       its variants.  For example, use "%lf", "%lg", or "%le"
 20:       when reading doubles, otherwise errors will occur.
 21: */
 22: #include <petscmat.h>
 23: #include "ex72mmio.h"

 25: int main(int argc, char **argv)
 26: {
 27:   MM_typecode  matcode;
 28:   FILE        *file;
 29:   PetscInt     M, N, ninput;
 30:   PetscInt    *ia, *ja;
 31:   Mat          A;
 32:   char         filein[PETSC_MAX_PATH_LEN], fileout[PETSC_MAX_PATH_LEN];
 33:   char         ordering[256] = MATORDERINGRCM;
 34:   PetscInt     i, j, nz, *rownz;
 35:   PetscScalar *val, zero = 0.0;
 36:   PetscViewer  view;
 37:   PetscBool    sametype, flag, symmetric = PETSC_FALSE, skew = PETSC_FALSE, real = PETSC_FALSE, pattern = PETSC_FALSE, aijonly = PETSC_FALSE, permute = PETSC_FALSE;
 38:   IS           rowperm = NULL, colperm = NULL;
 39:   PetscMPIInt  size;

 41:   PetscInitialize(&argc, &argv, (char *)0, help);
 42:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 45:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Matrix Market example options", "");
 46:   {
 47:     PetscOptionsString("-fin", "Input Matrix Market file", "", filein, filein, sizeof(filein), &flag);
 49:     PetscOptionsString("-fout", "Output file in petsc sparse binary format", "", fileout, fileout, sizeof(fileout), &flag);
 51:     PetscOptionsBool("-aij_only", "Use MATAIJ for all cases", "", aijonly, &aijonly, NULL);
 52:     PetscOptionsFList("-permute", "Permute matrix and vector to solving in new ordering", "", MatOrderingList, ordering, ordering, sizeof(ordering), &permute);
 53:   }
 54:   PetscOptionsEnd();

 56:   /* Read in matrix */
 57:   PetscFOpen(PETSC_COMM_SELF, filein, "r", &file);


 61:   /*  This is how one can screen matrix types if their application */
 62:   /*  only supports a subset of the Matrix Market data types.      */

 65:   if (mm_is_symmetric(matcode)) symmetric = PETSC_TRUE;
 66:   if (mm_is_skew(matcode)) skew = PETSC_TRUE;
 67:   if (mm_is_real(matcode)) real = PETSC_TRUE;
 68:   if (mm_is_pattern(matcode)) pattern = PETSC_TRUE;

 70:   /* Find out size of sparse matrix .... */

 73:   mm_write_banner(stdout, matcode);
 74:   PetscPrintf(PETSC_COMM_SELF, "M: %d, N: %d, nnz: %d\n", M, N, nz);

 76:   /* Reseve memory for matrices */
 77:   PetscMalloc4(nz, &ia, nz, &ja, nz, &val, M, &rownz);
 78:   for (i = 0; i < M; i++) rownz[i] = 1; /* Since we will add 0.0 to diagonal entries */

 80:   /* NOTE: when reading in doubles, ANSI C requires the use of the "l"  */
 81:   /*   specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
 82:   /*  (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15)            */

 84:   for (i = 0; i < nz; i++) {
 85:     if (pattern) {
 86:       ninput = fscanf(file, "%d %d\n", &ia[i], &ja[i]);
 88:       val[i] = 1.0;
 89:     } else if (real) {
 90:       ninput = fscanf(file, "%d %d %lg\n", &ia[i], &ja[i], &val[i]);
 92:     }
 93:     ia[i]--;
 94:     ja[i]--;                                /* adjust from 1-based to 0-based */
 95:     if (ia[i] != ja[i]) {                   /* already counted the diagonals above */
 96:       if ((symmetric && aijonly) || skew) { /* transpose */
 97:         rownz[ia[i]]++;
 98:         rownz[ja[i]]++;
 99:       } else rownz[ia[i]]++;
100:     }
101:   }
102:   PetscFClose(PETSC_COMM_SELF, file);
103:   PetscPrintf(PETSC_COMM_SELF, "Reading matrix completes.\n");

105:   /* Create, preallocate, and then assemble the matrix */
106:   MatCreate(PETSC_COMM_SELF, &A);
107:   MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N);

109:   if (symmetric && !aijonly) {
110:     MatSetType(A, MATSEQSBAIJ);
111:     MatSetFromOptions(A);
112:     MatSetUp(A);
113:     MatSeqSBAIJSetPreallocation(A, 1, 0, rownz);
114:     PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &sametype);
116:   } else {
117:     MatSetType(A, MATSEQAIJ);
118:     MatSetFromOptions(A);
119:     MatSetUp(A);
120:     MatSeqAIJSetPreallocation(A, 0, rownz);
121:     PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &sametype);
123:   }

125:   /* Add zero to diagonals, in case the matrix missing diagonals */
126:   for (j = 0; j < M; j++) MatSetValues(A, 1, &j, 1, &j, &zero, INSERT_VALUES);
127:   /* Add values to the matrix, these correspond to lower triangular part for symmetric or skew matrices */
128:   for (j = 0; j < nz; j++) MatSetValues(A, 1, &ia[j], 1, &ja[j], &val[j], INSERT_VALUES);

130:   /* Add values to upper triangular part for some cases */
131:   if (symmetric && aijonly) {
132:     /* MatrixMarket matrix stores symm matrix in lower triangular part. Take its transpose */
133:     for (j = 0; j < nz; j++) MatSetValues(A, 1, &ja[j], 1, &ia[j], &val[j], INSERT_VALUES);
134:   }
135:   if (skew) {
136:     for (j = 0; j < nz; j++) {
137:       val[j] = -val[j];
138:       MatSetValues(A, 1, &ja[j], 1, &ia[j], &val[j], INSERT_VALUES);
139:     }
140:   }

142:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
143:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

145:   if (permute) {
146:     Mat Aperm;
147:     MatGetOrdering(A, ordering, &rowperm, &colperm);
148:     MatPermute(A, rowperm, colperm, &Aperm);
149:     MatDestroy(&A);
150:     A = Aperm; /* Replace original operator with permuted version */
151:   }

153:   /* Write out matrix */
154:   PetscPrintf(PETSC_COMM_SELF, "Writing matrix to binary file %s using PETSc %s format ...\n", fileout, (symmetric && !aijonly) ? "SBAIJ" : "AIJ");
155:   PetscViewerBinaryOpen(PETSC_COMM_SELF, fileout, FILE_MODE_WRITE, &view);
156:   MatView(A, view);
157:   PetscViewerDestroy(&view);
158:   PetscPrintf(PETSC_COMM_SELF, "Writing matrix completes.\n");

160:   PetscFree4(ia, ja, val, rownz);
161:   MatDestroy(&A);
162:   ISDestroy(&rowperm);
163:   ISDestroy(&colperm);
164:   PetscFinalize();
165:   return 0;
166: }

168: /*TEST

170:    build:
171:       requires:  !complex double !defined(PETSC_USE_64BIT_INDICES)
172:       depends: ex72mmio.c

174:    test:
175:       suffix: 1
176:       args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/amesos2_test_mat0.mtx -fout petscmat.aij
177:       output_file: output/ex72_1.out

179:    test:
180:       suffix: 2
181:       args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/LFAT5.mtx -fout petscmat.sbaij
182:       output_file: output/ex72_2.out

184:    test:
185:       suffix: 3
186:       args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/m_05_05_crk.mtx -fout petscmat2.aij
187:       output_file: output/ex72_3.out

189:    test:
190:       suffix: 4
191:       args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/amesos2_test_mat0.mtx -fout petscmat.aij -permute rcm
192:       output_file: output/ex72_4.out
193: TEST*/