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