Actual source code: ex53.c
petsc-3.12.0 2019-09-29
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;
12: KSP ksp;
13: PC pc;
14: PetscInt Istart,Iend,local_m,local_n,i;
15: PetscMPIInt rank;
16: PetscInt method=2,mat_size=40,block_size=2,*A_indices=NULL,*B_indices=NULL,A_size=0,B_size=0;
17: IS A_IS, B_IS;
19: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
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) {
27: PetscPrintf(PETSC_COMM_SELF," matrix size = %D, block size = %D\n",mat_size,block_size);
28: }
30: MatCreate(PETSC_COMM_WORLD,&A);
31: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,mat_size,mat_size);
32: MatSetType(A,MATMPIAIJ);
33: MatSetUp(A);
35: MatGetOwnershipRange(A,&Istart,&Iend);
37: for (i = Istart; i < Iend; ++i) {
38: MatSetValue(A,i,i,2,INSERT_VALUES);
39: if (i < mat_size-1) {
40: MatSetValue(A,i,i+1,-1,INSERT_VALUES);
41: }
42: if (i > 0) {
43: MatSetValue(A,i,i-1,-1,INSERT_VALUES);
44: }
45: }
47: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
48: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
50: MatGetLocalSize(A,&local_m,&local_n);
52: /* Create Index Sets */
53: if (rank == 0) {
54: if (method > 1) {
55: /* with method > 1, the fieldsplit B is set to zero */
56: A_size = Iend-Istart;
57: B_size = 0;
58: } else {
59: /* 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 */
60: A_size = (Iend-Istart)/2;
61: B_size = (Iend-Istart)/2;
62: }
63: PetscCalloc1(A_size,&A_indices);
64: PetscCalloc1(B_size,&B_indices);
65: for (i = 0; i < A_size; ++i) A_indices[i] = Istart + i;
66: for (i = 0; i < B_size; ++i) B_indices[i] = Istart + i + A_size;
67: } else if (rank == 1) {
68: A_size = (Iend-Istart)/2;
69: B_size = (Iend-Istart)/2;
70: PetscCalloc1(A_size,&A_indices);
71: PetscCalloc1(B_size,&B_indices);
72: for (i = 0; i < A_size; ++i) A_indices[i] = Istart + i;
73: for (i = 0; i < B_size; ++i) B_indices[i] = Istart + i + A_size;
74: }
76: ISCreateGeneral(PETSC_COMM_WORLD,A_size,A_indices,PETSC_OWN_POINTER,&A_IS);
77: ISCreateGeneral(PETSC_COMM_WORLD,B_size,B_indices,PETSC_OWN_POINTER,&B_IS);
78: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d]: A_size = %D, B_size = %D\n",rank,A_size,B_size);
79: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
81: /* Solve the system */
82: KSPCreate(PETSC_COMM_WORLD,&ksp);
83: KSPSetType(ksp,KSPGMRES);
84: KSPSetOperators(ksp,A,A);
86: /* Define the fieldsplit for the global matrix */
87: KSPGetPC(ksp,&pc);
88: PCSetType(pc,PCFIELDSPLIT);
89: PCFieldSplitSetIS(pc,"a",A_IS);
90: PCFieldSplitSetIS(pc,"b",B_IS);
91: ISSetBlockSize(A_IS,block_size);
92: ISSetBlockSize(B_IS,block_size);
94: KSPSetFromOptions(ksp);
95: KSPSetUp(ksp);
97: ISDestroy(&A_IS);
98: ISDestroy(&B_IS);
99: KSPDestroy(&ksp);
100: MatDestroy(&A);
101: PetscFinalize();
102: return ierr;
103: }
106: /*TEST
108: test:
109: nsize: 2
110: args: -method 1
112: test:
113: suffix: 2
114: nsize: 2
115: args: -method 2
117: TEST*/