Actual source code: ex4.c
2: static char help[] = "Solves a linear system in parallel with KSP and HMG.\n\
3: Input parameters include:\n\
4: -view_exact_sol : write exact solution vector to stdout\n\
5: -m <mesh_x> : number of mesh points in x-direction\n\
6: -n <mesh_y> : number of mesh points in y-direction\n\
7: -bs : number of variables on each mesh vertex \n\n";
9: /*
10: Simple example is used to test PCHMG
11: */
12: #include <petscksp.h>
14: int main(int argc,char **args)
15: {
16: Vec x,b,u; /* approx solution, RHS, exact solution */
17: Mat A; /* linear system matrix */
18: KSP ksp; /* linear solver context */
19: PetscReal norm; /* norm of solution error */
20: PetscInt i,j,Ii,J,Istart,Iend,m = 8,n = 7,its,bs=1,II,JJ,jj;
21: PetscBool flg,test=PETSC_FALSE,reuse=PETSC_FALSE,viewexpl=PETSC_FALSE;
22: PetscScalar v;
23: PC pc;
25: PetscInitialize(&argc,&args,(char*)0,help);
26: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
27: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
28: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
29: PetscOptionsGetBool(NULL,NULL,"-test_hmg_interface",&test,NULL);
30: PetscOptionsGetBool(NULL,NULL,"-test_reuse_interpolation",&reuse,NULL);
31: PetscOptionsGetBool(NULL,NULL,"-view_explicit_mat",&viewexpl,NULL);
33: MatCreate(PETSC_COMM_WORLD,&A);
34: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n*bs,m*n*bs);
35: MatSetBlockSize(A,bs);
36: MatSetFromOptions(A);
37: MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
38: MatSeqAIJSetPreallocation(A,5,NULL);
39: #if defined(PETSC_HAVE_HYPRE)
40: MatHYPRESetPreallocation(A,5,NULL,5,NULL);
41: #endif
43: MatGetOwnershipRange(A,&Istart,&Iend);
45: for (Ii=Istart/bs; Ii<Iend/bs; Ii++) {
46: v = -1.0; i = Ii/n; j = Ii - i*n;
47: if (i>0) {
48: J = Ii - n;
49: for (jj=0; jj<bs; jj++) {
50: II = Ii*bs + jj;
51: JJ = J*bs + jj;
52: MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
53: }
54: }
55: if (i<m-1) {
56: J = Ii + n;
57: for (jj=0; jj<bs; jj++) {
58: II = Ii*bs + jj;
59: JJ = J*bs + jj;
60: MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
61: }
62: }
63: if (j>0) {
64: J = Ii - 1;
65: for (jj=0; jj<bs; jj++) {
66: II = Ii*bs + jj;
67: JJ = J*bs + jj;
68: MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
69: }
70: }
71: if (j<n-1) {
72: J = Ii + 1;
73: for (jj=0; jj<bs; jj++) {
74: II = Ii*bs + jj;
75: JJ = J*bs + jj;
76: MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
77: }
78: }
79: v = 4.0;
80: for (jj=0; jj<bs; jj++) {
81: II = Ii*bs + jj;
82: MatSetValues(A,1,&II,1,&II,&v,ADD_VALUES);
83: }
84: }
86: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
87: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
88: if (viewexpl) {
89: Mat E;
90: MatComputeOperator(A,MATAIJ,&E);
91: MatView(E,NULL);
92: MatDestroy(&E);
93: }
95: MatCreateVecs(A,&u,NULL);
96: VecSetFromOptions(u);
97: VecDuplicate(u,&b);
98: VecDuplicate(b,&x);
100: VecSet(u,1.0);
101: MatMult(A,u,b);
103: flg = PETSC_FALSE;
104: PetscOptionsGetBool(NULL,NULL,"-view_exact_sol",&flg,NULL);
105: if (flg) VecView(u,PETSC_VIEWER_STDOUT_WORLD);
107: KSPCreate(PETSC_COMM_WORLD,&ksp);
108: KSPSetOperators(ksp,A,A);
109: KSPSetTolerances(ksp,1.e-2/((m+1)*(n+1)),1.e-50,PETSC_DEFAULT,PETSC_DEFAULT);
111: if (test) {
112: KSPGetPC(ksp,&pc);
113: PCSetType(pc,PCHMG);
114: PCHMGSetInnerPCType(pc,PCGAMG);
115: PCHMGSetReuseInterpolation(pc,PETSC_TRUE);
116: PCHMGSetUseSubspaceCoarsening(pc,PETSC_TRUE);
117: PCHMGUseMatMAIJ(pc,PETSC_FALSE);
118: PCHMGSetCoarseningComponent(pc,0);
119: }
121: KSPSetFromOptions(ksp);
122: KSPSolve(ksp,b,x);
124: if (reuse) {
125: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
126: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
127: KSPSolve(ksp,b,x);
128: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
129: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
130: KSPSolve(ksp,b,x);
131: /* Make sparsity pattern different and reuse interpolation */
132: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
133: MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);
134: MatGetSize(A,&m,NULL);
135: n = 0;
136: v = 0;
137: m--;
138: /* Connect the last element to the first element */
139: MatSetValue(A,m,n,v,ADD_VALUES);
140: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
141: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
142: KSPSolve(ksp,b,x);
143: }
145: VecAXPY(x,-1.0,u);
146: VecNorm(x,NORM_2,&norm);
147: KSPGetIterationNumber(ksp,&its);
149: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);
151: KSPDestroy(&ksp);
152: VecDestroy(&u);
153: VecDestroy(&x);
154: VecDestroy(&b);
155: MatDestroy(&A);
157: PetscFinalize();
158: return 0;
159: }
161: /*TEST
163: build:
164: requires: !complex !single
166: test:
167: suffix: hypre
168: nsize: 2
169: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
170: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre
172: test:
173: suffix: hypre_bs4
174: nsize: 2
175: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
176: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -pc_hmg_use_subspace_coarsening 1
178: test:
179: suffix: hypre_asm
180: nsize: 2
181: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
182: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -pc_hmg_use_subspace_coarsening 1 -mg_levels_3_pc_type asm
184: test:
185: suffix: hypre_fieldsplit
186: nsize: 2
187: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
188: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -mg_levels_4_pc_type fieldsplit
190: test:
191: suffix: gamg
192: nsize: 2
193: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg
195: test:
196: suffix: gamg_bs4
197: nsize: 2
198: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -pc_hmg_use_subspace_coarsening 1
200: test:
201: suffix: gamg_asm
202: nsize: 2
203: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -pc_hmg_use_subspace_coarsening 1 -mg_levels_1_pc_type asm
205: test:
206: suffix: gamg_fieldsplit
207: nsize: 2
208: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -mg_levels_1_pc_type fieldsplit
210: test:
211: suffix: interface
212: nsize: 2
213: args: -ksp_monitor -ksp_rtol 1e-6 -test_hmg_interface 1 -bs 4
215: test:
216: suffix: reuse
217: nsize: 2
218: args: -ksp_monitor -ksp_rtol 1e-6 -pc_type hmg -pc_hmg_reuse_interpolation 1 -test_reuse_interpolation 1 -hmg_inner_pc_type gamg
220: test:
221: suffix: component
222: nsize: 2
223: args: -ksp_monitor -ksp_rtol 1e-6 -pc_type hmg -pc_hmg_coarsening_component 2 -pc_hmg_use_subspace_coarsening 1 -bs 4 -hmg_inner_pc_type gamg
225: testset:
226: output_file: output/ex4_expl.out
227: nsize: {{1 2}}
228: filter: grep -v "MPI processes" | grep -v " type:" | grep -v "Mat Object"
229: args: -ksp_converged_reason -view_explicit_mat -pc_type none -ksp_type {{cg gmres}}
230: test:
231: suffix: expl_aij
232: args: -mat_type aij
233: test:
234: suffix: expl_hypre
235: requires: hypre
236: args: -mat_type hypre
238: test:
239: suffix: hypre_device
240: nsize: {{1 2}}
241: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
242: args: -mat_type hypre -ksp_converged_reason -pc_type hypre -m 13 -n 17
244: test:
245: suffix: hypre_device_cusparse
246: output_file: output/ex4_hypre_device.out
247: nsize: {{1 2}}
248: requires: hypre cuda defined(PETSC_HAVE_HYPRE_DEVICE)
249: args: -mat_type {{aij aijcusparse}} -vec_type cuda -ksp_converged_reason -pc_type hypre -m 13 -n 17
251: TEST*/