SUBROUTINE petsc_create(a,b,bs,petsc_num,dnnz,onnz,nloc,nglob,nhalo)
      IMPLICIT NONE
      INTEGER*4, INTENT(IN), DIMENSION(:) :: petsc_num,dnnz,onnz
      INTEGER*4, INTENT(IN)  :: bs,nloc,nglob,nhalo
      INTEGER*4              :: ierr
      Vec, INTENT(INOUT)     :: b
      Mat, INTENT(INOUT)     :: a
      ISLocalToGlobalMapping :: rowmap, colmap

      !Matrix creation
      CALL MatCreateBAIJ(PETSC_COMM_WORLD,bs,nloc*bs,nloc*bs,nglob*bs,nglob*bs,&
                         PETSC_NULL_INTEGER,dnnz,PETSC_NULL_INTEGER,onnz,a,ierr)

      CALL MatSetOption(a,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr)
      CALL MatSetOption(a,MAT_NO_OFF_PROC_ZERO_ROWS,PETSC_TRUE,ierr)

      !Mapping local numbering to global PETSc numbering
      !Row mapping
      CALL ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,bs,nloc,petsc_num(1:nloc),PETSC_COPY_VALUES,rowmap,ierr)
      !Column mapping
      CALL ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,bs,nloc+nhalo,petsc_num,PETSC_COPY_VALUES,colmap,ierr)

      !Matrix finalization
      CALL MatSetLocalToGlobalMapping(a,rowmap,colmap,ierr)

      !Vector creation (RHS)
      CALL VecCreate(PETSC_COMM_WORLD,b,ierr)
      CALL VecSetSizes(b,nloc*bs,nglob*bs,ierr)
      CALL VecSetBlockSize(b,bs,ierr)
      CALL VecSetType(b,"mpi",ierr)
      CALL VecSetLocalToGlobalMapping(b,rowmap,ierr)

      CALL ISLocalToGlobalMappingDestroy(rowmap,ierr)
      CALL ISLocalToGlobalMappingDestroy(colmap,ierr)

   ENDSUBROUTINE petsc_create