1 static char help[] = "Tests MatReorderForNonzeroDiagonal().\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** argv)5 int main(int argc, char **argv)
6 {
7 Mat mat, B, C;
8 PetscInt i, j;
9 PetscMPIInt size;
10 PetscScalar v;
11 IS isrow, iscol, identity;
12 PetscViewer viewer;
13
14 PetscFunctionBeginUser;
15 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
16
17 /* ------- Assemble matrix, --------- */
18
19 PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
20 PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, 4, 4));
21 PetscCall(MatSetFromOptions(mat));
22 PetscCall(MatSetUp(mat));
23
24 /* set anti-diagonal of matrix */
25 v = 1.0;
26 i = 0;
27 j = 3;
28 PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES));
29 v = 2.0;
30 i = 1;
31 j = 2;
32 PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES));
33 v = 3.0;
34 i = 2;
35 j = 1;
36 PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES));
37 v = 4.0;
38 i = 3;
39 j = 0;
40 PetscCall(MatSetValues(mat, 1, &i, 1, &j, &v, INSERT_VALUES));
41 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
42 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
43
44 PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD, &viewer));
45 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_DENSE));
46 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix\n"));
47 PetscCall(MatView(mat, viewer));
48
49 PetscCall(MatGetOrdering(mat, MATORDERINGNATURAL, &isrow, &iscol));
50
51 PetscCall(MatPermute(mat, isrow, iscol, &B));
52 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by identity\n"));
53 PetscCall(MatView(B, viewer));
54 PetscCall(MatDestroy(&B));
55
56 PetscCall(MatReorderForNonzeroDiagonal(mat, 1.e-8, isrow, iscol));
57 PetscCall(MatPermute(mat, isrow, iscol, &B));
58 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by identity + NonzeroDiagonal()\n"));
59 PetscCall(MatView(B, viewer));
60 PetscCall(PetscViewerASCIIPrintf(viewer, "Row permutation\n"));
61 PetscCall(ISView(isrow, viewer));
62 PetscCall(PetscViewerASCIIPrintf(viewer, "Column permutation\n"));
63 PetscCall(ISView(iscol, viewer));
64 PetscCall(MatDestroy(&B));
65
66 PetscCall(ISDestroy(&isrow));
67 PetscCall(ISDestroy(&iscol));
68
69 PetscCall(MatGetOrdering(mat, MATORDERINGND, &isrow, &iscol));
70 PetscCall(MatPermute(mat, isrow, iscol, &B));
71 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by ND\n"));
72 PetscCall(MatView(B, viewer));
73 PetscCall(MatDestroy(&B));
74 PetscCall(PetscViewerASCIIPrintf(viewer, "ND row permutation\n"));
75 PetscCall(ISView(isrow, viewer));
76 PetscCall(PetscViewerASCIIPrintf(viewer, "ND column permutation\n"));
77 PetscCall(ISView(iscol, viewer));
78
79 PetscCall(MatReorderForNonzeroDiagonal(mat, 1.e-8, isrow, iscol));
80 PetscCall(MatPermute(mat, isrow, iscol, &B));
81 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by ND + NonzeroDiagonal()\n"));
82 PetscCall(MatView(B, viewer));
83 PetscCall(MatDestroy(&B));
84 PetscCall(PetscViewerASCIIPrintf(viewer, "ND + NonzeroDiagonal() row permutation\n"));
85 PetscCall(ISView(isrow, viewer));
86 PetscCall(PetscViewerASCIIPrintf(viewer, "ND + NonzeroDiagonal() column permutation\n"));
87 PetscCall(ISView(iscol, viewer));
88
89 PetscCall(ISDestroy(&isrow));
90 PetscCall(ISDestroy(&iscol));
91
92 PetscCall(MatGetOrdering(mat, MATORDERINGRCM, &isrow, &iscol));
93 PetscCall(MatPermute(mat, isrow, iscol, &B));
94 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by RCM\n"));
95 PetscCall(MatView(B, viewer));
96 PetscCall(MatDestroy(&B));
97 PetscCall(PetscViewerASCIIPrintf(viewer, "RCM row permutation\n"));
98 PetscCall(ISView(isrow, viewer));
99 PetscCall(PetscViewerASCIIPrintf(viewer, "RCM column permutation\n"));
100 PetscCall(ISView(iscol, viewer));
101
102 PetscCall(MatReorderForNonzeroDiagonal(mat, 1.e-8, isrow, iscol));
103 PetscCall(MatPermute(mat, isrow, iscol, &B));
104 PetscCall(PetscViewerASCIIPrintf(viewer, "Original matrix permuted by RCM + NonzeroDiagonal()\n"));
105 PetscCall(MatView(B, viewer));
106 PetscCall(PetscViewerASCIIPrintf(viewer, "RCM + NonzeroDiagonal() row permutation\n"));
107 PetscCall(ISView(isrow, viewer));
108 PetscCall(PetscViewerASCIIPrintf(viewer, "RCM + NonzeroDiagonal() column permutation\n"));
109 PetscCall(ISView(iscol, viewer));
110 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
111 if (size == 1) {
112 PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE));
113 PetscCall(ISCreateStride(PETSC_COMM_SELF, 4, 0, 1, &identity));
114 PetscCall(MatPermute(B, identity, identity, &C));
115 PetscCall(MatConvert(C, MATSEQSBAIJ, MAT_INPLACE_MATRIX, &C));
116 PetscCall(MatDestroy(&C));
117 PetscCall(ISDestroy(&identity));
118 }
119 PetscCall(MatDestroy(&B));
120 /* Test MatLUFactor(); set diagonal as zeros as requested by PETSc matrix factorization */
121 for (i = 0; i < 4; i++) {
122 v = 0.0;
123 PetscCall(MatSetValues(mat, 1, &i, 1, &i, &v, INSERT_VALUES));
124 }
125 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
126 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
127 PetscCall(MatLUFactor(mat, isrow, iscol, NULL));
128
129 /* Free data structures */
130 PetscCall(ISDestroy(&isrow));
131 PetscCall(ISDestroy(&iscol));
132 PetscCall(MatDestroy(&mat));
133
134 PetscCall(PetscFinalize());
135 return 0;
136 }
137
138 /*TEST
139
140 test:
141
142 TEST*/
143