! Test MatCreateNest() with NULL index sets
#include <petsc/finclude/petscmat.h>
program main
  use petscmat
  implicit none

  Mat                   :: A, D, Id, Acopy
  Mat, dimension(4)      :: mats
  Vec                   :: v, w
  PetscInt              :: i, rstart, rend
  PetscInt, parameter    :: n = 6, nb = 2
  PetscScalar, parameter :: one = 1.0, two = 2.0
  PetscErrorCode        :: ierr

  PetscCallA(PetscInitialize(ierr))

  PetscCallA(MatCreateConstantDiagonal(PETSC_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, n, n, one, Id, ierr))
  PetscCallA(MatCreateVecs(Id, v, w, ierr))
  PetscCallA(VecGetOwnershipRange(v, rstart, rend, ierr))
  do i = rstart, rend - 1
    PetscCallA(VecSetValue(v, i, two/(i + 1), INSERT_VALUES, ierr))
  end do
  PetscCallA(VecAssemblyBegin(v, ierr))
  PetscCallA(VecAssemblyEnd(v, ierr))
  PetscCallA(MatCreateDiagonal(v, D, ierr))

  mats = [PETSC_NULL_MAT, D, Id, PETSC_NULL_MAT]
  PetscCallA(MatCreateNest(PETSC_COMM_WORLD, nb, PETSC_NULL_IS_ARRAY, nb, PETSC_NULL_IS_ARRAY, mats, A, ierr))
  PetscCallA(MatView(A, PETSC_VIEWER_STDOUT_WORLD, ierr))
  PetscCallA(MatNestSetSubMats(A, nb, PETSC_NULL_IS_ARRAY, nb, PETSC_NULL_IS_ARRAY, mats, ierr))
  PetscCallA(MatView(A, PETSC_VIEWER_STDOUT_WORLD, ierr))

  ! test MatCopy()
  PetscCallA(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, Acopy, ierr))
  PetscCallA(MatCopy(A, Acopy, DIFFERENT_NONZERO_PATTERN, ierr))

  PetscCallA(MatDestroy(Acopy, ierr))
  PetscCallA(MatDestroy(Id, ierr))
  PetscCallA(VecDestroy(v, ierr))
  PetscCallA(VecDestroy(w, ierr))
  PetscCallA(MatDestroy(D, ierr))
  PetscCallA(MatDestroy(A, ierr))
  PetscCallA(PetscFinalize(ierr))

end program

!/*TEST
!
!   test:
!      nsize: 2
!
!TEST*/