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:   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*/