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_info
  7:  */

  9: #include <petscmat.h>
 12: int main(int argc,char **args)
 13: {
 14:   Mat                   A,Asp;
 15:   PetscViewer           fd;               /* viewer */
 16:   char                  file[PETSC_MAX_PATH_LEN];     /* input file name */
 17:   PetscErrorCode        ierr;
 18:   PetscInt              m,n,rstart,rend;
 19:   PetscBool             flg;
 20:   PetscInt             row,ncols,j,nrows,nnzA=0,nnzAsp=0;
 21:   const PetscInt       *cols;
 22:   const PetscScalar    *vals;
 23:   PetscReal            norm,percent,val,dtol=1.e-16;
 24:   PetscMPIInt          rank;
 25:   MatInfo              matinfo;
 26:   PetscInt             Dnnz,Onnz;
 27: 

 29:   PetscInitialize(&argc,&args,(char *)0,help);
 30:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 32:   /* Determine files from which we read the linear systems. */
 33:   PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
 34:   if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");

 36:   /* Open binary file.  Note that we use FILE_MODE_READ to indicate
 37:      reading from this file. */
 38:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);

 40:   /* Load the matrix; then destroy the viewer. */
 41:   MatCreate(PETSC_COMM_WORLD,&A);
 42:   MatSetOptionsPrefix(A,"a_");
 43:   MatSetFromOptions(A);
 44:   MatLoad(A,fd);
 45:   PetscViewerDestroy(&fd);
 46:   MatGetSize(A,&m,&n);
 47:   MatGetInfo(A,MAT_LOCAL,&matinfo);
 48:   //printf("matinfo.nz_used %g\n",matinfo.nz_used);

 50:   /* Get a sparse matrix Asp by dumping zero entries of A */
 51:   MatCreate(PETSC_COMM_WORLD,&Asp);
 52:   MatSetSizes(Asp,m,n,PETSC_DECIDE,PETSC_DECIDE);
 53:   MatSetOptionsPrefix(Asp,"asp_");
 54:   MatSetFromOptions(Asp);
 55:   Dnnz  = (PetscInt)matinfo.nz_used/m + 1;
 56:   Onnz  = Dnnz/2;
 57:   printf("Dnnz %d %d\n",Dnnz,Onnz);
 58:   MatSeqAIJSetPreallocation(Asp,Dnnz,PETSC_NULL);
 59:   MatMPIAIJSetPreallocation(Asp,Dnnz,PETSC_NULL,Onnz,PETSC_NULL);
 60: 
 61:   /* Check zero rows */
 62:   MatGetOwnershipRange(A,&rstart,&rend);
 63:   nrows = 0;
 64:   for (row=rstart; row<rend; row++){
 65:     MatGetRow(A,row,&ncols,&cols,&vals);
 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:         MatSetValues(Asp,1,&row,1,&cols[j],&vals[j],INSERT_VALUES);
 73:         nnzAsp++;
 74:       }
 75:     }
 76:     if (!norm) nrows++;
 77:     MatRestoreRow(A,row,&ncols,&cols,&vals);
 78:   }
 79:   MatAssemblyBegin(Asp,MAT_FINAL_ASSEMBLY);
 80:   MatAssemblyEnd(Asp,MAT_FINAL_ASSEMBLY);
 81: 
 82:   percent=(PetscReal)nnzA*100/(m*n);
 83:   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);
 84:   percent=(PetscReal)nnzAsp*100/(m*n);
 85:   PetscPrintf(PETSC_COMM_SELF," [%d] Matrix Asp nnzAsp %d, %g percent\n",rank,nnzAsp,percent);

 87:   /* investigate matcoloring for Asp */
 88:   PetscBool     Asp_coloring = PETSC_FALSE;
 89:   PetscOptionsHasName(PETSC_NULL,"-Asp_color",&Asp_coloring);
 90:   if (Asp_coloring){
 91:     ISColoring    iscoloring;
 92:     MatFDColoring matfdcoloring;
 93:     PetscPrintf(PETSC_COMM_WORLD," Create coloring of Asp...\n");
 94:     MatGetColoring(Asp,MATCOLORINGSL,&iscoloring);
 95:     MatFDColoringCreate(Asp,iscoloring,&matfdcoloring);
 96:     MatFDColoringSetFromOptions(matfdcoloring);
 97:     //MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD);
 98:     ISColoringDestroy(&iscoloring);
 99:     MatFDColoringDestroy(&matfdcoloring);
100:   }

102:   /* Write Asp in binary for study - see ~petsc/src/mat/examples/tests/ex124.c */
103:   PetscBool Asp_write = PETSC_FALSE;
104:   PetscOptionsHasName(PETSC_NULL,"-Asp_write",&Asp_write);
105:   if (Asp_write){
106:     PetscViewer    viewer;
107:     PetscPrintf(PETSC_COMM_SELF,"Write Asp into file Asp.dat ...\n");
108:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,"Asp.dat",FILE_MODE_WRITE,&viewer);
109:     MatView(Asp,viewer);
110:     PetscViewerDestroy(&viewer);
111:   }

113:   MatDestroy(&A);
114:   MatDestroy(&Asp);
115:   PetscFinalize();
116:   return 0;
117: }