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:   num_cells = merge(1, 0, rank == 0)

 22:   PetscCallA(VecCreate(PETSC_COMM_WORLD, v, ierr))
 23:   PetscCallA(VecSetSizes(v, num_cells*blocksize, PETSC_DECIDE, ierr))
 24:   PetscCallA(VecSetBlockSize(v, blocksize, ierr))
 25:   PetscCallA(VecSetFromOptions(v, ierr))

 27:   subsize = num_cells
 28:   allocate (subindices(0:subsize - 1))
 29:   subindices = [(i, i=0, subsize - 1)]*blocksize + field
 30:   PetscCallA(ISCreateGeneral(PETSC_COMM_WORLD, subsize, subindices, PETSC_COPY_VALUES, index_set, ierr))
 31:   deallocate (subindices)

 33:   PetscCallA(VecGetSubVector(v, index_set, subv, ierr))
 34:   PetscCallA(VecRestoreSubVector(v, index_set, subv, ierr))
 35:   PetscCallA(ISDestroy(index_set, ierr))

 37:   PetscCallA(VecDestroy(v, ierr))
 38:   PetscCallA(PetscFinalize(ierr))
 39: end

 41: !/*TEST
 42: !
 43: !   test:
 44: !      nsize: 2
 45: !      filter: sort -b
 46: !      filter_output: sort -b
 47: !      output_file: output/empty.out
 48: !
 49: !TEST*/