Actual source code: ex3.c
1: static char help[] = "Bilinear elements on the unit square for Laplacian. To test the parallel\n\
2: matrix assembly, the matrix is intentionally laid out across processors\n\
3: differently from the way it is assembled. Input arguments are:\n\
4: -m <size> : problem size\n\n";
6: /* Addendum: piggy-backing on this example to test KSPChebyshev methods */
8: #include <petscksp.h>
10: PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
11: {
12: PetscFunctionBeginUser;
13: Ke[0] = H / 6.0;
14: Ke[1] = -.125 * H;
15: Ke[2] = H / 12.0;
16: Ke[3] = -.125 * H;
17: Ke[4] = -.125 * H;
18: Ke[5] = H / 6.0;
19: Ke[6] = -.125 * H;
20: Ke[7] = H / 12.0;
21: Ke[8] = H / 12.0;
22: Ke[9] = -.125 * H;
23: Ke[10] = H / 6.0;
24: Ke[11] = -.125 * H;
25: Ke[12] = -.125 * H;
26: Ke[13] = H / 12.0;
27: Ke[14] = -.125 * H;
28: Ke[15] = H / 6.0;
29: PetscFunctionReturn(PETSC_SUCCESS);
30: }
32: PetscErrorCode FormElementRhs(PetscReal x, PetscReal y, PetscReal H, PetscScalar *r)
33: {
34: PetscFunctionBeginUser;
35: r[0] = 0.;
36: r[1] = 0.;
37: r[2] = 0.;
38: r[3] = 0.0;
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: int main(int argc, char **args)
43: {
44: Mat C;
45: PetscMPIInt rank, size;
46: PetscInt i, m = 5, N, start, end, M, its;
47: PetscScalar val, Ke[16], r[4];
48: PetscReal x, y, h, norm;
49: PetscInt idx[4], count, *rows;
50: Vec u, ustar, b, build_sol;
51: KSP ksp;
52: PetscBool viewkspest = PETSC_FALSE, testbuildsolution = PETSC_FALSE, ishypre;
53: PC pc;
55: PetscFunctionBeginUser;
56: PetscCall(PetscInitialize(&argc, &args, NULL, help));
57: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
58: PetscCall(PetscOptionsGetBool(NULL, NULL, "-ksp_est_view", &viewkspest, NULL));
59: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_build_solution", &testbuildsolution, NULL));
60: N = (m + 1) * (m + 1); /* dimension of matrix */
61: M = m * m; /* number of elements */
62: h = 1.0 / m; /* mesh width */
63: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
64: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
66: /* Create stiffness matrix */
67: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
68: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N));
69: PetscCall(MatSetFromOptions(C));
70: PetscCall(MatSetUp(C));
71: start = rank * (M / size) + ((M % size) < rank ? (M % size) : rank);
72: end = start + M / size + ((M % size) > rank);
74: /* Assemble matrix */
75: PetscCall(FormElementStiffness(h * h, Ke)); /* element stiffness for Laplacian */
76: for (i = start; i < end; i++) {
77: /* node numbers for the four corners of element */
78: idx[0] = (m + 1) * (i / m) + (i % m);
79: idx[1] = idx[0] + 1;
80: idx[2] = idx[1] + m + 1;
81: idx[3] = idx[2] - 1;
82: PetscCall(MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES));
83: }
84: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
85: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
87: /* Create right-hand side and solution vectors */
88: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
89: PetscCall(VecSetSizes(u, PETSC_DECIDE, N));
90: PetscCall(VecSetFromOptions(u));
91: PetscCall(PetscObjectSetName((PetscObject)u, "Approx. Solution"));
92: PetscCall(VecDuplicate(u, &b));
93: PetscCall(PetscObjectSetName((PetscObject)b, "Right hand side"));
94: PetscCall(VecDuplicate(b, &ustar));
96: /* Assemble right-hand-side vector */
97: for (i = start; i < end; i++) {
98: /* location of lower left corner of element */
99: x = h * (i % m);
100: y = h * (i / m);
101: /* node numbers for the four corners of element */
102: idx[0] = (m + 1) * (i / m) + (i % m);
103: idx[1] = idx[0] + 1;
104: idx[2] = idx[1] + m + 1;
105: idx[3] = idx[2] - 1;
106: PetscCall(FormElementRhs(x, y, h * h, r));
107: PetscCall(VecSetValues(b, 4, idx, r, ADD_VALUES));
108: }
109: PetscCall(VecAssemblyBegin(b));
110: PetscCall(VecAssemblyEnd(b));
112: /* Modify matrix and right-hand side for Dirichlet boundary conditions */
113: PetscCall(PetscMalloc1(4 * m, &rows));
114: for (i = 0; i < m + 1; i++) {
115: rows[i] = i; /* bottom */
116: rows[3 * m - 1 + i] = m * (m + 1) + i; /* top */
117: }
118: count = m + 1; /* left side */
119: for (i = m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
121: count = 2 * m; /* left side */
122: for (i = 2 * m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
123: for (i = 0; i < 4 * m; i++) {
124: val = h * (rows[i] / (m + 1));
125: PetscCall(VecSetValues(u, 1, &rows[i], &val, INSERT_VALUES));
126: PetscCall(VecSetValues(b, 1, &rows[i], &val, INSERT_VALUES));
127: }
128: PetscCall(MatZeroRows(C, 4 * m, rows, 1.0, NULL, NULL));
130: PetscCall(PetscFree(rows));
131: PetscCall(VecAssemblyBegin(u));
132: PetscCall(VecAssemblyEnd(u));
133: PetscCall(VecAssemblyBegin(b));
134: PetscCall(VecAssemblyEnd(b));
136: {
137: Mat A;
138: PetscCall(MatConvert(C, MATSAME, MAT_INITIAL_MATRIX, &A));
139: PetscCall(MatDestroy(&C));
140: PetscCall(MatConvert(A, MATSAME, MAT_INITIAL_MATRIX, &C));
141: PetscCall(MatDestroy(&A));
142: }
144: /* Solve linear system */
145: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
146: PetscCall(KSPSetOperators(ksp, C, C));
147: PetscCall(KSPSetFromOptions(ksp));
148: PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
150: /* verify that PCView_HYPRE() handles PETSC_DECIDE parameters correctly */
151: PetscCall(KSPGetPC(ksp, &pc));
152: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCHYPRE, &ishypre));
153: if (ishypre) PetscCall(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD));
155: PetscCall(KSPSolve(ksp, b, u));
157: if (testbuildsolution) {
158: PetscBool ok;
160: PetscCall(VecDuplicate(u, &build_sol));
161: PetscCall(KSPBuildSolution(ksp, build_sol, NULL));
162: PetscCall(VecEqual(u, build_sol, &ok));
163: PetscCheck(ok, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "KSPBuildSolution() returned incorrect solution");
164: PetscCall(VecDestroy(&build_sol));
165: }
167: if (viewkspest) {
168: KSP kspest;
170: PetscCall(KSPChebyshevEstEigGetKSP(ksp, &kspest));
171: if (kspest) PetscCall(KSPView(kspest, PETSC_VIEWER_STDOUT_WORLD));
172: }
174: /* Check error */
175: PetscCall(VecGetOwnershipRange(ustar, &start, &end));
176: for (i = start; i < end; i++) {
177: val = h * (i / (m + 1));
178: PetscCall(VecSetValues(ustar, 1, &i, &val, INSERT_VALUES));
179: }
180: PetscCall(VecAssemblyBegin(ustar));
181: PetscCall(VecAssemblyEnd(ustar));
182: PetscCall(VecAXPY(u, -1.0, ustar));
183: PetscCall(VecNorm(u, NORM_2, &norm));
184: PetscCall(KSPGetIterationNumber(ksp, &its));
185: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)(norm * h), its));
187: /* Free work space */
188: PetscCall(KSPDestroy(&ksp));
189: PetscCall(VecDestroy(&ustar));
190: PetscCall(VecDestroy(&u));
191: PetscCall(VecDestroy(&b));
192: PetscCall(MatDestroy(&C));
193: PetscCall(PetscFinalize());
194: return 0;
195: }
197: /*TEST
199: test:
200: suffix: vector_allocate_1
201: filter: grep -w -A15 "Destructions" | grep -w -m 1 "Vector" | awk "{print ((\$2 + 0) <= 20)}"
202: requires: defined(PETSC_USE_LOG)
203: args: -pc_type jacobi -m 5 -ksp_orthogonalization_cgs_refinement_type refine_always -ksp_gmres_vector_allocation_size 1 -log_view
205: test:
206: suffix: vector_allocate_2
207: filter: grep -w -A15 "Destructions" | grep -w -m 1 "Vector" | awk "{print ((\$2 + 0) > 20)}"
208: requires: defined(PETSC_USE_LOG)
209: args: -pc_type jacobi -m 5 -ksp_orthogonalization_cgs_refinement_type refine_always -ksp_gmres_vector_allocation_size 20 -log_view
211: test:
212: args: -pc_type jacobi -ksp_monitor -m 5 -ksp_orthogonalization_cgs_refinement_type refine_always
214: test:
215: suffix: 2
216: nsize: 2
217: args: -pc_type jacobi -ksp_monitor -m 5 -ksp_orthogonalization_cgs_refinement_type refine_always
219: test:
220: suffix: 2_kokkos
221: nsize: 2
222: args: -vec_mdot_use_gemv {{0 1}} -vec_maxpy_use_gemv {{0 1}}
223: args: -pc_type jacobi -ksp_monitor -m 5 -ksp_orthogonalization_cgs_refinement_type refine_always -mat_type aijkokkos -vec_type kokkos
224: output_file: output/ex3_2.out
225: requires: kokkos_kernels
227: test:
228: suffix: nocheby
229: args: -ksp_est_view
231: test:
232: suffix: chebynoest
233: args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_eigenvalues 0.1,1.0
235: test:
236: suffix: chebyest
237: args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_esteig
238: filter: sed -e "s/Iterations 19/Iterations 20/g"
240: test:
241: suffix: gamg_provided_not_ok
242: filter: grep -v "variant HERMITIAN" | sed -e "s/Iterations 4/Iterations 5/g"
243: args: -pc_type gamg -mg_levels_pc_type sor -mg_levels_esteig_ksp_type cg -ksp_view
245: test:
246: suffix: build_solution
247: requires: !complex
248: filter: grep -v Norm
249: args: -ksp_type {{chebyshev cg groppcg pipecg pipecgrr pipelcg pipeprcg cgne nash stcg gltr fcg pipefcg gmres fgmres lgmres dgmres pgmres tcqmr bcgs ibcgs qmrcgs fbcgs fbcgsr bcgsl pipebcgs cgs tfqmr cr pipecr bicg minres lcd gcr cgls richardson}} -test_build_solution
250: output_file: output/empty.out
252: test:
253: suffix: hypre
254: requires: hypre
255: args: -pc_type hypre
257: test:
258: suffix: caliper
259: nsize: 2
260: requires: hypre caliper
261: env: CALI_CONFIG=runtime-report,max_column_width=200,calc.inclusive,mpi-report,output=stdout
262: args: -pc_type hypre
263: filter: grep "Min time/rank Max time/rank"
265: test:
266: suffix: pflare
267: requires: !complex !single pflare
268: args: -pc_type air -ksp_max_it 5 -m 10
269: filter: grep -v .
270: output_file: output/empty.out
272: test:
273: suffix: 2_pflare
274: nsize: 2
275: requires: !complex !single pflare
276: args: -pc_type air -ksp_max_it 5 -m 10
277: filter: grep -v .
278: output_file: output/empty.out
280: TEST*/