1: static char help[] = "Reads a matrix from PETSc binary file. Use for view or investigating matrix data structure. \n\n";
  2: /*
  3:  Example:
  4:       ./ex16 -f <matrix file> -a_mat_view draw -draw_pause -1
  5:       ./ex16 -f <matrix file> -a_mat_view ascii::ascii_info
  6:  */

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

 24:   PetscFunctionBeginUser;
 25:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 26:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

 28:   /* Determine files from which we read the linear systems. */
 29:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &flg));
 30:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");

 32:   /* Open binary file.  Note that we use FILE_MODE_READ to indicate
 33:      reading from this file. */
 34:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));

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

 46:   /* Get a sparse matrix Asp by dumping zero entries of A */
 47:   PetscCall(MatCreate(PETSC_COMM_WORLD, &Asp));
 48:   PetscCall(MatSetSizes(Asp, m, n, PETSC_DECIDE, PETSC_DECIDE));
 49:   PetscCall(MatSetOptionsPrefix(Asp, "asp_"));
 50:   PetscCall(MatSetFromOptions(Asp));
 51:   Dnnz = (PetscInt)matinfo.nz_used / m + 1;
 52:   Onnz = Dnnz / 2;
 53:   printf("Dnnz %d %d\n", Dnnz, Onnz);
 54:   PetscCall(MatSeqAIJSetPreallocation(Asp, Dnnz, NULL));
 55:   PetscCall(MatMPIAIJSetPreallocation(Asp, Dnnz, NULL, Onnz, NULL));
 56:   /* The allocation above is approximate so we must set this option to be permissive.
 57:    * Real code should preallocate exactly. */
 58:   PetscCall(MatSetOption(Asp, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));

 60:   /* Check zero rows */
 61:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 62:   nrows = 0;
 63:   for (row = rstart; row < rend; row++) {
 64:     PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
 65:     nnzA += ncols;
 66:     norm = 0.0;
 67:     for (j = 0; j < ncols; j++) {
 68:       val = PetscAbsScalar(vals[j]);
 69:       if (norm < val) norm = norm;
 70:       if (val > dtol) {
 71:         PetscCall(MatSetValues(Asp, 1, &row, 1, &cols[j], &vals[j], INSERT_VALUES));
 72:         nnzAsp++;
 73:       }
 74:     }
 75:     if (!norm) nrows++;
 76:     PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
 77:   }
 78:   PetscCall(MatAssemblyBegin(Asp, MAT_FINAL_ASSEMBLY));
 79:   PetscCall(MatAssemblyEnd(Asp, MAT_FINAL_ASSEMBLY));

 81:   percent = (PetscReal)nnzA * 100 / (m * n);
 82:   PetscCall(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));
 83:   percent = (PetscReal)nnzAsp * 100 / (m * n);
 84:   PetscCall(PetscPrintf(PETSC_COMM_SELF, " [%d] Matrix Asp nnzAsp %d, %g percent\n", rank, nnzAsp, percent));

 86:   /* investigate matcoloring for Asp */
 87:   PetscBool Asp_coloring = PETSC_FALSE;
 88:   PetscCall(PetscOptionsHasName(NULL, NULL, "-Asp_color", &Asp_coloring));
 89:   if (Asp_coloring) {
 90:     MatColoring   mc;
 91:     ISColoring    iscoloring;
 92:     MatFDColoring matfdcoloring;
 93:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create coloring of Asp...\n"));
 94:     PetscCall(MatColoringCreate(Asp, &mc));
 95:     PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
 96:     PetscCall(MatColoringSetFromOptions(mc));
 97:     PetscCall(MatColoringApply(mc, &iscoloring));
 98:     PetscCall(MatColoringDestroy(&mc));
 99:     PetscCall(MatFDColoringCreate(Asp, iscoloring, &matfdcoloring));
100:     PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
101:     PetscCall(MatFDColoringSetUp(Asp, iscoloring, matfdcoloring));
102:     /*PetscCall(MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD));*/
103:     PetscCall(ISColoringDestroy(&iscoloring));
104:     PetscCall(MatFDColoringDestroy(&matfdcoloring));
105:   }

107:   /* Write Asp in binary for study - see ~petsc/src/mat/tests/ex124.c */
108:   PetscBool Asp_write = PETSC_FALSE;
109:   PetscCall(PetscOptionsHasName(NULL, NULL, "-Asp_write", &Asp_write));
110:   if (Asp_write) {
111:     PetscViewer viewer;
112:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "Write Asp into file Asp.dat ...\n"));
113:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Asp.dat", FILE_MODE_WRITE, &viewer));
114:     PetscCall(MatView(Asp, viewer));
115:     PetscCall(PetscViewerDestroy(&viewer));
116:   }

118:   PetscCall(MatDestroy(&A));
119:   PetscCall(MatDestroy(&Asp));
120:   PetscCall(PetscFinalize());
121:   return 0;
122: }