Actual source code: ex17.c
petsc-3.12.0 2019-09-29
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 {TEST_1,TEST_2,TEST_3,HELMHOLTZ_1,HELMHOLTZ_2} TestType;
8: extern PetscErrorCode FormTestMatrix(Mat,PetscInt,TestType);
10: int main(int argc,char **args)
11: {
12: Vec x,b,u; /* approx solution, RHS, exact solution */
13: Mat A; /* linear system matrix */
14: KSP ksp; /* KSP context */
16: PetscInt n = 10,its, dim,p = 1,use_random;
17: PetscScalar none = -1.0,pfive = 0.5;
18: PetscReal norm;
19: PetscRandom rctx;
20: TestType type;
21: PetscBool flg;
23: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
24: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
25: PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
26: switch (p) {
27: case 1: type = TEST_1; dim = n; break;
28: case 2: type = TEST_2; dim = n; break;
29: case 3: type = TEST_3; dim = n; break;
30: case 4: type = HELMHOLTZ_1; dim = n*n; break;
31: case 5: type = HELMHOLTZ_2; dim = n*n; break;
32: default: type = TEST_1; dim = n;
33: }
35: /* Create vectors */
36: VecCreate(PETSC_COMM_WORLD,&x);
37: VecSetSizes(x,PETSC_DECIDE,dim);
38: VecSetFromOptions(x);
39: VecDuplicate(x,&b);
40: VecDuplicate(x,&u);
42: use_random = 1;
43: flg = PETSC_FALSE;
44: PetscOptionsGetBool(NULL,NULL,"-norandom",&flg,NULL);
45: if (flg) {
46: use_random = 0;
47: VecSet(u,pfive);
48: } else {
49: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
50: PetscRandomSetFromOptions(rctx);
51: VecSetRandom(u,rctx);
52: }
54: /* Create and assemble matrix */
55: MatCreate(PETSC_COMM_WORLD,&A);
56: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
57: MatSetFromOptions(A);
58: MatSetUp(A);
59: FormTestMatrix(A,n,type);
60: MatMult(A,u,b);
61: flg = PETSC_FALSE;
62: PetscOptionsGetBool(NULL,NULL,"-printout",&flg,NULL);
63: if (flg) {
64: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
65: VecView(u,PETSC_VIEWER_STDOUT_WORLD);
66: VecView(b,PETSC_VIEWER_STDOUT_WORLD);
67: }
69: /* Create KSP context; set operators and options; solve linear system */
70: KSPCreate(PETSC_COMM_WORLD,&ksp);
71: KSPSetOperators(ksp,A,A);
72: KSPSetFromOptions(ksp);
73: KSPSolve(ksp,b,x);
74: /* KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD); */
76: /* Check error */
77: VecAXPY(x,none,u);
78: VecNorm(x,NORM_2,&norm);
79: KSPGetIterationNumber(ksp,&its);
80: if (norm >= 1.e-12) {
81: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
82: } else {
83: PetscPrintf(PETSC_COMM_WORLD,"Norm of error < 1.e-12, Iterations %D\n",its);
84: }
86: /* Free work space */
87: VecDestroy(&x); VecDestroy(&u);
88: VecDestroy(&b); MatDestroy(&A);
89: if (use_random) {PetscRandomDestroy(&rctx);}
90: KSPDestroy(&ksp);
91: PetscFinalize();
92: return ierr;
93: }
95: PetscErrorCode FormTestMatrix(Mat A,PetscInt n,TestType type)
96: {
97: PetscScalar val[5];
99: PetscInt i,j,Ii,J,col[5],Istart,Iend;
101: MatGetOwnershipRange(A,&Istart,&Iend);
102: if (type == TEST_1) {
103: val[0] = 1.0; val[1] = 4.0; val[2] = -2.0;
104: for (i=1; i<n-1; i++) {
105: col[0] = i-1; col[1] = i; col[2] = i+1;
106: MatSetValues(A,1,&i,3,col,val,INSERT_VALUES);
107: }
108: i = n-1; col[0] = n-2; col[1] = n-1;
109: MatSetValues(A,1,&i,2,col,val,INSERT_VALUES);
110: i = 0; col[0] = 0; col[1] = 1; val[0] = 4.0; val[1] = -2.0;
111: MatSetValues(A,1,&i,2,col,val,INSERT_VALUES);
112: } else if (type == TEST_2) {
113: val[0] = 1.0; val[1] = 0.0; val[2] = 2.0; val[3] = 1.0;
114: for (i=2; i<n-1; i++) {
115: col[0] = i-2; col[1] = i-1; col[2] = i; col[3] = i+1;
116: MatSetValues(A,1,&i,4,col,val,INSERT_VALUES);
117: }
118: i = n-1; col[0] = n-3; col[1] = n-2; col[2] = n-1;
119: MatSetValues(A,1,&i,3,col,val,INSERT_VALUES);
120: i = 1; col[0] = 0; col[1] = 1; col[2] = 2;
121: MatSetValues(A,1,&i,3,col,&val[1],INSERT_VALUES);
122: i = 0;
123: MatSetValues(A,1,&i,2,col,&val[2],INSERT_VALUES);
124: } else if (type == TEST_3) {
125: val[0] = PETSC_i * 2.0;
126: val[1] = 4.0; val[2] = 0.0; val[3] = 1.0; val[4] = 0.7;
127: for (i=1; i<n-3; i++) {
128: col[0] = i-1; col[1] = i; col[2] = i+1; col[3] = i+2; col[4] = i+3;
129: MatSetValues(A,1,&i,5,col,val,INSERT_VALUES);
130: }
131: i = n-3; col[0] = n-4; col[1] = n-3; col[2] = n-2; col[3] = n-1;
132: MatSetValues(A,1,&i,4,col,val,INSERT_VALUES);
133: i = n-2; col[0] = n-3; col[1] = n-2; col[2] = n-1;
134: MatSetValues(A,1,&i,3,col,val,INSERT_VALUES);
135: i = n-1; col[0] = n-2; col[1] = n-1;
136: MatSetValues(A,1,&i,2,col,val,INSERT_VALUES);
137: i = 0; col[0] = 0; col[1] = 1; col[2] = 2; col[3] = 3;
138: MatSetValues(A,1,&i,4,col,&val[1],INSERT_VALUES);
139: } else if (type == HELMHOLTZ_1) {
140: /* Problem domain: unit square: (0,1) x (0,1)
141: Solve Helmholtz equation:
142: -delta u - sigma1*u + i*sigma2*u = f,
143: where delta = Laplace operator
144: Dirichlet b.c.'s on all sides
145: */
146: PetscRandom rctx;
147: PetscReal h2,sigma1 = 5.0;
148: PetscScalar sigma2;
149: PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);
150: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
151: PetscRandomSetFromOptions(rctx);
152: PetscRandomSetInterval(rctx,0.0,PETSC_i);
153: h2 = 1.0/((n+1)*(n+1));
154: for (Ii=Istart; Ii<Iend; Ii++) {
155: *val = -1.0; i = Ii/n; j = Ii - i*n;
156: if (i>0) {
157: J = Ii-n; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
158: }
159: if (i<n-1) {
160: J = Ii+n; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
161: }
162: if (j>0) {
163: J = Ii-1; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
164: }
165: if (j<n-1) {
166: J = Ii+1; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
167: }
168: PetscRandomGetValue(rctx,&sigma2);
169: *val = 4.0 - sigma1*h2 + sigma2*h2;
170: MatSetValues(A,1,&Ii,1,&Ii,val,ADD_VALUES);
171: }
172: PetscRandomDestroy(&rctx);
173: } else if (type == HELMHOLTZ_2) {
174: /* Problem domain: unit square: (0,1) x (0,1)
175: Solve Helmholtz equation:
176: -delta u - sigma1*u = f,
177: where delta = Laplace operator
178: Dirichlet b.c.'s on 3 sides
179: du/dn = i*alpha*u on (1,y), 0<y<1
180: */
181: PetscReal h2,sigma1 = 200.0;
182: PetscScalar alpha_h;
183: PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);
184: h2 = 1.0/((n+1)*(n+1));
185: alpha_h = (PETSC_i * 10.0) / (PetscReal)(n+1); /* alpha_h = alpha * h */
186: for (Ii=Istart; Ii<Iend; Ii++) {
187: *val = -1.0; i = Ii/n; j = Ii - i*n;
188: if (i>0) {
189: J = Ii-n; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
190: }
191: if (i<n-1) {
192: J = Ii+n; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
193: }
194: if (j>0) {
195: J = Ii-1; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
196: }
197: if (j<n-1) {
198: J = Ii+1; MatSetValues(A,1,&Ii,1,&J,val,ADD_VALUES);
199: }
200: *val = 4.0 - sigma1*h2;
201: if (!((Ii+1)%n)) *val += alpha_h;
202: MatSetValues(A,1,&Ii,1,&Ii,val,ADD_VALUES);
203: }
204: } else SETERRQ(PetscObjectComm((PetscObject)A),1,"FormTestMatrix: unknown test matrix type");
206: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
207: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
209: return 0;
210: }
212: /*TEST
214: build:
215: requires: complex
217: test:
218: 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
219: requires: complex
221: test:
222: suffix: 2
223: nsize: 3
224: requires: complex
225: 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
226: output_file: output/ex17_1.out
228: test:
229: suffix: superlu_dist
230: requires: superlu_dist complex
231: args: -n 6 -p 5 -norandom -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_colperm MMD_ATA
233: test:
234: suffix: superlu_dist_2
235: requires: superlu_dist complex
236: nsize: 3
237: output_file: output/ex17_superlu_dist.out
238: args: -n 6 -p 5 -norandom -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_colperm MMD_ATA
240: TEST*/