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
 15:   PetscInt zero, one, rend
 16:   PetscBool done, bb
 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:   zero = 0
 36:   one = 1
 37:   bb = .true.
 38:   PetscCallA(MatGetRowIJ(Ad, one, bb, bb, n, ia, ja, done, ierr))
 39:   PetscCallA(MatSeqAIJGetArray(Ad, aa, ierr))
 40:   do i = 1, n
 41:     write (7 + rank, *) 'row ', i + rstart, ' number nonzeros ', ia(i + 1) - ia(i)
 42:     do j = ia(i), ia(i + 1) - 1
 43:       write (7 + rank, *) '  ', j, ja(j) + rstart, aa(j)
 44:     end do
 45:   end do
 46:   PetscCallA(MatRestoreRowIJ(Ad, one, bb, bb, n, ia, ja, done, ierr))
 47:   PetscCallA(MatSeqAIJRestoreArray(Ad, aa, ierr))
 48: !
 49: !   Print off-diagonal portion of matrix
 50: !
 51:   PetscCallA(MatGetRowIJ(Ao, one, bb, bb, n, ia, ja, done, ierr))
 52:   PetscCallA(MatSeqAIJGetArray(Ao, aa, ierr))
 53:   do i = 1, n
 54:     write (7 + rank, *) 'row ', i + rstart, ' number nonzeros ', ia(i + 1) - ia(i)
 55:     do j = ia(i), ia(i + 1) - 1
 56:       write (7 + rank, *) '  ', j, icol(ja(j)) + 1, aa(j)
 57:     end do
 58:   end do
 59:   PetscCallA(MatMPIAIJRestoreSeqAIJ(A, Ad, Ao, icol, ierr))
 60:   PetscCallA(MatRestoreRowIJ(Ao, one, bb, bb, n, ia, ja, done, ierr))
 61:   PetscCallA(MatSeqAIJRestoreArray(Ao, aa, ierr))

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

 65:   PetscCallA(MatGetDiagonalBlock(A, Ad, ierr))
 66:   PetscCallA(MatView(Ad, PETSC_VIEWER_STDOUT_WORLD, ierr))

 68:   PetscCallA(MatView(A, PETSC_VIEWER_STDOUT_WORLD, ierr))
 69:   PetscCallA(MatDestroy(A, ierr))
 70:   PetscCallA(PetscViewerDestroy(v, ierr))
 71:   PetscCallA(PetscFinalize(ierr))
 72: end

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