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