1! 2! 3! Fortran kernel for sparse matrix-vector product in the AIJ matrix format 4! 5#include <petsc/finclude/petscsys.h> 6! 7pure subroutine FortranMultTransposeAddAIJ(n,x,ii,jj,a,y) 8 implicit none (type, external) 9 PetscScalar, intent(in) :: x(0:*),a(0:*) 10 PetscScalar, intent(inout) :: y(0:*) 11 PetscInt, intent(in) :: n,ii(*),jj(0:*) 12 13 PetscInt :: i,jstart,jend 14 15 jend = ii(1) 16 do i=1,n 17 jstart = jend 18 jend = ii(i+1) 19 y(jj(jstart:jend-1)) = y(jj(jstart:jend-1)) + x(i-1)*a(jstart:jend-1) 20 end do 21end subroutine FortranMultTransposeAddAIJ 22 23pure subroutine FortranMultAIJ(n,x,ii,jj,a,y) 24 implicit none (type, external) 25 PetscScalar, intent(in) :: x(0:*),a(0:*) 26 PetscScalar, intent(inout) :: y(*) 27 PetscInt, intent(in) :: n,ii(*),jj(0:*) 28 29 PetscInt :: i,jstart,jend 30 31#ifdef PETSC_USE_OPENMP_KERNELS 32 !omp parallel do private(jstart,jend) 33#endif 34 do i=1,n 35 jstart = ii(i) 36 jend = ii(i+1) 37 y(i) = sum(a(jstart:jend-1)*x(jj(jstart:jend-1))) 38 end do 39end subroutine FortranMultAIJ 40