1: !
  2: !   This program demonstrates use of MatGetRowIJ() from Fortran
  3: !
  4: #include <petsc/finclude/petscmat.h>
  5: program main
  6:   use petscmat
  7:   implicit none

  9:   Mat A, Ad, Ao
 10:   PetscErrorCode ierr
 11:   PetscMPIInt rank
 12:   PetscViewer v
 13:   PetscInt i, j
 14:   PetscInt n, rstart, rend
 15:   PetscInt, parameter :: one = 1
 16:   PetscBool done
 17:   PetscScalar, pointer :: aa(:)
 18:   PetscInt, pointer :: ia(:), ja(:), icol(:)

 20:   PetscCallA(PetscInitialize(ierr))
 21:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))

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

 29:   PetscCallA(MatMPIAIJGetSeqAIJ(A, Ad, Ao, icol, ierr))
 30:   PetscCallA(MatGetOwnershipRange(A, rstart, rend, ierr))
 31: !
 32: !   Print diagonal portion of matrix
 33: !
 34:   PetscCallA(PetscSequentialPhaseBegin(PETSC_COMM_WORLD, 1, ierr))
 35:   PetscCallA(MatGetRowIJ(Ad, one, PETSC_TRUE, PETSC_TRUE, n, ia, ja, done, ierr))
 36:   PetscCallA(MatSeqAIJGetArray(Ad, aa, ierr))
 37:   do i = 1, n
 38:     write (7 + rank, *) 'row ', i + rstart, ' number nonzeros ', ia(i + 1) - ia(i)
 39:     do j = ia(i), ia(i + 1) - 1
 40:       write (7 + rank, *) '  ', j, ja(j) + rstart, aa(j)
 41:     end do
 42:   end do
 43:   PetscCallA(MatRestoreRowIJ(Ad, one, PETSC_TRUE, PETSC_TRUE, n, ia, ja, done, ierr))
 44:   PetscCallA(MatSeqAIJRestoreArray(Ad, aa, ierr))
 45: !
 46: !   Print off-diagonal portion of matrix
 47: !
 48:   PetscCallA(MatGetRowIJ(Ao, one, PETSC_TRUE, PETSC_TRUE, n, ia, ja, done, ierr))
 49:   PetscCallA(MatSeqAIJGetArray(Ao, aa, ierr))
 50:   do i = 1, n
 51:     write (7 + rank, *) 'row ', i + rstart, ' number nonzeros ', ia(i + 1) - ia(i)
 52:     do j = ia(i), ia(i + 1) - 1
 53:       write (7 + rank, *) '  ', j, icol(ja(j)) + 1, aa(j)
 54:     end do
 55:   end do
 56:   PetscCallA(MatMPIAIJRestoreSeqAIJ(A, Ad, Ao, icol, ierr))
 57:   PetscCallA(MatRestoreRowIJ(Ao, one, PETSC_TRUE, PETSC_TRUE, n, ia, ja, done, ierr))
 58:   PetscCallA(MatSeqAIJRestoreArray(Ao, aa, ierr))

 60:   PetscCallA(PetscSequentialPhaseEnd(PETSC_COMM_WORLD, 1, ierr))

 62:   PetscCallA(MatGetDiagonalBlock(A, Ad, ierr))
 63:   PetscCallA(MatView(Ad, PETSC_VIEWER_STDOUT_WORLD, ierr))

 65:   PetscCallA(MatView(A, PETSC_VIEWER_STDOUT_WORLD, ierr))
 66:   PetscCallA(MatDestroy(A, ierr))
 67:   PetscCallA(PetscViewerDestroy(v, ierr))
 68:   PetscCallA(PetscFinalize(ierr))
 69: end

 71: !/*TEST
 72: !
 73: !     build:
 74: !       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
 75: !
 76: !     test:
 77: !        args: -binary_read_double -options_left false
 78: !
 79: !TEST*/