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 PetscCallA(PetscInitialize(ierr)) 22 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)) 23 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)) 29 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)) 61 62 PetscCallA(PetscSequentialPhaseEnd(PETSC_COMM_WORLD,1,ierr)) 63 64 PetscCallA(MatGetDiagonalBlock(A,Ad,ierr)) 65 PetscCallA(MatView(Ad,PETSC_VIEWER_STDOUT_WORLD,ierr)) 66 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 72 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*/ 82