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;
18: PetscInitialize(&argc,&argv,(char*)0,help);
19: MPI_Comm_rank(MPI_COMM_WORLD,&rank);
21: PetscOptionsGetInt(PETSC_NULL,PETSC_NULL,"-mat_size",&mat_size,PETSC_NULL);
22: PetscOptionsGetInt(PETSC_NULL,PETSC_NULL,"-method",&method,PETSC_NULL);
23: PetscOptionsGetInt(PETSC_NULL,PETSC_NULL,"-block_size",&block_size,PETSC_NULL);
25: if (rank == 0) {
26: PetscPrintf(PETSC_COMM_SELF," matrix size = %D, block size = %D\n",mat_size,block_size);
27: }
29: MatCreate(PETSC_COMM_WORLD,&A);
30: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,mat_size,mat_size);
31: MatSetType(A,MATMPIAIJ);
32: MatSetUp(A);
34: MatGetOwnershipRange(A,&Istart,&Iend);
36: for (i = Istart; i < Iend; ++i) {
37: MatSetValue(A,i,i,2,INSERT_VALUES);
38: if (i < mat_size-1) {
39: MatSetValue(A,i,i+1,-1,INSERT_VALUES);
40: }
41: if (i > 0) {
42: MatSetValue(A,i,i-1,-1,INSERT_VALUES);
43: }
44: }
46: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
47: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
49: MatGetLocalSize(A,&local_m,&local_n);
51: /* Create Index Sets */
52: if (rank == 0) {
53: if (method > 1) {
54: /* with method > 1, the fieldsplit B is set to zero */
55: A_size = Iend-Istart;
56: B_size = 0;
57: } else {
58: /* 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 */
59: A_size = (Iend-Istart)/2;
60: B_size = (Iend-Istart)/2;
61: }
62: PetscCalloc1(A_size,&A_indices);
63: PetscCalloc1(B_size,&B_indices);
64: for (i = 0; i < A_size; ++i) A_indices[i] = Istart + i;
65: for (i = 0; i < B_size; ++i) B_indices[i] = Istart + i + A_size;
66: } else if (rank == 1) {
67: A_size = (Iend-Istart)/2;
68: B_size = (Iend-Istart)/2;
69: PetscCalloc1(A_size,&A_indices);
70: PetscCalloc1(B_size,&B_indices);
71: for (i = 0; i < A_size; ++i) A_indices[i] = Istart + i;
72: for (i = 0; i < B_size; ++i) B_indices[i] = Istart + i + A_size;
73: }
75: ISCreateGeneral(PETSC_COMM_WORLD,A_size,A_indices,PETSC_OWN_POINTER,&A_IS);
76: ISCreateGeneral(PETSC_COMM_WORLD,B_size,B_indices,PETSC_OWN_POINTER,&B_IS);
77: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d]: A_size = %D, B_size = %D\n",rank,A_size,B_size);
78: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
80: /* Solve the system */
81: KSPCreate(PETSC_COMM_WORLD,&ksp);
82: KSPSetType(ksp,KSPGMRES);
83: KSPSetOperators(ksp,A,A);
85: /* Define the fieldsplit for the global matrix */
86: KSPGetPC(ksp,&pc);
87: PCSetType(pc,PCFIELDSPLIT);
88: PCFieldSplitSetIS(pc,"a",A_IS);
89: PCFieldSplitSetIS(pc,"b",B_IS);
90: ISSetBlockSize(A_IS,block_size);
91: ISSetBlockSize(B_IS,block_size);
93: KSPSetFromOptions(ksp);
94: KSPSetUp(ksp);
96: ISDestroy(&A_IS);
97: ISDestroy(&B_IS);
98: KSPDestroy(&ksp);
99: MatDestroy(&A);
100: PetscFinalize();
101: return 0;
102: }
104: /*TEST
106: test:
107: nsize: 2
108: args: -method 1
110: test:
111: suffix: 2
112: nsize: 2
113: args: -method 2
115: TEST*/