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 */
7
8 #include <petscmat.h>
main(int argc,char ** args)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;
23
24 PetscFunctionBeginUser;
25 PetscCall(PetscInitialize(&argc, &args, NULL, help));
26 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
27
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");
31
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));
35
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);*/
45
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));
59
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));
80
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));
85
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 }
106
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 }
117
118 PetscCall(MatDestroy(&A));
119 PetscCall(MatDestroy(&Asp));
120 PetscCall(PetscFinalize());
121 return 0;
122 }
123