Actual source code: ex53.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

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