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, 0, 0));
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: args: -pc_type jacobi -ksp_monitor -m 5 -ksp_orthogonalization_cgs_refinement_type refine_always
202: test:
203: suffix: 2
204: nsize: 2
205: args: -pc_type jacobi -ksp_monitor -m 5 -ksp_orthogonalization_cgs_refinement_type refine_always
207: test:
208: suffix: 2_kokkos
209: nsize: 2
210: args: -vec_mdot_use_gemv {{0 1}} -vec_maxpy_use_gemv {{0 1}}
211: args: -pc_type jacobi -ksp_monitor -m 5 -ksp_orthogonalization_cgs_refinement_type refine_always -mat_type aijkokkos -vec_type kokkos
212: output_file: output/ex3_2.out
213: requires: kokkos_kernels
215: test:
216: suffix: nocheby
217: args: -ksp_est_view
219: test:
220: suffix: chebynoest
221: args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_eigenvalues 0.1,1.0
223: test:
224: suffix: chebyest
225: args: -ksp_est_view -ksp_type chebyshev -ksp_chebyshev_esteig
226: filter: sed -e "s/Iterations 19/Iterations 20/g"
228: test:
229: suffix: gamg_provided_not_ok
230: filter: grep -v "variant HERMITIAN" | sed -e "s/Iterations 4/Iterations 5/g"
231: args: -pc_type gamg -mg_levels_pc_type sor -mg_levels_esteig_ksp_type cg -ksp_view
233: test:
234: suffix: build_solution
235: requires: !complex
236: filter: grep -v Norm
237: 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
238: output_file: output/empty.out
240: test:
241: suffix: hypre
242: requires: hypre
243: args: -pc_type hypre
245: test:
246: suffix: caliper
247: nsize: 2
248: requires: hypre caliper
249: env: CALI_CONFIG=runtime-report,max_column_width=200,calc.inclusive,mpi-report,output=stdout
250: args: -pc_type hypre
251: filter: grep "Min time/rank Max time/rank"
253: test:
254: suffix: pflare
255: requires: !complex !single pflare
256: args: -pc_type air -ksp_max_it 5 -m 10
257: filter: grep -v .
258: output_file: output/empty.out
260: test:
261: suffix: 2_pflare
262: nsize: 2
263: requires: !complex !single pflare
264: args: -pc_type air -ksp_max_it 5 -m 10
265: filter: grep -v .
266: output_file: output/empty.out
268: TEST*/