1 #include <petscksp.h>
2 #include <petscpc.h>
3 #include <petscviewer.h>
4
5 typedef struct {
6 PetscInt num_levels;
7 PetscInt *n_per_level;
8 Mat stiff;
9 Mat *ProlongationOps;
10 PetscBT *CFMarkers;
11 KSP kspHypre;
12 } *DataCompression;
13
14 PetscErrorCode Create1dLaplacian(PetscInt, Mat *);
15 PetscErrorCode DataCompExportMats(DataCompression);
16 PetscErrorCode DataCompDestroy(DataCompression);
17
main(int Argc,char ** Args)18 int main(int Argc, char **Args)
19 {
20 PetscInt n_nodes = 33;
21 Vec x, b;
22 PC pcHypre;
23 DataCompression data_comp;
24
25 PetscFunctionBeginUser;
26 PetscCall(PetscInitialize(&Argc, &Args, NULL, NULL));
27
28 PetscCall(PetscNew(&data_comp));
29
30 // Creating stiffness matrix
31 PetscCall(Create1dLaplacian(n_nodes, &data_comp->stiff));
32 PetscCall(PetscObjectSetName((PetscObject)data_comp->stiff, "Stiffness"));
33
34 // Set-up BoomerAMG PC to get Prolongation Operators and Coarse/Fine splittings
35 PetscCall(KSPCreate(PETSC_COMM_WORLD, &data_comp->kspHypre));
36 PetscCall(KSPSetType(data_comp->kspHypre, KSPRICHARDSON));
37 PetscCall(KSPGetPC(data_comp->kspHypre, &pcHypre));
38 PetscCall(PCSetType(pcHypre, PCHYPRE));
39 PetscCall(PCHYPRESetType(pcHypre, "boomeramg"));
40 PetscCall(PCSetFromOptions(pcHypre));
41 PetscCall(PCSetOperators(pcHypre, data_comp->stiff, data_comp->stiff));
42 PetscCall(PCSetUp(pcHypre));
43
44 PetscCall(MatCreateVecs(data_comp->stiff, &x, &b));
45 PetscCall(PCApply(pcHypre, x, b));
46 PetscCall(VecDestroy(&x));
47 PetscCall(VecDestroy(&b));
48
49 //Viewing the PC and Extracting the Prolongation Operators and CFMarkers
50 PetscCall(PCView(pcHypre, NULL));
51 PetscCall(PCGetInterpolations(pcHypre, &data_comp->num_levels, &data_comp->ProlongationOps));
52 PetscCall(PCHYPREGetCFMarkers(pcHypre, &data_comp->n_per_level, &data_comp->CFMarkers));
53
54 PetscCall(DataCompExportMats(data_comp));
55 PetscCall(DataCompDestroy(data_comp));
56 PetscCall(PetscFinalize());
57 return 0;
58 }
59
Create1dLaplacian(PetscInt n,Mat * mat)60 PetscErrorCode Create1dLaplacian(PetscInt n, Mat *mat)
61 {
62 PetscFunctionBeginUser;
63 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, mat));
64 PetscCall(MatSetValue(*mat, n - 1, n - 1, 2.0, INSERT_VALUES));
65 for (PetscInt i = 0; i < n - 1; i++) {
66 PetscCall(MatSetValue(*mat, i, i, 2.0, INSERT_VALUES));
67 PetscCall(MatSetValue(*mat, i + 1, i, -1.0, INSERT_VALUES));
68 PetscCall(MatSetValue(*mat, i, i + 1, -1.0, INSERT_VALUES));
69 }
70 PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
71 PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
72 PetscFunctionReturn(PETSC_SUCCESS);
73 }
74
DataCompExportMats(DataCompression data_comp)75 PetscErrorCode DataCompExportMats(DataCompression data_comp)
76 {
77 PetscFunctionBeginUser;
78 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Num levels: %" PetscInt_FMT "\n", data_comp->num_levels));
79 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " -- Nodes per level --\n"));
80 for (PetscInt i = 0; i < data_comp->num_levels; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Level %" PetscInt_FMT ": %" PetscInt_FMT "\n", i, data_comp->n_per_level[i]));
81
82 for (PetscInt i = 0; i < data_comp->num_levels - 1; i++) {
83 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Prolongation Operator - Level %" PetscInt_FMT "\n", i));
84 PetscCall(PetscObjectSetName((PetscObject)data_comp->ProlongationOps[i], "P"));
85 PetscCall(MatView(data_comp->ProlongationOps[i], PETSC_VIEWER_STDOUT_WORLD));
86 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n"));
87 }
88
89 for (PetscInt i = 0; i < data_comp->num_levels - 1; i++) {
90 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Coarse/Fine splitting - Level %" PetscInt_FMT "\n", i + 1));
91 PetscCall(PetscBTView(data_comp->n_per_level[i + 1], data_comp->CFMarkers[i], PETSC_VIEWER_STDOUT_WORLD));
92 }
93 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Stiffness matrix, sparse format:\n"));
94 PetscCall(MatViewFromOptions(data_comp->stiff, NULL, "-mat_view_stiff"));
95 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Finished calling the Viewer functions\n"));
96 PetscFunctionReturn(PETSC_SUCCESS);
97 }
98
DataCompDestroy(DataCompression data_comp)99 PetscErrorCode DataCompDestroy(DataCompression data_comp)
100 {
101 PetscFunctionBeginUser;
102 if (data_comp == NULL) PetscFunctionReturn(PETSC_SUCCESS);
103 PetscCall(MatDestroy(&data_comp->stiff));
104 PetscCall(KSPDestroy(&data_comp->kspHypre));
105 for (PetscInt i = 0; i < data_comp->num_levels - 1; i++) {
106 PetscCall(MatDestroy(&data_comp->ProlongationOps[i]));
107 PetscCall(PetscBTDestroy(&data_comp->CFMarkers[i]));
108 }
109 PetscCall(PetscFree(data_comp->ProlongationOps));
110 PetscCall(PetscFree(data_comp->n_per_level));
111 PetscCall(PetscFree(data_comp->CFMarkers));
112 PetscCall(PetscFree(data_comp));
113 PetscFunctionReturn(PETSC_SUCCESS);
114 }
115
116 /*TEST
117
118 test:
119 requires: hypre !complex !defined(PETSC_HAVE_HYPRE_DEVICE)
120 args: -pc_hypre_boomeramg_coarsen_type modifiedRuge-Stueben -pc_hypre_boomeramg_interp_type classical -pc_hypre_boomeramg_strong_threshold 0.25 pc_hypre_boomeramg_numfunctions 1 -pc_hypre_boomeramg_max_row_sum 1.0 -mat_view_stiff
121
122 TEST*/
123