xref: /petsc/src/mat/tutorials/ex16.c (revision 84ff8c8b54fd7c9cb88641c01dfe6357ec5f72d0)
1 
2 static char help[] = "Reads a matrix from PETSc binary file. Use for view or investigating matrix data structure. \n\n";
3 /*
4  Example:
5       ./ex16 -f <matrix file> -a_mat_view draw -draw_pause -1
6       ./ex16 -f <matrix file> -a_mat_view ascii::ascii_info
7  */
8 
9 #include <petscmat.h>
10 int main(int argc,char **args)
11 {
12   Mat               A,Asp;
13   PetscViewer       fd;                        /* viewer */
14   char              file[PETSC_MAX_PATH_LEN];  /* input file name */
15   PetscErrorCode    ierr;
16   PetscInt          m,n,rstart,rend;
17   PetscBool         flg;
18   PetscInt          row,ncols,j,nrows,nnzA=0,nnzAsp=0;
19   const PetscInt    *cols;
20   const PetscScalar *vals;
21   PetscReal         norm,percent,val,dtol=1.e-16;
22   PetscMPIInt       rank;
23   MatInfo           matinfo;
24   PetscInt          Dnnz,Onnz;
25 
26   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
27   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr);
28 
29   /* Determine files from which we read the linear systems. */
30   ierr = PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);CHKERRQ(ierr);
31   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option");
32 
33   /* Open binary file.  Note that we use FILE_MODE_READ to indicate
34      reading from this file. */
35   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);
36 
37   /* Load the matrix; then destroy the viewer. */
38   ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
39   ierr = MatSetOptionsPrefix(A,"a_");CHKERRQ(ierr);
40   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
41   ierr = MatLoad(A,fd);CHKERRQ(ierr);
42   ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);
43   ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr);
44   ierr = MatGetInfo(A,MAT_LOCAL,&matinfo);CHKERRQ(ierr);
45   /*printf("matinfo.nz_used %g\n",matinfo.nz_used);*/
46 
47   /* Get a sparse matrix Asp by dumping zero entries of A */
48   ierr = MatCreate(PETSC_COMM_WORLD,&Asp);CHKERRQ(ierr);
49   ierr = MatSetSizes(Asp,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
50   ierr = MatSetOptionsPrefix(Asp,"asp_");CHKERRQ(ierr);
51   ierr = MatSetFromOptions(Asp);CHKERRQ(ierr);
52   Dnnz = (PetscInt)matinfo.nz_used/m + 1;
53   Onnz = Dnnz/2;
54   printf("Dnnz %d %d\n",Dnnz,Onnz);
55   ierr = MatSeqAIJSetPreallocation(Asp,Dnnz,NULL);CHKERRQ(ierr);
56   ierr = MatMPIAIJSetPreallocation(Asp,Dnnz,NULL,Onnz,NULL);CHKERRQ(ierr);
57   /* The allocation above is approximate so we must set this option to be permissive.
58    * Real code should preallocate exactly. */
59   ierr = MatSetOption(Asp,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
60 
61   /* Check zero rows */
62   ierr  = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
63   nrows = 0;
64   for (row=rstart; row<rend; row++) {
65     ierr  = MatGetRow(A,row,&ncols,&cols,&vals);CHKERRQ(ierr);
66     nnzA += ncols;
67     norm  = 0.0;
68     for (j=0; j<ncols; j++) {
69       val = PetscAbsScalar(vals[j]);
70       if (norm < val) norm = norm;
71       if (val > dtol) {
72         ierr = MatSetValues(Asp,1,&row,1,&cols[j],&vals[j],INSERT_VALUES);CHKERRQ(ierr);
73         nnzAsp++;
74       }
75     }
76     if (!norm) nrows++;
77     ierr = MatRestoreRow(A,row,&ncols,&cols,&vals);CHKERRQ(ierr);
78   }
79   ierr = MatAssemblyBegin(Asp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
80   ierr = MatAssemblyEnd(Asp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
81 
82   percent=(PetscReal)nnzA*100/(m*n);
83   ierr   = PetscPrintf(PETSC_COMM_SELF," [%d] Matrix A local size %d,%d; nnzA %d, %g percent; No. of zero rows: %d\n",rank,m,n,nnzA,percent,nrows);CHKERRQ(ierr);
84   percent=(PetscReal)nnzAsp*100/(m*n);
85   ierr   = PetscPrintf(PETSC_COMM_SELF," [%d] Matrix Asp nnzAsp %d, %g percent\n",rank,nnzAsp,percent);CHKERRQ(ierr);
86 
87   /* investigate matcoloring for Asp */
88   PetscBool Asp_coloring = PETSC_FALSE;
89   ierr = PetscOptionsHasName(NULL,NULL,"-Asp_color",&Asp_coloring);CHKERRQ(ierr);
90   if (Asp_coloring) {
91     MatColoring   mc;
92     ISColoring    iscoloring;
93     MatFDColoring matfdcoloring;
94     ierr = PetscPrintf(PETSC_COMM_WORLD," Create coloring of Asp...\n");CHKERRQ(ierr);
95     ierr = MatColoringCreate(Asp,&mc);CHKERRQ(ierr);
96     ierr = MatColoringSetType(mc,MATCOLORINGSL);CHKERRQ(ierr);
97     ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr);
98     ierr = MatColoringApply(mc,&iscoloring);CHKERRQ(ierr);
99     ierr = MatColoringDestroy(&mc);CHKERRQ(ierr);
100     ierr = MatFDColoringCreate(Asp,iscoloring,&matfdcoloring);CHKERRQ(ierr);
101     ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr);
102     ierr = MatFDColoringSetUp(Asp,iscoloring,matfdcoloring);CHKERRQ(ierr);
103     /*ierr = MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);*/
104     ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr);
105     ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr);
106   }
107 
108   /* Write Asp in binary for study - see ~petsc/src/mat/tests/ex124.c */
109   PetscBool Asp_write = PETSC_FALSE;
110   ierr = PetscOptionsHasName(NULL,NULL,"-Asp_write",&Asp_write);CHKERRQ(ierr);
111   if (Asp_write) {
112     PetscViewer viewer;
113     ierr = PetscPrintf(PETSC_COMM_SELF,"Write Asp into file Asp.dat ...\n");CHKERRQ(ierr);
114     ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Asp.dat",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
115     ierr = MatView(Asp,viewer);CHKERRQ(ierr);
116     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
117   }
118 
119   ierr = MatDestroy(&A);CHKERRQ(ierr);
120   ierr = MatDestroy(&Asp);CHKERRQ(ierr);
121   ierr = PetscFinalize();
122   return ierr;
123 }
124 
125