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;
29: PetscInitialize(&argc, &args, (char *)0, help);
30: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
31: 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: VecCreate(PETSC_COMM_WORLD, &x);
60: VecSetSizes(x, PETSC_DECIDE, dim);
61: VecSetFromOptions(x);
62: VecDuplicate(x, &b);
63: VecDuplicate(x, &u);
65: use_random = 1;
66: flg = PETSC_FALSE;
67: PetscOptionsGetBool(NULL, NULL, "-norandom", &flg, NULL);
68: if (flg) {
69: use_random = 0;
70: VecSet(u, pfive);
71: } else {
72: PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
73: PetscRandomSetFromOptions(rctx);
74: VecSetRandom(u, rctx);
75: }
77: /* Create and assemble matrix */
78: MatCreate(PETSC_COMM_WORLD, &A);
79: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim);
80: MatSetFromOptions(A);
81: MatSetUp(A);
82: FormTestMatrix(A, n, type);
83: MatMult(A, u, b);
84: flg = PETSC_FALSE;
85: PetscOptionsGetBool(NULL, NULL, "-printout", &flg, NULL);
86: if (flg) {
87: MatView(A, PETSC_VIEWER_STDOUT_WORLD);
88: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
89: VecView(b, PETSC_VIEWER_STDOUT_WORLD);
90: }
92: /* Create KSP context; set operators and options; solve linear system */
93: KSPCreate(PETSC_COMM_WORLD, &ksp);
94: KSPSetOperators(ksp, A, A);
95: KSPSetFromOptions(ksp);
96: KSPSolve(ksp, b, x);
97: /* KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD); */
99: /* Check error */
100: VecAXPY(x, none, u);
101: VecNorm(x, NORM_2, &norm);
102: KSPGetIterationNumber(ksp, &its);
103: if (norm >= 1.e-12) {
104: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
105: } else {
106: PetscPrintf(PETSC_COMM_WORLD, "Norm of error < 1.e-12, Iterations %" PetscInt_FMT "\n", its);
107: }
109: /* Free work space */
110: VecDestroy(&x);
111: VecDestroy(&u);
112: VecDestroy(&b);
113: MatDestroy(&A);
114: if (use_random) PetscRandomDestroy(&rctx);
115: KSPDestroy(&ksp);
116: 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: 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: 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: 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: 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: 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: 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: MatSetValues(A, 1, &i, 3, col, &val[1], INSERT_VALUES);
168: i = 0;
169: 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: 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: 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: MatSetValues(A, 1, &i, 3, col, val, INSERT_VALUES);
195: i = n - 1;
196: col[0] = n - 2;
197: col[1] = n - 1;
198: 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: 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: PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL);
216: PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
217: PetscRandomSetFromOptions(rctx);
218: 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: MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES);
227: }
228: if (i < n - 1) {
229: J = Ii + n;
230: MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES);
231: }
232: if (j > 0) {
233: J = Ii - 1;
234: MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES);
235: }
236: if (j < n - 1) {
237: J = Ii + 1;
238: MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES);
239: }
240: PetscRandomGetValue(rctx, &sigma2);
241: *val = 4.0 - sigma1 * h2 + sigma2 * h2;
242: MatSetValues(A, 1, &Ii, 1, &Ii, val, ADD_VALUES);
243: }
244: 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: 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: MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES);
265: }
266: if (i < n - 1) {
267: J = Ii + n;
268: MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES);
269: }
270: if (j > 0) {
271: J = Ii - 1;
272: MatSetValues(A, 1, &Ii, 1, &J, val, ADD_VALUES);
273: }
274: if (j < n - 1) {
275: J = Ii + 1;
276: 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: 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: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
285: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
287: return 0;
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*/