Actual source code: ex60.c
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;
16: PetscInitialize(&argc, &args, NULL,help);
17: PetscOptionsGetInt(NULL,NULL, "-M", &M, NULL);
18: PetscOptionsGetInt(NULL,NULL, "-N", &N, NULL);
19: PetscOptionsGetBool(NULL,NULL, "-user_subdomains", &userSubdomains, NULL);
20: /* Do parallel decomposition */
21: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
22: MPI_Comm_size(PETSC_COMM_WORLD, &size);
23: sized = (PetscMPIInt) PetscSqrtReal((PetscReal) size);
27: /* Assemble the matrix for the five point stencil, YET AGAIN
28: Every other process will be empty */
29: MatCreate(PETSC_COMM_WORLD, &A);
30: m = (sized > 1) ? (rank % 2) ? 0 : 2*M/sized : M;
31: n = N/sized;
32: MatSetSizes(A, m*n, m*n, M*N, M*N);
33: MatSetFromOptions(A);
34: MatSetUp(A);
35: MatGetOwnershipRange(A, &rstart, &rend);
36: for (r = rstart; r < rend; ++r) {
37: const PetscScalar diag = 4.0, offdiag = -1.0;
38: const PetscInt i = r/N;
39: const PetscInt j = r - i*N;
40: PetscInt c;
42: if (i > 0) {c = r - n; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
43: if (i < M-1) {c = r + n; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
44: if (j > 0) {c = r - 1; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
45: if (j < N-1) {c = r + 1; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
46: MatSetValues(A, 1, &r, 1, &r, &diag, INSERT_VALUES);
47: }
48: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
49: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
50: /* Setup Solve */
51: MatCreateVecs(A, &x, &b);
52: VecDuplicate(x, &u);
53: VecSet(u, 1.0);
54: MatMult(A, u, b);
55: KSPCreate(PETSC_COMM_WORLD, &ksp);
56: KSPSetOperators(ksp, A, A);
57: KSPGetPC(ksp, &pc);
58: PCSetType(pc, PCASM);
59: /* Setup ASM by hand */
60: if (userSubdomains) {
61: IS is;
62: PetscInt *rows;
64: /* Use no overlap for now */
65: PetscMalloc1(rend-rstart, &rows);
66: for (r = rstart; r < rend; ++r) rows[r-rstart] = r;
67: ISCreateGeneral(PETSC_COMM_SELF, rend-rstart, rows, PETSC_OWN_POINTER, &is);
68: PCASMSetLocalSubdomains(pc, 1, &is, &is);
69: ISDestroy(&is);
70: }
71: KSPSetFromOptions(ksp);
72: /* Solve and Compare */
73: KSPSolve(ksp, b, x);
74: VecAXPY(x, -1.0, u);
75: VecNorm(x, NORM_INFINITY, &error);
76: PetscPrintf(PETSC_COMM_WORLD, "Infinity norm of the error: %g\n", (double) error);
77: /* Cleanup */
78: KSPDestroy(&ksp);
79: MatDestroy(&A);
80: VecDestroy(&u);
81: VecDestroy(&x);
82: VecDestroy(&b);
83: PetscFinalize();
84: return 0;
85: }
87: /*TEST
89: test:
90: suffix: 0
91: args: -ksp_view
93: test:
94: requires: !sycl kokkos_kernels
95: suffix: 0_kokkos
96: args: -ksp_view -mat_type aijkokkos
98: test:
99: requires: cuda
100: suffix: 0_cuda
101: args: -ksp_view -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse
103: test:
104: suffix: 1
105: nsize: 4
106: args: -ksp_view
108: test:
109: requires: !sycl kokkos_kernels
110: suffix: 1_kokkos
111: nsize: 4
112: args: -ksp_view -mat_type aijkokkos
114: test:
115: requires: cuda
116: suffix: 1_cuda
117: nsize: 4
118: args: -ksp_view -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse
120: test:
121: suffix: 2
122: nsize: 4
123: args: -user_subdomains -ksp_view
125: test:
126: requires: !sycl kokkos_kernels
127: suffix: 2_kokkos
128: nsize: 4
129: args: -user_subdomains -ksp_view -mat_type aijkokkos
131: test:
132: requires: cuda
133: suffix: 2_cuda
134: nsize: 4
135: args: -user_subdomains -ksp_view -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse
137: TEST*/