Actual source code: ex72.c
petsc-3.8.4 2018-03-24
2: #include <petscmat.h>
4: static char help[] = "Read in a Symmetric matrix in MatrixMarket format (only the lower triangle). \n\
5: Assemble it to a PETSc sparse SBAIJ (upper triangle) matrix. \n\
6: Write it in a AIJ matrix (entire matrix) to a file. \n\
7: Input parameters are: \n\
8: -fin <filename> : input file \n\
9: -fout <filename> : output file \n\n";
11: int main(int argc,char **args)
12: {
13: Mat A;
14: char filein[PETSC_MAX_PATH_LEN],fileout[PETSC_MAX_PATH_LEN],buf[PETSC_MAX_PATH_LEN];
15: PetscInt i,m,n,nnz;
17: PetscMPIInt size;
18: PetscScalar *val,zero=0.0;
19: FILE *file;
20: PetscViewer view;
21: int *row,*col,*rownz;
22: PetscBool flg;
24: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
25: #if defined(PETSC_USE_COMPLEX)
26: SETERRQ(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
27: PetscFinalize();
28: return 0;
29: #endif
31: MPI_Comm_size(PETSC_COMM_WORLD,&size);
32: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Uniprocessor Example only\n");
34: /* Read in matrix and RHS */
35: PetscOptionsGetString(NULL,NULL,"-fin",filein,PETSC_MAX_PATH_LEN,&flg);
36: if (!flg) SETERRQ(PETSC_COMM_SELF,1,"Must indicate input file with -fin option");
37: PetscFOpen(PETSC_COMM_SELF,filein,"r",&file);
39: /* process header with comments */
40: do {
41: char *str = fgets(buf,PETSC_MAX_PATH_LEN-1,file);
42: if (!str) SETERRQ(PETSC_COMM_SELF,1,"Incorrect format in file");
43: }while (buf[0] == '%');
45: /* The first non-comment line has the matrix dimensions */
46: sscanf(buf,"%d %d %d\n",&m,&n,&nnz);
47: PetscPrintf (PETSC_COMM_SELF,"m = %d, n = %d, nnz = %d\n",m,n,nnz);
49: /* reseve memory for matrices */
50: PetscMalloc4(nnz,&row,nnz,&col,nnz,&val,m,&rownz);
51: for (i=0; i<m; i++) rownz[i] = 1; /* add 0.0 to diagonal entries */
53: for (i=0; i<nnz; i++) {
54: fscanf(file,"%d %d %le\n",&row[i],&col[i],(double*)&val[i]);
55: if (ierr == EOF) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"i=%d, reach EOF\n",i);
56: row[i]--; col[i]--; /* adjust from 1-based to 0-based */
57: rownz[col[i]]++;
58: }
59: fclose(file);
60: PetscPrintf(PETSC_COMM_SELF,"Read file completes.\n");
62: /* Creat and asseble SBAIJ matrix */
63: MatCreate(PETSC_COMM_SELF,&A);
64: MatSetType(A,MATSBAIJ);
65: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);
66: MatSetFromOptions(A);
67: MatSeqSBAIJSetPreallocation(A,1,0,rownz);
69: /* Add zero to diagonals, in case the matrix missing diagonals */
70: for (i=0; i<m; i++){
71: MatSetValues(A,1,&i,1,&i,&zero,INSERT_VALUES);
72: }
73: for (i=0; i<nnz; i++) {
74: MatSetValues(A,1,&col[i],1,&row[i],&val[i],INSERT_VALUES);
75: }
76: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
77: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
78: PetscPrintf(PETSC_COMM_SELF,"Assemble SBAIJ matrix completes.\n");
80: /* Write the entire matrix in AIJ format to a file */
81: PetscOptionsGetString(NULL,NULL,"-fout",fileout,PETSC_MAX_PATH_LEN,&flg);
82: if (flg) {
83: PetscPrintf(PETSC_COMM_SELF,"Write the entire matrix in AIJ format to file %s\n",fileout);
84: PetscViewerBinaryOpen(PETSC_COMM_WORLD,fileout,FILE_MODE_WRITE,&view);
85: MatView(A,view);
86: PetscViewerDestroy(&view);
87: }
89: PetscFree4(row,col,val,rownz);
90: MatDestroy(&A);
91: PetscFinalize();
92: return ierr;
93: }