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 PetscScalar :: alpha 12 PetscInt, intent(in) :: n,ii(*),jj(0:*) 13 14 PetscInt i,j,jstart,jend 15 16 jend = ii(1) 17 do i=1,n 18 jstart = jend 19 jend = ii(i+1) 20 alpha = x(i-1) 21 do j=jstart,jend-1 22 y(jj(j)) = y(jj(j)) + alpha*a(j) 23 end do 24 end do 25end subroutine FortranMultTransposeAddAIJ 26 27pure subroutine FortranMultAIJ(n,x,ii,jj,a,y) 28 implicit none (type, external) 29 PetscScalar, intent(in) :: x(0:*),a(0:*) 30 PetscScalar, intent(inout) :: y(*) 31 PetscInt, intent(in) :: n,ii(*),jj(0:*) 32 33 PetscInt i,j,jstart,jend 34 PetscScalar sum 35 36#ifdef PETSC_USE_OPENMP_KERNELS 37 !omp parallel do private(j,jstart,jend,sum) 38#endif 39 do i=1,n 40 jstart = ii(i) 41 jend = ii(i+1) 42 sum = 0.d0 43 do j=jstart,jend-1 44 sum = sum + a(j)*x(jj(j)) 45 end do 46 y(i) = sum 47 end do 48end subroutine FortranMultAIJ 49