1 #include "mmloader.h" 2 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; 12 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)); 20 21 if (mm_is_symmetric(matcode)) symmetric = PETSC_TRUE; 22 if (mm_is_skew(matcode)) skew = PETSC_TRUE; 23 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."); 26 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; 30 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 rownz[ja[i]]++; 42 } else rownz[ia[i]]++; 43 } 44 PetscCall(PetscFClose(PETSC_COMM_SELF, file)); 45 46 /* Create, preallocate, and then assemble the matrix */ 47 PetscCall(MatCreate(PETSC_COMM_SELF, A)); 48 PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, M, N)); 49 50 if (symmetric && !aijonly) { 51 PetscCall(MatSetType(*A, MATSEQSBAIJ)); 52 PetscCall(MatSetFromOptions(*A)); 53 PetscCall(MatSetUp(*A)); 54 PetscCall(MatSeqSBAIJSetPreallocation(*A, 1, 0, rownz)); 55 PetscCall(PetscObjectTypeCompare((PetscObject)(*A), MATSEQSBAIJ, &sametype)); 56 PetscCheck(sametype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Only AIJ and SBAIJ are supported. Your mattype is not supported"); 57 } else { 58 PetscCall(MatSetType(*A, MATSEQAIJ)); 59 PetscCall(MatSetFromOptions(*A)); 60 PetscCall(MatSetUp(*A)); 61 PetscCall(MatSeqAIJSetPreallocation(*A, 0, rownz)); 62 PetscCall(PetscObjectTypeCompare((PetscObject)(*A), MATSEQAIJ, &sametype)); 63 PetscCheck(sametype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Only AIJ and SBAIJ are supported. Your mattype is not supported"); 64 } 65 /* Add values to the matrix, these correspond to lower triangular part for symmetric or skew matrices */ 66 for (j = 0; j < nz; j++) PetscCall(MatSetValues(*A, 1, &ia[j], 1, &ja[j], &val[j], INSERT_VALUES)); 67 68 /* Add values to upper triangular part for some cases */ 69 if (symmetric && aijonly) { 70 /* MatrixMarket matrix stores symm matrix in lower triangular part. Take its transpose */ 71 for (j = 0; j < nz; j++) PetscCall(MatSetValues(*A, 1, &ja[j], 1, &ia[j], &val[j], INSERT_VALUES)); 72 } 73 if (skew) { 74 for (j = 0; j < nz; j++) { 75 val[j] = -val[j]; 76 PetscCall(MatSetValues(*A, 1, &ja[j], 1, &ia[j], &val[j], INSERT_VALUES)); 77 } 78 } 79 PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY)); 80 PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY)); 81 PetscCall(PetscFree4(ia, ja, val, rownz)); 82 PetscFunctionReturn(PETSC_SUCCESS); 83 } 84