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;
22: PetscBool flg,test=PETSC_FALSE,reuse=PETSC_FALSE,viewexpl=PETSC_FALSE;
23: PetscScalar v;
24: PC pc;
26: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
27: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
28: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
29: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
30: PetscOptionsGetBool(NULL,NULL,"-test_hmg_interface",&test,NULL);
31: PetscOptionsGetBool(NULL,NULL,"-test_reuse_interpolation",&reuse,NULL);
32: PetscOptionsGetBool(NULL,NULL,"-view_explicit_mat",&viewexpl,NULL);
34: MatCreate(PETSC_COMM_WORLD,&A);
35: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n*bs,m*n*bs);
36: MatSetBlockSize(A,bs);
37: MatSetFromOptions(A);
38: MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
39: MatSeqAIJSetPreallocation(A,5,NULL);
40: #if defined(PETSC_HAVE_HYPRE)
41: MatHYPRESetPreallocation(A,5,NULL,5,NULL);
42: #endif
44: MatGetOwnershipRange(A,&Istart,&Iend);
46: for (Ii=Istart/bs; Ii<Iend/bs; Ii++) {
47: v = -1.0; i = Ii/n; j = Ii - i*n;
48: if (i>0) {
49: J = Ii - n;
50: for (jj=0; jj<bs; jj++) {
51: II = Ii*bs + jj;
52: JJ = J*bs + jj;
53: MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
54: }
55: }
56: if (i<m-1) {
57: J = Ii + n;
58: for (jj=0; jj<bs; jj++) {
59: II = Ii*bs + jj;
60: JJ = J*bs + jj;
61: MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
62: }
63: }
64: if (j>0) {
65: J = Ii - 1;
66: for (jj=0; jj<bs; jj++) {
67: II = Ii*bs + jj;
68: JJ = J*bs + jj;
69: MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
70: }
71: }
72: if (j<n-1) {
73: J = Ii + 1;
74: for (jj=0; jj<bs; jj++) {
75: II = Ii*bs + jj;
76: JJ = J*bs + jj;
77: MatSetValues(A,1,&II,1,&JJ,&v,ADD_VALUES);
78: }
79: }
80: v = 4.0;
81: for (jj=0; jj<bs; jj++) {
82: II = Ii*bs + jj;
83: MatSetValues(A,1,&II,1,&II,&v,ADD_VALUES);
84: }
85: }
87: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
89: if (viewexpl) {
90: Mat E;
91: MatComputeOperator(A,MATAIJ,&E);
92: MatView(E,NULL);
93: MatDestroy(&E);
94: }
96: MatCreateVecs(A,&u,NULL);
97: VecSetFromOptions(u);
98: VecDuplicate(u,&b);
99: VecDuplicate(b,&x);
101: VecSet(u,1.0);
102: MatMult(A,u,b);
104: flg = PETSC_FALSE;
105: PetscOptionsGetBool(NULL,NULL,"-view_exact_sol",&flg,NULL);
106: if (flg) {VecView(u,PETSC_VIEWER_STDOUT_WORLD);}
108: KSPCreate(PETSC_COMM_WORLD,&ksp);
109: KSPSetOperators(ksp,A,A);
110: KSPSetTolerances(ksp,1.e-2/((m+1)*(n+1)),1.e-50,PETSC_DEFAULT,PETSC_DEFAULT);
112: if (test) {
113: KSPGetPC(ksp,&pc);
114: PCSetType(pc,PCHMG);
115: PCHMGSetInnerPCType(pc,PCGAMG);
116: PCHMGSetReuseInterpolation(pc,PETSC_TRUE);
117: PCHMGSetUseSubspaceCoarsening(pc,PETSC_TRUE);
118: PCHMGUseMatMAIJ(pc,PETSC_FALSE);
119: PCHMGSetCoarseningComponent(pc,0);
120: }
122: KSPSetFromOptions(ksp);
123: KSPSolve(ksp,b,x);
125: if (reuse) {
126: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
127: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
128: KSPSolve(ksp,b,x);
129: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
130: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
131: KSPSolve(ksp,b,x);
132: /* Make sparsity pattern different and reuse interpolation */
133: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
134: MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);
135: MatGetSize(A,&m,NULL);
136: n = 0;
137: v = 0;
138: m--;
139: /* Connect the last element to the first element */
140: MatSetValue(A,m,n,v,ADD_VALUES);
141: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
142: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
143: KSPSolve(ksp,b,x);
144: }
146: VecAXPY(x,-1.0,u);
147: VecNorm(x,NORM_2,&norm);
148: KSPGetIterationNumber(ksp,&its);
150: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);
152: KSPDestroy(&ksp);
153: VecDestroy(&u);
154: VecDestroy(&x);
155: VecDestroy(&b);
156: MatDestroy(&A);
158: PetscFinalize();
159: return ierr;
160: }
162: /*TEST
164: build:
165: requires: !complex !single
167: test:
168: suffix: hypre
169: nsize: 2
170: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
171: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre
173: test:
174: suffix: hypre_bs4
175: nsize: 2
176: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
177: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -pc_hmg_use_subspace_coarsening 1
179: test:
180: suffix: hypre_asm
181: nsize: 2
182: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
183: 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
185: test:
186: suffix: hypre_fieldsplit
187: nsize: 2
188: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
189: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type hypre -bs 4 -mg_levels_4_pc_type fieldsplit
191: test:
192: suffix: gamg
193: nsize: 2
194: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg
196: test:
197: suffix: gamg_bs4
198: nsize: 2
199: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -pc_hmg_use_subspace_coarsening 1
201: test:
202: suffix: gamg_asm
203: nsize: 2
204: 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
206: test:
207: suffix: gamg_fieldsplit
208: nsize: 2
209: args: -ksp_monitor -pc_type hmg -ksp_rtol 1e-6 -hmg_inner_pc_type gamg -bs 4 -mg_levels_1_pc_type fieldsplit
211: test:
212: suffix: interface
213: nsize: 2
214: args: -ksp_monitor -ksp_rtol 1e-6 -test_hmg_interface 1 -bs 4
216: test:
217: suffix: reuse
218: nsize: 2
219: args: -ksp_monitor -ksp_rtol 1e-6 -pc_type hmg -pc_hmg_reuse_interpolation 1 -test_reuse_interpolation 1 -hmg_inner_pc_type gamg
221: test:
222: suffix: component
223: nsize: 2
224: 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
226: testset:
227: output_file: output/ex4_expl.out
228: nsize: {{1 2}}
229: filter: grep -v "MPI processes" | grep -v " type:" | grep -v "Mat Object"
230: args: -ksp_converged_reason -view_explicit_mat -pc_type none -ksp_type {{cg gmres}}
231: test:
232: suffix: expl_aij
233: args: -mat_type aij
234: test:
235: suffix: expl_hypre
236: requires: hypre
237: args: -mat_type hypre
239: test:
240: suffix: hypre_device
241: nsize: {{1 2}}
242: requires: hypre defined(PETSC_HAVE_HYPRE_DEVICE)
243: args: -mat_type hypre -ksp_converged_reason -pc_type hypre -m 13 -n 17
245: test:
246: suffix: hypre_device_cusparse
247: output_file: output/ex4_hypre_device.out
248: nsize: {{1 2}}
249: requires: hypre cuda defined(PETSC_HAVE_HYPRE_DEVICE)
250: args: -mat_type {{aij aijcusparse}} -vec_type cuda -ksp_converged_reason -pc_type hypre -m 13 -n 17
252: TEST*/