xref: /petsc/src/mat/tests/ex72.c (revision 3c859ba3a04a72e8efcb87bc7ffc046a6cbab413)
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   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr);
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","");CHKERRQ(ierr);
45   {
46     ierr = PetscOptionsString("-fin","Input Matrix Market file","",filein,filein,sizeof(filein),&flag);CHKERRQ(ierr);
47     PetscCheck(flag,PETSC_COMM_SELF,PETSC_ERR_USER_INPUT,"Please use -fin <filename> to specify the input file name!");
48     ierr = PetscOptionsString("-fout","Output file in petsc sparse binary format","",fileout,fileout,sizeof(fileout),&flag);CHKERRQ(ierr);
49     PetscCheck(flag,PETSC_COMM_SELF,PETSC_ERR_USER_INPUT,"Please use -fout <filename> to specify the output file name!");
50     ierr = PetscOptionsBool("-aij_only","Use MATAIJ for all cases","",aijonly,&aijonly,NULL);CHKERRQ(ierr);
51     ierr = PetscOptionsFList("-permute","Permute matrix and vector to solving in new ordering","",MatOrderingList,ordering,ordering,sizeof(ordering),&permute);CHKERRQ(ierr);
52   }
53   ierr = PetscOptionsEnd();CHKERRQ(ierr);
54 
55   /* Read in matrix */
56   ierr = PetscFOpen(PETSC_COMM_SELF,filein,"r",&file);CHKERRQ(ierr);
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   ierr = mm_write_banner(stdout, matcode);CHKERRQ(ierr);
73   ierr = PetscPrintf(PETSC_COMM_SELF,"M: %d, N: %d, nnz: %d\n",M,N,nz);CHKERRQ(ierr);
74 
75   /* Reseve memory for matrices */
76   ierr = PetscMalloc4(nz,&ia,nz,&ja,nz,&val,M,&rownz);CHKERRQ(ierr);
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       PetscCheckFalse(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       PetscCheckFalse(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   ierr = PetscFClose(PETSC_COMM_SELF,file);CHKERRQ(ierr);
101   ierr = PetscPrintf(PETSC_COMM_SELF,"Reading matrix completes.\n");CHKERRQ(ierr);
102 
103   /* Create, preallocate, and then assemble the matrix */
104   ierr = MatCreate(PETSC_COMM_SELF,&A);CHKERRQ(ierr);
105   ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr);
106 
107   if (symmetric && !aijonly) {
108     ierr = MatSetType(A,MATSEQSBAIJ);CHKERRQ(ierr);
109     ierr = MatSetFromOptions(A);CHKERRQ(ierr);
110     ierr = MatSetUp(A);CHKERRQ(ierr);
111     ierr = MatSeqSBAIJSetPreallocation(A,1,0,rownz);CHKERRQ(ierr);
112     ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQSBAIJ,&sametype);CHKERRQ(ierr);
113     PetscCheckFalse(!sametype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Only AIJ and SBAIJ are supported. Your mattype is not supported");
114   } else {
115     ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr);
116     ierr = MatSetFromOptions(A);CHKERRQ(ierr);
117     ierr = MatSetUp(A);CHKERRQ(ierr);
118     ierr = MatSeqAIJSetPreallocation(A,0,rownz);CHKERRQ(ierr);
119     ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&sametype);CHKERRQ(ierr);
120     PetscCheckFalse(!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++)  {
125     ierr = MatSetValues(A,1,&j,1,&j,&zero,INSERT_VALUES);CHKERRQ(ierr);
126   }
127   /* Add values to the matrix, these correspond to lower triangular part for symmetric or skew matrices */
128   for (j=0; j<nz; j++) {
129     ierr = MatSetValues(A,1,&ia[j],1,&ja[j],&val[j],INSERT_VALUES);CHKERRQ(ierr);
130   }
131 
132   /* Add values to upper triangular part for some cases */
133   if (symmetric && aijonly) {
134     /* MatrixMarket matrix stores symm matrix in lower triangular part. Take its transpose */
135     for (j=0; j<nz; j++) {
136       ierr = MatSetValues(A,1,&ja[j],1,&ia[j],&val[j],INSERT_VALUES);CHKERRQ(ierr);
137     }
138   }
139   if (skew) {
140     for (j=0; j<nz; j++) {
141       val[j] = -val[j];
142       ierr = MatSetValues(A,1,&ja[j],1,&ia[j],&val[j],INSERT_VALUES);CHKERRQ(ierr);
143     }
144   }
145 
146   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
147   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
148 
149   if (permute) {
150     Mat Aperm;
151     ierr = MatGetOrdering(A,ordering,&rowperm,&colperm);CHKERRQ(ierr);
152     ierr = MatPermute(A,rowperm,colperm,&Aperm);CHKERRQ(ierr);
153     ierr = MatDestroy(&A);CHKERRQ(ierr);
154     A    = Aperm;               /* Replace original operator with permuted version */
155   }
156 
157   /* Write out matrix */
158   ierr = PetscPrintf(PETSC_COMM_SELF,"Writing matrix to binary file %s using PETSc %s format ...\n",fileout,(symmetric && !aijonly)?"SBAIJ":"AIJ");CHKERRQ(ierr);
159   ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,fileout,FILE_MODE_WRITE,&view);CHKERRQ(ierr);
160   ierr = MatView(A,view);CHKERRQ(ierr);
161   ierr = PetscViewerDestroy(&view);CHKERRQ(ierr);
162   ierr = PetscPrintf(PETSC_COMM_SELF,"Writing matrix completes.\n");CHKERRQ(ierr);
163 
164   ierr = PetscFree4(ia,ja,val,rownz);CHKERRQ(ierr);
165   ierr = MatDestroy(&A);CHKERRQ(ierr);
166   ierr = ISDestroy(&rowperm);CHKERRQ(ierr);
167   ierr = ISDestroy(&colperm);CHKERRQ(ierr);
168   ierr = PetscFinalize();CHKERRQ(ierr);
169   return 0;
170 }
171 
172 /*TEST
173 
174    build:
175       requires:  !complex double !defined(PETSC_USE_64BIT_INDICES)
176       depends: ex72mmio.c
177 
178    test:
179       suffix: 1
180       args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/amesos2_test_mat0.mtx -fout petscmat.aij
181       output_file: output/ex72_1.out
182 
183    test:
184       suffix: 2
185       args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/LFAT5.mtx -fout petscmat.sbaij
186       output_file: output/ex72_2.out
187 
188    test:
189       suffix: 3
190       args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/m_05_05_crk.mtx -fout petscmat2.aij
191       output_file: output/ex72_3.out
192 
193    test:
194       suffix: 4
195       args: -fin ${wPETSC_DIR}/share/petsc/datafiles/matrices/amesos2_test_mat0.mtx -fout petscmat.aij -permute rcm
196       output_file: output/ex72_4.out
197 TEST*/
198