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,ierr,size,*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 40 PetscInitialize(&argc,&argv,(char *)0,help); 41 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 42 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 43 44 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Matrix Market example options","");PetscCall(ierr); 45 { 46 PetscCall(PetscOptionsString("-fin","Input Matrix Market file","",filein,filein,sizeof(filein),&flag)); 47 PetscCheck(flag,PETSC_COMM_SELF,PETSC_ERR_USER_INPUT,"Please use -fin <filename> to specify the input file name!"); 48 PetscCall(PetscOptionsString("-fout","Output file in petsc sparse binary format","",fileout,fileout,sizeof(fileout),&flag)); 49 PetscCheck(flag,PETSC_COMM_SELF,PETSC_ERR_USER_INPUT,"Please use -fout <filename> to specify the output file name!"); 50 PetscCall(PetscOptionsBool("-aij_only","Use MATAIJ for all cases","",aijonly,&aijonly,NULL)); 51 PetscCall(PetscOptionsFList("-permute","Permute matrix and vector to solving in new ordering","",MatOrderingList,ordering,ordering,sizeof(ordering),&permute)); 52 } 53 ierr = PetscOptionsEnd();PetscCall(ierr); 54 55 /* Read in matrix */ 56 PetscCall(PetscFOpen(PETSC_COMM_SELF,filein,"r",&file)); 57 58 PetscCheck(mm_read_banner(file, &matcode) == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Could not process Matrix Market banner."); 59 60 /* This is how one can screen matrix types if their application */ 61 /* only supports a subset of the Matrix Market data types. */ 62 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)); 63 64 if (mm_is_symmetric(matcode)) symmetric = PETSC_TRUE; 65 if (mm_is_skew(matcode)) skew = PETSC_TRUE; 66 if (mm_is_real(matcode)) real = PETSC_TRUE; 67 if (mm_is_pattern(matcode)) pattern = PETSC_TRUE; 68 69 /* Find out size of sparse matrix .... */ 70 PetscCheck(mm_read_mtx_crd_size(file, &M, &N, &nz) == 0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Size of sparse matrix is wrong."); 71 72 PetscCall(mm_write_banner(stdout, matcode)); 73 PetscCall(PetscPrintf(PETSC_COMM_SELF,"M: %d, N: %d, nnz: %d\n",M,N,nz)); 74 75 /* Reseve memory for matrices */ 76 PetscCall(PetscMalloc4(nz,&ia,nz,&ja,nz,&val,M,&rownz)); 77 for (i=0; i<M; i++) rownz[i] = 1; /* Since we will add 0.0 to diagonal entries */ 78 79 /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */ 80 /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */ 81 /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */ 82 83 for (i=0; i<nz; i++) { 84 if (pattern) { 85 ninput = fscanf(file, "%d %d\n", &ia[i], &ja[i]); 86 PetscCheck(ninput >= 2,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file"); 87 val[i] = 1.0; 88 } else if (real) { 89 ninput = fscanf(file, "%d %d %lg\n", &ia[i], &ja[i], &val[i]); 90 PetscCheck(ninput >= 3,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file"); 91 } 92 ia[i]--; ja[i]--; /* adjust from 1-based to 0-based */ 93 if (ia[i] != ja[i]) { /* already counted the diagonals above */ 94 if ((symmetric && aijonly) || skew) { /* transpose */ 95 rownz[ia[i]]++; 96 rownz[ja[i]]++; 97 } else rownz[ia[i]]++; 98 } 99 } 100 PetscCall(PetscFClose(PETSC_COMM_SELF,file)); 101 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Reading matrix completes.\n")); 102 103 /* Create, preallocate, and then assemble the matrix */ 104 PetscCall(MatCreate(PETSC_COMM_SELF,&A)); 105 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N)); 106 107 if (symmetric && !aijonly) { 108 PetscCall(MatSetType(A,MATSEQSBAIJ)); 109 PetscCall(MatSetFromOptions(A)); 110 PetscCall(MatSetUp(A)); 111 PetscCall(MatSeqSBAIJSetPreallocation(A,1,0,rownz)); 112 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&sametype)); 113 PetscCheck(sametype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Only AIJ and SBAIJ are supported. Your mattype is not supported"); 114 } else { 115 PetscCall(MatSetType(A,MATSEQAIJ)); 116 PetscCall(MatSetFromOptions(A)); 117 PetscCall(MatSetUp(A)); 118 PetscCall(MatSeqAIJSetPreallocation(A,0,rownz)); 119 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&sametype)); 120 PetscCheck(sametype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Only AIJ and SBAIJ are supported. Your mattype is not supported"); 121 } 122 123 /* Add zero to diagonals, in case the matrix missing diagonals */ 124 for (j=0; j<M; j++) PetscCall(MatSetValues(A,1,&j,1,&j,&zero,INSERT_VALUES)); 125 /* Add values to the matrix, these correspond to lower triangular part for symmetric or skew matrices */ 126 for (j=0; j<nz; j++) PetscCall(MatSetValues(A,1,&ia[j],1,&ja[j],&val[j],INSERT_VALUES)); 127 128 /* Add values to upper triangular part for some cases */ 129 if (symmetric && aijonly) { 130 /* MatrixMarket matrix stores symm matrix in lower triangular part. Take its transpose */ 131 for (j=0; j<nz; j++) PetscCall(MatSetValues(A,1,&ja[j],1,&ia[j],&val[j],INSERT_VALUES)); 132 } 133 if (skew) { 134 for (j=0; j<nz; j++) { 135 val[j] = -val[j]; 136 PetscCall(MatSetValues(A,1,&ja[j],1,&ia[j],&val[j],INSERT_VALUES)); 137 } 138 } 139 140 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 141 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 142 143 if (permute) { 144 Mat Aperm; 145 PetscCall(MatGetOrdering(A,ordering,&rowperm,&colperm)); 146 PetscCall(MatPermute(A,rowperm,colperm,&Aperm)); 147 PetscCall(MatDestroy(&A)); 148 A = Aperm; /* Replace original operator with permuted version */ 149 } 150 151 /* Write out matrix */ 152 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Writing matrix to binary file %s using PETSc %s format ...\n",fileout,(symmetric && !aijonly)?"SBAIJ":"AIJ")); 153 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF,fileout,FILE_MODE_WRITE,&view)); 154 PetscCall(MatView(A,view)); 155 PetscCall(PetscViewerDestroy(&view)); 156 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Writing matrix completes.\n")); 157 158 PetscCall(PetscFree4(ia,ja,val,rownz)); 159 PetscCall(MatDestroy(&A)); 160 PetscCall(ISDestroy(&rowperm)); 161 PetscCall(ISDestroy(&colperm)); 162 PetscCall(PetscFinalize()); 163 return 0; 164 } 165 166 /*TEST 167 168 build: 169 requires: !complex double !defined(PETSC_USE_64BIT_INDICES) 170 depends: ex72mmio.c 171 172 test: 173 suffix: 1 174 args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/amesos2_test_mat0.mtx -fout petscmat.aij 175 output_file: output/ex72_1.out 176 177 test: 178 suffix: 2 179 args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/LFAT5.mtx -fout petscmat.sbaij 180 output_file: output/ex72_2.out 181 182 test: 183 suffix: 3 184 args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/m_05_05_crk.mtx -fout petscmat2.aij 185 output_file: output/ex72_3.out 186 187 test: 188 suffix: 4 189 args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/amesos2_test_mat0.mtx -fout petscmat.aij -permute rcm 190 output_file: output/ex72_4.out 191 TEST*/ 192