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