Actual source code: ex17.c
1: static char help[] = "Solves a linear system with KSP. This problem is\n\
2: intended to test the complex numbers version of various solvers.\n\n";
4: #include <petscksp.h>
6: typedef enum {
7: TEST_1,
8: TEST_2,
9: TEST_3,
10: HELMHOLTZ_1,
11: HELMHOLTZ_2
12: } TestType;
13: extern PetscErrorCode FormTestMatrix(Mat, PetscInt, TestType);
15: int main(int argc, char **args)
16: {
17: Vec x, b, u; /* approx solution, RHS, exact solution */
18: Mat A; /* linear system matrix */
19: KSP ksp; /* KSP context */
20: PetscInt n = 10, its, dim, p = 1, use_random;
21: PetscScalar none = -1.0, pfive = 0.5;
22: PetscReal norm;
23: PetscRandom rctx;
24: TestType type;
25: PetscBool flg;
27: PetscFunctionBeginUser;
28: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
29: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
30: PetscCall(PetscOptionsGetInt(NULL, NULL, "-p", &p, NULL));
31: switch (p) {
32: case 1:
33: type = TEST_1;
34: dim = n;
35: break;
36: case 2:
37: type = TEST_2;
38: dim = n;
39: break;
40: case 3:
41: type = TEST_3;
42: dim = n;
43: break;
44: case 4:
45: type = HELMHOLTZ_1;
46: dim = n * n;
47: break;
48: case 5:
49: type = HELMHOLTZ_2;
50: dim = n * n;
51: break;
52: default:
53: type = TEST_1;
54: dim = n;
55: }
57: /* Create vectors */
58: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
59: PetscCall(VecSetSizes(x, PETSC_DECIDE, dim));
60: PetscCall(VecSetFromOptions(x));
61: PetscCall(VecDuplicate(x, &b));
62: PetscCall(VecDuplicate(x, &u));
64: use_random = 1;
65: flg = PETSC_FALSE;
66: PetscCall(PetscOptionsGetBool(NULL, NULL, "-norandom", &flg, NULL));
67: if (flg) {
68: use_random = 0;
69: PetscCall(VecSet(u, pfive));
70: } else {
71: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
72: PetscCall(PetscRandomSetFromOptions(rctx));
73: PetscCall(VecSetRandom(u, rctx));
74: }
76: /* Create and assemble matrix */
77: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
78: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim));
79: PetscCall(MatSetFromOptions(A));
80: PetscCall(MatSetUp(A));
81: PetscCall(FormTestMatrix(A, n, type));
82: PetscCall(MatMult(A, u, b));
83: flg = PETSC_FALSE;
84: PetscCall(PetscOptionsGetBool(NULL, NULL, "-printout", &flg, NULL));
85: if (flg) {
86: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
87: PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));
88: PetscCall(VecView(b, PETSC_VIEWER_STDOUT_WORLD));
89: }
91: /* Create KSP context; set operators and options; solve linear system */
92: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
93: PetscCall(KSPSetOperators(ksp, A, A));
94: PetscCall(KSPSetFromOptions(ksp));
95: PetscCall(KSPSolve(ksp, b, x));
96: /* PetscCall(KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD)); */
98: /* Check error */
99: PetscCall(VecAXPY(x, none, u));
100: PetscCall(VecNorm(x, NORM_2, &norm));
101: PetscCall(KSPGetIterationNumber(ksp, &its));
102: if (norm >= 1.e-12) {
103: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
104: } else {
105: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error < 1.e-12, Iterations %" PetscInt_FMT "\n", its));
106: }
108: /* Free work space */
109: PetscCall(VecDestroy(&x));
110: PetscCall(VecDestroy(&u));
111: PetscCall(VecDestroy(&b));
112: PetscCall(MatDestroy(&A));
113: if (use_random) PetscCall(PetscRandomDestroy(&rctx));
114: PetscCall(KSPDestroy(&ksp));
115: PetscCall(PetscFinalize());
116: return 0;
117: }
119: PetscErrorCode FormTestMatrix(Mat A, PetscInt n, TestType type)
120: {
121: PetscScalar val[5];
122: PetscInt i, j, Ii, J, col[5], Istart, Iend;
124: PetscFunctionBeginUser;
125: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
126: if (type == TEST_1) {
127: val[0] = 1.0;
128: val[1] = 4.0;
129: val[2] = -2.0;
130: for (i = 1; i < n - 1; i++) {
131: col[0] = i - 1;
132: col[1] = i;
133: col[2] = i + 1;
134: PetscCall(MatSetValues(A, 1, &i, 3, col, val, INSERT_VALUES));
135: }
136: i = n - 1;
137: col[0] = n - 2;
138: col[1] = n - 1;
139: PetscCall(MatSetValues(A, 1, &i, 2, col, val, INSERT_VALUES));
140: i = 0;
141: col[0] = 0;
142: col[1] = 1;
143: val[0] = 4.0;
144: val[1] = -2.0;
145: PetscCall(MatSetValues(A, 1, &i, 2, col, val, INSERT_VALUES));
146: } else if (type == TEST_2) {
147: val[0] = 1.0;
148: val[1] = 0.0;
149: val[2] = 2.0;
150: val[3] = 1.0;
151: for (i = 2; i < n - 1; i++) {
152: col[0] = i - 2;
153: col[1] = i - 1;
154: col[2] = i;
155: col[3] = i + 1;
156: PetscCall(MatSetValues(A, 1, &i, 4, col, val, INSERT_VALUES));
157: }
158: i = n - 1;
159: col[0] = n - 3;
160: col[1] = n - 2;
161: col[2] = n - 1;
162: PetscCall(MatSetValues(A, 1, &i, 3, col, val, INSERT_VALUES));
163: i = 1;
164: col[0] = 0;
165: col[1] = 1;
166: col[2] = 2;
167: PetscCall(MatSetValues(A, 1, &i, 3, col, &val[1], INSERT_VALUES));
168: i = 0;
169: PetscCall(MatSetValues(A, 1, &i, 2, col, &val[2], INSERT_VALUES));
170: } else if (type == TEST_3) {
171: val[0] = PETSC_i * 2.0;
172: val[1] = 4.0;
173: val[2] = 0.0;
174: val[3] = 1.0;
175: val[4] = 0.7;
176: for (i = 1; i < n - 3; i++) {
177: col[0] = i - 1;
178: col[1] = i;
179: col[2] = i + 1;
180: col[3] = i + 2;
181: col[4] = i + 3;
182: PetscCall(MatSetValues(A, 1, &i, 5, col, val, INSERT_VALUES));
183: }
184: i = n - 3;
185: col[0] = n - 4;
186: col[1] = n - 3;
187: col[2] = n - 2;
188: col[3] = n - 1;
189: PetscCall(MatSetValues(A, 1, &i, 4, col, val, INSERT_VALUES));
190: i = n - 2;
191: col[0] = n - 3;
192: col[1] = n - 2;
193: col[2] = n - 1;
194: PetscCall(MatSetValues(A, 1, &i, 3, col, val, INSERT_VALUES));
195: i = n - 1;
196: col[0] = n - 2;
197: col[1] = n - 1;
198: PetscCall(MatSetValues(A, 1, &i, 2, col, val, INSERT_VALUES));
199: i = 0;
200: col[0] = 0;
201: col[1] = 1;
202: col[2] = 2;
203: col[3] = 3;
204: PetscCall(MatSetValues(A, 1, &i, 4, col, &val[1], INSERT_VALUES));
205: } else if (type == HELMHOLTZ_1) {
206: /* Problem domain: unit square: (0,1) x (0,1)
207: Solve Helmholtz equation:
208: -delta u - sigma1*u + i*sigma2*u = f,
209: where delta = Laplace operator
210: Dirichlet b.c.'s on all sides
211: */
212: PetscRandom rctx;
213: PetscReal h2, sigma1 = 5.0;
214: PetscScalar sigma2;
215: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL));
216: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
217: PetscCall(PetscRandomSetFromOptions(rctx));
218: PetscCall(PetscRandomSetInterval(rctx, 0.0, PETSC_i));
219: h2 = 1.0 / ((n + 1) * (n + 1));
220: for (Ii = Istart; Ii < Iend; Ii++) {
221: *val = -1.0;
222: i = Ii / n;
223: j = Ii - i * n;
224: if (i > 0) {
225: J = Ii - n;
226: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
227: }
228: if (i < n - 1) {
229: J = Ii + n;
230: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
231: }
232: if (j > 0) {
233: J = Ii - 1;
234: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
235: }
236: if (j < n - 1) {
237: J = Ii + 1;
238: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
239: }
240: PetscCall(PetscRandomGetValue(rctx, &sigma2));
241: *val = 4.0 - sigma1 * h2 + sigma2 * h2;
242: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, val, ADD_VALUES));
243: }
244: PetscCall(PetscRandomDestroy(&rctx));
245: } else if (type == HELMHOLTZ_2) {
246: /* Problem domain: unit square: (0,1) x (0,1)
247: Solve Helmholtz equation:
248: -delta u - sigma1*u = f,
249: where delta = Laplace operator
250: Dirichlet b.c.'s on 3 sides
251: du/dn = i*alpha*u on (1,y), 0<y<1
252: */
253: PetscReal h2, sigma1 = 200.0;
254: PetscScalar alpha_h;
255: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL));
256: h2 = 1.0 / ((n + 1) * (n + 1));
257: alpha_h = (PETSC_i * 10.0) / (PetscReal)(n + 1); /* alpha_h = alpha * h */
258: for (Ii = Istart; Ii < Iend; Ii++) {
259: *val = -1.0;
260: i = Ii / n;
261: j = Ii - i * n;
262: if (i > 0) {
263: J = Ii - n;
264: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
265: }
266: if (i < n - 1) {
267: J = Ii + n;
268: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
269: }
270: if (j > 0) {
271: J = Ii - 1;
272: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
273: }
274: if (j < n - 1) {
275: J = Ii + 1;
276: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
277: }
278: *val = 4.0 - sigma1 * h2;
279: if (!((Ii + 1) % n)) *val += alpha_h;
280: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, val, ADD_VALUES));
281: }
282: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_USER_INPUT, "FormTestMatrix: unknown test matrix type");
284: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
285: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*TEST
291: build:
292: requires: complex
294: test:
295: args: -ksp_gmres_cgs_refinement_type refine_always -n 6 -ksp_monitor_short -p 5 -norandom -ksp_type gmres -pc_type jacobi -ksp_max_it 15
296: requires: complex
298: test:
299: suffix: 2
300: nsize: 3
301: requires: complex
302: args: -ksp_gmres_cgs_refinement_type refine_always -n 6 -ksp_monitor_short -p 5 -norandom -ksp_type gmres -pc_type jacobi -ksp_max_it 15
303: output_file: output/ex17_1.out
305: test:
306: suffix: superlu_dist
307: requires: superlu_dist complex
308: args: -n 6 -p 5 -norandom -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_colperm MMD_ATA
310: test:
311: suffix: superlu_dist_2
312: requires: superlu_dist complex
313: nsize: 3
314: output_file: output/ex17_superlu_dist.out
315: args: -n 6 -p 5 -norandom -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_colperm MMD_ATA
317: TEST*/