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*/