Actual source code: ex79f.F90

petsc-3.12.0 2019-09-29
Report Typos and Errors
  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:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 22:       if (ierr .ne. 0) then
 23:         print*,'Unable to initialize PETSc'
 24:         stop
 25:       endif
 26:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

 28:       call PetscViewerBinaryOpen(PETSC_COMM_WORLD,                          &
 29:      & '${PETSC_DIR}/share/petsc/datafiles/matrices/' //                       &
 30:      & 'ns-real-int32-float64',                                               &
 31:      &                          FILE_MODE_READ,v,ierr)
 32:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
 33:       call MatSetType(A, MATMPIAIJ,ierr)
 34:       call MatLoad(A,v,ierr)
 35:       CHKERRA(ierr)
 36:       call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr)

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

 72:       call PetscSequentialPhaseEnd(PETSC_COMM_WORLD,1,ierr)

 74:       call MatGetDiagonalBlock(A,Ad,ierr)
 75:       call MatView(Ad,PETSC_VIEWER_STDOUT_WORLD,ierr)

 77:       call MatView(A,PETSC_VIEWER_STDOUT_WORLD,ierr)
 78:       call MatDestroy(A,ierr)
 79:       call PetscViewerDestroy(v,ierr)
 80:       call PetscFinalize(ierr)
 81:       end

 83: !/*TEST
 84: !
 85: !     build:
 86: !       requires: double !complex !define(PETSC_USE_64BIT_INDICES)
 87: !
 88: !     test:
 89: !        args: -binary_read_double -options_left false
 90: !
 91: !TEST*/