1! 2! This program demonstrates use of MatGetRowIJ() from Fortran 3! 4 program main 5 6#include <petsc/finclude/petscmat.h> 7 use petscmat 8 implicit none 9 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) 20 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) 27 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) 37 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) 71 72 call PetscSequentialPhaseEnd(PETSC_COMM_WORLD,1,ierr) 73 74 call MatGetDiagonalBlock(A,Ad,ierr) 75 call MatView(Ad,PETSC_VIEWER_STDOUT_WORLD,ierr) 76 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 82 83!/*TEST 84! 85! build: 86! requires: double !complex !defined(PETSC_USE_64BIT_INDICES) 87! 88! test: 89! args: -binary_read_double -options_left false 90! 91!TEST*/ 92