Actual source code: ex79f.F90

  1: !
  2: !   This program demonstrates use of MatGetRowIJ() from Fortran
  3: !
  4:       program main

  6: #include <petsc/finclude/petscmat.h>
  7:       use petscmat
  8:       implicit none

 10:       Mat         A,Ad,Ao
 11:       PetscErrorCode ierr
 12:       PetscMPIInt rank
 13:       PetscViewer v
 14:       PetscInt i,j,ia(1),ja(1)
 15:       PetscInt n,icol(1),rstart
 16:       PetscInt zero,one,rend
 17:       PetscBool   done
 18:       PetscOffset iia,jja,aaa,iicol
 19:       PetscScalar aa(1)

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

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

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

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

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

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

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