!
!   This program demonstrates use of MatGetRowIJ() from Fortran
!
#include <petsc/finclude/petscmat.h>
program main
  use petscmat
  implicit none

  Mat A, Ad, Ao
  PetscErrorCode ierr
  PetscMPIInt rank
  PetscViewer v
  PetscInt i, j
  PetscInt n, rstart
  PetscInt zero, one, rend
  PetscBool done, bb
  PetscScalar, pointer :: aa(:)
  PetscInt, pointer :: ia(:), ja(:), icol(:)

  PetscCallA(PetscInitialize(ierr))
  PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))

  PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, '${PETSC_DIR}/share/petsc/datafiles/matrices/'//'ns-real-int32-float64', FILE_MODE_READ, v, ierr))
  PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
  PetscCallA(MatSetType(A, MATMPIAIJ, ierr))
  PetscCallA(MatLoad(A, v, ierr))
  PetscCallA(MatView(A, PETSC_VIEWER_STDOUT_WORLD, ierr))

  PetscCallA(MatMPIAIJGetSeqAIJ(A, Ad, Ao, icol, ierr))
  PetscCallA(MatGetOwnershipRange(A, rstart, rend, ierr))
!
!   Print diagonal portion of matrix
!
  PetscCallA(PetscSequentialPhaseBegin(PETSC_COMM_WORLD, 1, ierr))
  zero = 0
  one = 1
  bb = .true.
  PetscCallA(MatGetRowIJ(Ad, one, bb, bb, n, ia, ja, done, ierr))
  PetscCallA(MatSeqAIJGetArray(Ad, aa, ierr))
  do i = 1, n
    write (7 + rank, *) 'row ', i + rstart, ' number nonzeros ', ia(i + 1) - ia(i)
    do j = ia(i), ia(i + 1) - 1
      write (7 + rank, *) '  ', j, ja(j) + rstart, aa(j)
    end do
  end do
  PetscCallA(MatRestoreRowIJ(Ad, one, bb, bb, n, ia, ja, done, ierr))
  PetscCallA(MatSeqAIJRestoreArray(Ad, aa, ierr))
!
!   Print off-diagonal portion of matrix
!
  PetscCallA(MatGetRowIJ(Ao, one, bb, bb, n, ia, ja, done, ierr))
  PetscCallA(MatSeqAIJGetArray(Ao, aa, ierr))
  do i = 1, n
    write (7 + rank, *) 'row ', i + rstart, ' number nonzeros ', ia(i + 1) - ia(i)
    do j = ia(i), ia(i + 1) - 1
      write (7 + rank, *) '  ', j, icol(ja(j)) + 1, aa(j)
    end do
  end do
  PetscCallA(MatMPIAIJRestoreSeqAIJ(A, Ad, Ao, icol, ierr))
  PetscCallA(MatRestoreRowIJ(Ao, one, bb, bb, n, ia, ja, done, ierr))
  PetscCallA(MatSeqAIJRestoreArray(Ao, aa, ierr))

  PetscCallA(PetscSequentialPhaseEnd(PETSC_COMM_WORLD, 1, ierr))

  PetscCallA(MatGetDiagonalBlock(A, Ad, ierr))
  PetscCallA(MatView(Ad, PETSC_VIEWER_STDOUT_WORLD, ierr))

  PetscCallA(MatView(A, PETSC_VIEWER_STDOUT_WORLD, ierr))
  PetscCallA(MatDestroy(A, ierr))
  PetscCallA(PetscViewerDestroy(v, ierr))
  PetscCallA(PetscFinalize(ierr))
end

!/*TEST
!
!     build:
!       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
!
!     test:
!        args: -binary_read_double -options_left false
!
!TEST*/