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