Actual source code: ex6.c
petsc-3.11.4 2019-09-28
1: static char help[] = "Working out corner cases of the ASM preconditioner.\n\n";
3: #include <petscksp.h>
5: int main(int argc,char **args)
6: {
7: KSP ksp;
8: PC pc;
9: Mat A;
10: Vec u, x, b;
11: PetscReal error;
12: PetscMPIInt rank, size, sized;
13: PetscInt M = 8, N = 8, m, n, rstart, rend, r;
14: PetscBool userSubdomains = PETSC_FALSE;
17: PetscInitialize(&argc, &args, NULL,help);if (ierr) return ierr;
18: PetscOptionsGetInt(NULL,NULL, "-M", &M, NULL);
19: PetscOptionsGetInt(NULL,NULL, "-N", &N, NULL);
20: PetscOptionsGetBool(NULL,NULL, "-user_subdomains", &userSubdomains, NULL);
21: /* Do parallel decomposition */
22: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
23: MPI_Comm_size(PETSC_COMM_WORLD, &size);
24: sized = (PetscMPIInt) PetscSqrtReal((PetscReal) size);
25: if (PetscSqr(sized) != size) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "This test may only be run on a nubmer of processes which is a perfect square, not %d", (int) size);
26: if (M % sized) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "The number of x-vertices %D does not divide the number of x-processes %d", M, (int) sized);
27: if (N % sized) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "The number of y-vertices %D does not divide the number of y-processes %d", N, (int) sized);
28: /* Assemble the matrix for the five point stencil, YET AGAIN
29: Every other process will be empty */
30: MatCreate(PETSC_COMM_WORLD, &A);
31: m = (sized > 1) ? (rank % 2) ? 0 : 2*M/sized : M;
32: n = N/sized;
33: MatSetSizes(A, m*n, m*n, M*N, M*N);
34: MatSetFromOptions(A);
35: MatSetUp(A);
36: MatGetOwnershipRange(A, &rstart, &rend);
37: for (r = rstart; r < rend; ++r) {
38: const PetscScalar diag = 4.0, offdiag = -1.0;
39: const PetscInt i = r/N;
40: const PetscInt j = r - i*N;
41: PetscInt c;
43: if (i > 0) {c = r - n; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
44: if (i < M-1) {c = r + n; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
45: if (j > 0) {c = r - 1; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
46: if (j < N-1) {c = r + 1; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
47: MatSetValues(A, 1, &r, 1, &r, &diag, INSERT_VALUES);
48: }
49: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
50: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
51: /* Setup Solve */
52: VecCreate(PETSC_COMM_WORLD, &b);
53: VecSetSizes(b, m*n, PETSC_DETERMINE);
54: VecSetFromOptions(b);
55: VecDuplicate(b, &u);
56: VecDuplicate(b, &x);
57: VecSet(u, 1.0);
58: MatMult(A, u, b);
59: KSPCreate(PETSC_COMM_WORLD, &ksp);
60: KSPSetOperators(ksp, A, A);
61: KSPGetPC(ksp, &pc);
62: PCSetType(pc, PCASM);
63: /* Setup ASM by hand */
64: if (userSubdomains) {
65: IS is;
66: PetscInt *rows;
68: /* Use no overlap for now */
69: PetscMalloc1(rend-rstart, &rows);
70: for (r = rstart; r < rend; ++r) rows[r-rstart] = r;
71: ISCreateGeneral(PETSC_COMM_SELF, rend-rstart, rows, PETSC_OWN_POINTER, &is);
72: PCASMSetLocalSubdomains(pc, 1, &is, &is);
73: ISDestroy(&is);
74: }
75: KSPSetFromOptions(ksp);
76: /* Solve and Compare */
77: KSPSolve(ksp, b, x);
78: VecAXPY(x, -1.0, u);
79: VecNorm(x, NORM_INFINITY, &error);
80: PetscPrintf(PETSC_COMM_WORLD, "Infinity norm of the error: %g\n", (double) error);
81: /* Cleanup */
82: KSPDestroy(&ksp);
83: MatDestroy(&A);
84: VecDestroy(&u);
85: VecDestroy(&x);
86: VecDestroy(&b);
87: PetscFinalize();
88: return ierr;
89: }
92: /*TEST
94: test:
95: suffix: 0
96: args: -ksp_view
98: test:
99: suffix: 1
100: nsize: 4
101: args: -ksp_view
103: test:
104: suffix: 2
105: nsize: 4
106: args: -user_subdomains -ksp_view ${ARGS}
108: TEST*/