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, NULL, 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: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
104: /* Free work space */
105: PetscCall(VecDestroy(&x));
106: PetscCall(VecDestroy(&u));
107: PetscCall(VecDestroy(&b));
108: PetscCall(MatDestroy(&A));
109: if (use_random) PetscCall(PetscRandomDestroy(&rctx));
110: PetscCall(KSPDestroy(&ksp));
111: PetscCall(PetscFinalize());
112: return 0;
113: }
115: PetscErrorCode FormTestMatrix(Mat A, PetscInt n, TestType type)
116: {
117: PetscScalar val[5];
118: PetscInt i, j, Ii, J, col[5], Istart, Iend;
120: PetscFunctionBeginUser;
121: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
122: if (type == TEST_1) {
123: val[0] = 1.0;
124: val[1] = 4.0;
125: val[2] = -2.0;
126: for (i = 1; i < n - 1; i++) {
127: col[0] = i - 1;
128: col[1] = i;
129: col[2] = i + 1;
130: PetscCall(MatSetValues(A, 1, &i, 3, col, val, INSERT_VALUES));
131: }
132: i = n - 1;
133: col[0] = n - 2;
134: col[1] = n - 1;
135: PetscCall(MatSetValues(A, 1, &i, 2, col, val, INSERT_VALUES));
136: i = 0;
137: col[0] = 0;
138: col[1] = 1;
139: val[0] = 4.0;
140: val[1] = -2.0;
141: PetscCall(MatSetValues(A, 1, &i, 2, col, val, INSERT_VALUES));
142: } else if (type == TEST_2) {
143: val[0] = 1.0;
144: val[1] = 0.0;
145: val[2] = 2.0;
146: val[3] = 1.0;
147: for (i = 2; i < n - 1; i++) {
148: col[0] = i - 2;
149: col[1] = i - 1;
150: col[2] = i;
151: col[3] = i + 1;
152: PetscCall(MatSetValues(A, 1, &i, 4, col, val, INSERT_VALUES));
153: }
154: i = n - 1;
155: col[0] = n - 3;
156: col[1] = n - 2;
157: col[2] = n - 1;
158: PetscCall(MatSetValues(A, 1, &i, 3, col, val, INSERT_VALUES));
159: i = 1;
160: col[0] = 0;
161: col[1] = 1;
162: col[2] = 2;
163: PetscCall(MatSetValues(A, 1, &i, 3, col, &val[1], INSERT_VALUES));
164: i = 0;
165: PetscCall(MatSetValues(A, 1, &i, 2, col, &val[2], INSERT_VALUES));
166: } else if (type == TEST_3) {
167: val[0] = PETSC_i * 2.0;
168: val[1] = 4.0;
169: val[2] = 0.0;
170: val[3] = 1.0;
171: val[4] = 0.7;
172: for (i = 1; i < n - 3; i++) {
173: col[0] = i - 1;
174: col[1] = i;
175: col[2] = i + 1;
176: col[3] = i + 2;
177: col[4] = i + 3;
178: PetscCall(MatSetValues(A, 1, &i, 5, col, val, INSERT_VALUES));
179: }
180: i = n - 3;
181: col[0] = n - 4;
182: col[1] = n - 3;
183: col[2] = n - 2;
184: col[3] = n - 1;
185: PetscCall(MatSetValues(A, 1, &i, 4, col, val, INSERT_VALUES));
186: i = n - 2;
187: col[0] = n - 3;
188: col[1] = n - 2;
189: col[2] = n - 1;
190: PetscCall(MatSetValues(A, 1, &i, 3, col, val, INSERT_VALUES));
191: i = n - 1;
192: col[0] = n - 2;
193: col[1] = n - 1;
194: PetscCall(MatSetValues(A, 1, &i, 2, col, val, INSERT_VALUES));
195: i = 0;
196: col[0] = 0;
197: col[1] = 1;
198: col[2] = 2;
199: col[3] = 3;
200: PetscCall(MatSetValues(A, 1, &i, 4, col, &val[1], INSERT_VALUES));
201: } else if (type == HELMHOLTZ_1) {
202: /* Problem domain: unit square: (0,1) x (0,1)
203: Solve Helmholtz equation:
204: -delta u - sigma1*u + i*sigma2*u = f,
205: where delta = Laplace operator
206: Dirichlet b.c.'s on all sides
207: */
208: PetscRandom rctx;
209: PetscReal h2, sigma1 = 5.0;
210: PetscScalar sigma2;
211: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL));
212: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
213: PetscCall(PetscRandomSetFromOptions(rctx));
214: PetscCall(PetscRandomSetInterval(rctx, 0.0, PETSC_i));
215: h2 = 1.0 / ((n + 1) * (n + 1));
216: for (Ii = Istart; Ii < Iend; Ii++) {
217: *val = -1.0;
218: i = Ii / n;
219: j = Ii - i * n;
220: if (i > 0) {
221: J = Ii - n;
222: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
223: }
224: if (i < n - 1) {
225: J = Ii + n;
226: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
227: }
228: if (j > 0) {
229: J = Ii - 1;
230: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
231: }
232: if (j < n - 1) {
233: J = Ii + 1;
234: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
235: }
236: PetscCall(PetscRandomGetValue(rctx, &sigma2));
237: *val = 4.0 - sigma1 * h2 + sigma2 * h2;
238: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, val, ADD_VALUES));
239: }
240: PetscCall(PetscRandomDestroy(&rctx));
241: } else if (type == HELMHOLTZ_2) {
242: /* Problem domain: unit square: (0,1) x (0,1)
243: Solve Helmholtz equation:
244: -delta u - sigma1*u = f,
245: where delta = Laplace operator
246: Dirichlet b.c.'s on 3 sides
247: du/dn = i*alpha*u on (1,y), 0<y<1
248: */
249: PetscReal h2, sigma1 = 200.0;
250: PetscScalar alpha_h;
251: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL));
252: h2 = 1.0 / ((n + 1) * (n + 1));
253: alpha_h = (PETSC_i * 10.0) / (PetscReal)(n + 1); /* alpha_h = alpha * h */
254: for (Ii = Istart; Ii < Iend; Ii++) {
255: *val = -1.0;
256: i = Ii / n;
257: j = Ii - i * n;
258: if (i > 0) {
259: J = Ii - n;
260: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
261: }
262: if (i < n - 1) {
263: J = Ii + n;
264: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
265: }
266: if (j > 0) {
267: J = Ii - 1;
268: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
269: }
270: if (j < n - 1) {
271: J = Ii + 1;
272: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES));
273: }
274: *val = 4.0 - sigma1 * h2;
275: if (!((Ii + 1) % n)) *val += alpha_h;
276: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, val, ADD_VALUES));
277: }
278: } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_USER_INPUT, "FormTestMatrix: unknown test matrix type");
280: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
281: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /*TEST
287: build:
288: requires: complex
290: test:
291: args: -ksp_orthogonalization_cgs_refinement_type refine_always -n 6 -ksp_monitor -p 5 -norandom -ksp_type gmres -pc_type jacobi -ksp_max_it 15
292: requires: complex
294: test:
295: suffix: 2
296: nsize: 3
297: requires: complex
298: args: -ksp_orthogonalization_cgs_refinement_type refine_always -n 6 -ksp_monitor -p 5 -norandom -ksp_type gmres -pc_type jacobi -ksp_max_it 15
299: output_file: output/ex17_1.out
301: test:
302: suffix: superlu_dist
303: requires: superlu_dist complex
304: args: -n 6 -p 5 -norandom -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_colperm MMD_ATA
306: test:
307: suffix: superlu_dist_2
308: requires: superlu_dist complex
309: nsize: 3
310: output_file: output/ex17_superlu_dist.out
311: args: -n 6 -p 5 -norandom -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_colperm MMD_ATA
313: TEST*/