Actual source code: ex53.c
2: static char help[] = "Tests setup PCFIELDSPLIT with blocked IS.\n\n";
3: /*
4: Contributed by Hoang Giang Bui, June 2017.
5: */
6: #include <petscksp.h>
8: int main(int argc, char *argv[])
9: {
10: Mat A;
11: KSP ksp;
12: PC pc;
13: PetscInt Istart, Iend, local_m, local_n, i;
14: PetscMPIInt rank;
15: PetscInt method = 2, mat_size = 40, block_size = 2, *A_indices = NULL, *B_indices = NULL, A_size = 0, B_size = 0;
16: IS A_IS, B_IS;
19: PetscInitialize(&argc, &argv, (char *)0, help);
20: MPI_Comm_rank(MPI_COMM_WORLD, &rank);
22: PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-mat_size", &mat_size, PETSC_NULL);
23: PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-method", &method, PETSC_NULL);
24: PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-block_size", &block_size, PETSC_NULL);
26: if (rank == 0) PetscPrintf(PETSC_COMM_SELF, " matrix size = %" PetscInt_FMT ", block size = %" PetscInt_FMT "\n", mat_size, block_size);
28: MatCreate(PETSC_COMM_WORLD, &A);
29: MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, mat_size, mat_size);
30: MatSetType(A, MATMPIAIJ);
31: MatSetUp(A);
33: MatGetOwnershipRange(A, &Istart, &Iend);
35: for (i = Istart; i < Iend; ++i) {
36: MatSetValue(A, i, i, 2, INSERT_VALUES);
37: if (i < mat_size - 1) MatSetValue(A, i, i + 1, -1, INSERT_VALUES);
38: if (i > 0) MatSetValue(A, i, i - 1, -1, INSERT_VALUES);
39: }
41: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
42: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
44: MatGetLocalSize(A, &local_m, &local_n);
46: /* Create Index Sets */
47: if (rank == 0) {
48: if (method > 1) {
49: /* with method > 1, the fieldsplit B is set to zero */
50: A_size = Iend - Istart;
51: B_size = 0;
52: } else {
53: /* with method = 1, the fieldsplit A and B is equal. It is noticed that A_size, or N/4, must be divided by block_size */
54: A_size = (Iend - Istart) / 2;
55: B_size = (Iend - Istart) / 2;
56: }
57: PetscCalloc1(A_size, &A_indices);
58: PetscCalloc1(B_size, &B_indices);
59: for (i = 0; i < A_size; ++i) A_indices[i] = Istart + i;
60: for (i = 0; i < B_size; ++i) B_indices[i] = Istart + i + A_size;
61: } else if (rank == 1) {
62: A_size = (Iend - Istart) / 2;
63: B_size = (Iend - Istart) / 2;
64: PetscCalloc1(A_size, &A_indices);
65: PetscCalloc1(B_size, &B_indices);
66: for (i = 0; i < A_size; ++i) A_indices[i] = Istart + i;
67: for (i = 0; i < B_size; ++i) B_indices[i] = Istart + i + A_size;
68: }
70: ISCreateGeneral(PETSC_COMM_WORLD, A_size, A_indices, PETSC_OWN_POINTER, &A_IS);
71: ISCreateGeneral(PETSC_COMM_WORLD, B_size, B_indices, PETSC_OWN_POINTER, &B_IS);
72: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]: A_size = %" PetscInt_FMT ", B_size = %" PetscInt_FMT "\n", rank, A_size, B_size);
73: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
75: /* Solve the system */
76: KSPCreate(PETSC_COMM_WORLD, &ksp);
77: KSPSetType(ksp, KSPGMRES);
78: KSPSetOperators(ksp, A, A);
80: /* Define the fieldsplit for the global matrix */
81: KSPGetPC(ksp, &pc);
82: PCSetType(pc, PCFIELDSPLIT);
83: PCFieldSplitSetIS(pc, "a", A_IS);
84: PCFieldSplitSetIS(pc, "b", B_IS);
85: ISSetBlockSize(A_IS, block_size);
86: ISSetBlockSize(B_IS, block_size);
88: KSPSetFromOptions(ksp);
89: KSPSetUp(ksp);
91: ISDestroy(&A_IS);
92: ISDestroy(&B_IS);
93: KSPDestroy(&ksp);
94: MatDestroy(&A);
95: PetscFinalize();
96: return 0;
97: }
99: /*TEST
101: test:
102: nsize: 2
103: args: -method 1
105: test:
106: suffix: 2
107: nsize: 2
108: args: -method 2
110: TEST*/