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

  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:   PetscInt           m, n, rstart, rend;
 16:   PetscBool          flg;
 17:   PetscInt           row, ncols, j, nrows, nnzA = 0, nnzAsp = 0;
 18:   const PetscInt    *cols;
 19:   const PetscScalar *vals;
 20:   PetscReal          norm, percent, val, dtol = 1.e-16;
 21:   PetscMPIInt        rank;
 22:   MatInfo            matinfo;
 23:   PetscInt           Dnnz, Onnz;

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

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

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

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

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

 61:   /* Check zero rows */
 62:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
 63:   nrows = 0;
 64:   for (row = rstart; row < rend; row++) {
 65:     PetscCall(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:         PetscCall(MatSetValues(Asp, 1, &row, 1, &cols[j], &vals[j], INSERT_VALUES));
 73:         nnzAsp++;
 74:       }
 75:     }
 76:     if (!norm) nrows++;
 77:     PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
 78:   }
 79:   PetscCall(MatAssemblyBegin(Asp, MAT_FINAL_ASSEMBLY));
 80:   PetscCall(MatAssemblyEnd(Asp, MAT_FINAL_ASSEMBLY));

 82:   percent = (PetscReal)nnzA * 100 / (m * n);
 83:   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));
 84:   percent = (PetscReal)nnzAsp * 100 / (m * n);
 85:   PetscCall(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:   PetscCall(PetscOptionsHasName(NULL, NULL, "-Asp_color", &Asp_coloring));
 90:   if (Asp_coloring) {
 91:     MatColoring   mc;
 92:     ISColoring    iscoloring;
 93:     MatFDColoring matfdcoloring;
 94:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Create coloring of Asp...\n"));
 95:     PetscCall(MatColoringCreate(Asp, &mc));
 96:     PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
 97:     PetscCall(MatColoringSetFromOptions(mc));
 98:     PetscCall(MatColoringApply(mc, &iscoloring));
 99:     PetscCall(MatColoringDestroy(&mc));
100:     PetscCall(MatFDColoringCreate(Asp, iscoloring, &matfdcoloring));
101:     PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
102:     PetscCall(MatFDColoringSetUp(Asp, iscoloring, matfdcoloring));
103:     /*PetscCall(MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD));*/
104:     PetscCall(ISColoringDestroy(&iscoloring));
105:     PetscCall(MatFDColoringDestroy(&matfdcoloring));
106:   }

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

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