1: !
  2: !  Test VecGetSubVector()
  3: !  Contributed-by: Adrian Croucher <gitlab@mg.gitlab.com>
  4: #include <petsc/finclude/petsc.h>
  5: program main
  6:   use petsc
  7:   implicit none

  9:   PetscMPIInt :: rank
 10:   PetscErrorCode :: ierr
 11:   PetscInt :: num_cells, subsize, i
 12:   PetscInt, parameter :: blocksize = 3, field = 0
 13:   Vec :: v, subv
 14:   IS :: index_set
 15:   PetscInt, allocatable :: subindices(:)

 17:   PetscCallA(PetscInitialize(ierr))
 18:   PetscCallMPIA(MPI_COMM_RANK(PETSC_COMM_WORLD, rank, ierr))

 20:   if (rank == 0) then
 21:     num_cells = 1
 22:   else
 23:     num_cells = 0
 24:   end if

 26:   PetscCallA(VecCreate(PETSC_COMM_WORLD, v, ierr))
 27:   PetscCallA(VecSetSizes(v, num_cells*blocksize, PETSC_DECIDE, ierr))
 28:   PetscCallA(VecSetBlockSize(v, blocksize, ierr))
 29:   PetscCallA(VecSetFromOptions(v, ierr))

 31:   subsize = num_cells
 32:   allocate (subindices(0:subsize - 1))
 33:   subindices = [(i, i=0, subsize - 1)]*blocksize + field
 34:   PetscCallA(ISCreateGeneral(PETSC_COMM_WORLD, subsize, subindices, PETSC_COPY_VALUES, index_set, ierr))
 35:   deallocate (subindices)

 37:   PetscCallA(VecGetSubVector(v, index_set, subv, ierr))
 38:   PetscCallA(VecRestoreSubVector(v, index_set, subv, ierr))
 39:   PetscCallA(ISDestroy(index_set, ierr))

 41:   PetscCallA(VecDestroy(v, ierr))
 42:   PetscCallA(PetscFinalize(ierr))
 43: end

 45: !/*TEST
 46: !
 47: !   test:
 48: !      nsize: 2
 49: !      filter: sort -b
 50: !      filter_output: sort -b
 51: !      output_file: output/empty.out
 52: !
 53: !TEST*/