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"; 8 9 /* 10 NOTES: 11 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. 16 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" 24 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; 40 41 PetscInitialize(&argc, &argv, (char *)0, help); 42 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 43 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 44 45 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Matrix Market example options", ""); 46 { 47 PetscCall(PetscOptionsString("-fin", "Input Matrix Market file", "", filein, filein, sizeof(filein), &flag)); 48 PetscCheck(flag, PETSC_COMM_SELF, PETSC_ERR_USER_INPUT, "Please use -fin <filename> to specify the input file name!"); 49 PetscCall(PetscOptionsString("-fout", "Output file in petsc sparse binary format", "", fileout, fileout, sizeof(fileout), &flag)); 50 PetscCheck(flag, PETSC_COMM_SELF, PETSC_ERR_USER_INPUT, "Please use -fout <filename> to specify the output file name!"); 51 PetscCall(PetscOptionsBool("-aij_only", "Use MATAIJ for all cases", "", aijonly, &aijonly, NULL)); 52 PetscCall(PetscOptionsFList("-permute", "Permute matrix and vector to solving in new ordering", "", MatOrderingList, ordering, ordering, sizeof(ordering), &permute)); 53 } 54 PetscOptionsEnd(); 55 56 /* Read in matrix */ 57 PetscCall(PetscFOpen(PETSC_COMM_SELF, filein, "r", &file)); 58 59 PetscCheck(mm_read_banner(file, &matcode) == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not process Matrix Market banner."); 60 61 /* This is how one can screen matrix types if their application */ 62 /* only supports a subset of the Matrix Market data types. */ 63 PetscCheck(mm_is_matrix(matcode) && mm_is_sparse(matcode), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input must be a sparse matrix. Market Market type: [%s]", mm_typecode_to_str(matcode)); 64 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; 69 70 /* Find out size of sparse matrix .... */ 71 PetscCheck(mm_read_mtx_crd_size(file, &M, &N, &nz) == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Size of sparse matrix is wrong."); 72 73 PetscCall(mm_write_banner(stdout, matcode)); 74 PetscCall(PetscPrintf(PETSC_COMM_SELF, "M: %d, N: %d, nnz: %d\n", M, N, nz)); 75 76 /* Reseve memory for matrices */ 77 PetscCall(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 */ 79 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) */ 83 84 for (i = 0; i < nz; i++) { 85 if (pattern) { 86 ninput = fscanf(file, "%d %d\n", &ia[i], &ja[i]); 87 PetscCheck(ninput >= 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file"); 88 val[i] = 1.0; 89 } else if (real) { 90 ninput = fscanf(file, "%d %d %lg\n", &ia[i], &ja[i], &val[i]); 91 PetscCheck(ninput >= 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file"); 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 PetscCall(PetscFClose(PETSC_COMM_SELF, file)); 103 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Reading matrix completes.\n")); 104 105 /* Create, preallocate, and then assemble the matrix */ 106 PetscCall(MatCreate(PETSC_COMM_SELF, &A)); 107 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N)); 108 109 if (symmetric && !aijonly) { 110 PetscCall(MatSetType(A, MATSEQSBAIJ)); 111 PetscCall(MatSetFromOptions(A)); 112 PetscCall(MatSetUp(A)); 113 PetscCall(MatSeqSBAIJSetPreallocation(A, 1, 0, rownz)); 114 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &sametype)); 115 PetscCheck(sametype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Only AIJ and SBAIJ are supported. Your mattype is not supported"); 116 } else { 117 PetscCall(MatSetType(A, MATSEQAIJ)); 118 PetscCall(MatSetFromOptions(A)); 119 PetscCall(MatSetUp(A)); 120 PetscCall(MatSeqAIJSetPreallocation(A, 0, rownz)); 121 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &sametype)); 122 PetscCheck(sametype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Only AIJ and SBAIJ are supported. Your mattype is not supported"); 123 } 124 125 /* Add zero to diagonals, in case the matrix missing diagonals */ 126 for (j = 0; j < M; j++) PetscCall(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++) PetscCall(MatSetValues(A, 1, &ia[j], 1, &ja[j], &val[j], INSERT_VALUES)); 129 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++) PetscCall(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 PetscCall(MatSetValues(A, 1, &ja[j], 1, &ia[j], &val[j], INSERT_VALUES)); 139 } 140 } 141 142 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 143 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 144 145 if (permute) { 146 Mat Aperm; 147 PetscCall(MatGetOrdering(A, ordering, &rowperm, &colperm)); 148 PetscCall(MatPermute(A, rowperm, colperm, &Aperm)); 149 PetscCall(MatDestroy(&A)); 150 A = Aperm; /* Replace original operator with permuted version */ 151 } 152 153 /* Write out matrix */ 154 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Writing matrix to binary file %s using PETSc %s format ...\n", fileout, (symmetric && !aijonly) ? "SBAIJ" : "AIJ")); 155 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, fileout, FILE_MODE_WRITE, &view)); 156 PetscCall(MatView(A, view)); 157 PetscCall(PetscViewerDestroy(&view)); 158 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Writing matrix completes.\n")); 159 160 PetscCall(PetscFree4(ia, ja, val, rownz)); 161 PetscCall(MatDestroy(&A)); 162 PetscCall(ISDestroy(&rowperm)); 163 PetscCall(ISDestroy(&colperm)); 164 PetscCall(PetscFinalize()); 165 return 0; 166 } 167 168 /*TEST 169 170 build: 171 requires: !complex double !defined(PETSC_USE_64BIT_INDICES) 172 depends: ex72mmio.c 173 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 178 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 183 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 188 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*/ 194