Actual source code: ex68.c
1: static char help[] = "Tests MatReorderForNonzeroDiagonal().\n\n";
3: #include <petscmat.h>
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;
14: PetscFunctionBeginUser;
15: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
17: /* ------- Assemble matrix, --------- */
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));
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));
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));
49: PetscCall(MatGetOrdering(mat, MATORDERINGNATURAL, &isrow, &iscol));
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));
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));
66: PetscCall(ISDestroy(&isrow));
67: PetscCall(ISDestroy(&iscol));
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));
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));
89: PetscCall(ISDestroy(&isrow));
90: PetscCall(ISDestroy(&iscol));
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));
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));
129: /* Free data structures */
130: PetscCall(ISDestroy(&isrow));
131: PetscCall(ISDestroy(&iscol));
132: PetscCall(MatDestroy(&mat));
134: PetscCall(PetscFinalize());
135: return 0;
136: }
138: /*TEST
140: test:
142: TEST*/