Actual source code: ex16.c

      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, NULL, 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: }