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