#include "petsc/finclude/petscdmda.h"
#include "petsc/finclude/petscdm.h"
#include "petsc/finclude/petscmat.h"
#include "petsc/finclude/petscdmcomposite.h"
program matnest_bug_reproducer
  use petscsys
  use petscmat
  use petscdm
  use petscdmda
  use petscdmcomposite

  implicit none (type,external)

  type(tDM)          :: dm, dmc
  type(tMat)         :: A, A_null, N, matarray(4)
  type(tIS), pointer :: isc(:)
  PetscErrorCode     :: ierr

  ! Initialize
  PetscCallA(PetscInitialize(PETSC_NULL_CHARACTER,ierr))

  ! Create DM
  PetscCallA(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE_INTEGER, PETSC_DECIDE_INTEGER, 1,0, PETSC_NULL_INTEGER_ARRAY, PETSC_NULL_INTEGER_ARRAY, dm,ierr))
  PetscCallA(DMSetup(dm,ierr))

  ! Create non-null matrix
  PetscCallA(DMCreateMatrix(dm,A,ierr))

  ! Create composite DM
  PetscCallA(DMCompositeCreate(PETSC_COMM_WORLD,dmc,ierr))
  PetscCallA(DMCompositeAddDM(dmc,dm,ierr))
  PetscCallA(DMCompositeAddDM(dmc,dm,ierr))
  PetscCallA(DMCompositeAddDM(dmc,dm,ierr))
  PetscCallA(DMCompositeAddDM(dmc,dm,ierr))
  PetscCallA(DMSetup(dmc,ierr))

  ! This works: matarray does not have null entries
  ! matarray = [A,A,A,A]
  ! PetscCallA(DMCompositeGetGlobalISs(dmc,isc,ierr))
  ! PetscCallA(MatCreateNest(PETSC_COMM_WORLD,2,isc,2,isc,matArray,N,ierr))
  ! PetscCallA(DMCompositeRestoreGlobalISs(dmc,isc,ierr))
  ! PetscCallA(MatView(N,PETSC_VIEWER_STDOUT_WORLD,ierr))
  ! PetscCallA(MatDestroy(N,ierr))

  ! This does not work: matarray has null entries
  matarray = [A,A,A,A_null]
  PetscCallA(DMCompositeGetGlobalISs(dmc,isc,ierr))
  PetscCallA(MatCreateNest(PETSC_COMM_WORLD,2,isc,2,isc,matArray,N,ierr))
  PetscCallA(DMCompositeRestoreGlobalISs(dmc,isc,ierr))
  PetscCallA(MatView(N,PETSC_VIEWER_STDOUT_WORLD,ierr))
  PetscCallA(MatDestroy(N,ierr))

  ! Destroy
  PetscCallA(MatDestroy(A,ierr))
  PetscCallA(MatDestroy(A_null,ierr))
  PetscCallA(DMDestroy(dmc,ierr))
  PetscCallA(DMDestroy(dm,ierr))

  ! Finalize
  PetscCallA(PetscFinalize(ierr))

end program