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