Actual source code: fmult.F90

      1: !
      2: !
      3: !    Fortran kernel for sparse matrix-vector product in the AIJ matrix format
      4: !
      5: #include <petsc/finclude/petscsys.h>
      6: !
      7: pure subroutine FortranMultTransposeAddAIJ(n, x, ii, jj, a, y)
      8:   use, intrinsic :: ISO_C_binding
      9:   implicit none(type, external)
     10:   PetscScalar, intent(in) :: x(0:*), a(0:*)
     11:   PetscScalar, intent(inout) :: y(0:*)
     12:   PetscInt, intent(in) :: n, ii(*), jj(0:*)

     14:   PetscInt :: i, jstart, jend

     16:   jend = ii(1)
     17:   do i = 1, n
     18:     jstart = jend
     19:     jend = ii(i + 1)
     20:     y(jj(jstart:jend - 1)) = y(jj(jstart:jend - 1)) + x(i - 1)*a(jstart:jend - 1)
     21:   end do
     22: end subroutine FortranMultTransposeAddAIJ

     24: pure subroutine FortranMultAIJ(n, x, ii, jj, a, y)
     25:   use, intrinsic :: ISO_C_binding
     26:   implicit none(type, external)
     27:   PetscScalar, intent(in) :: x(0:*), a(0:*)
     28:   PetscScalar, intent(inout) :: y(*)
     29:   PetscInt, intent(in) :: n, ii(*), jj(0:*)

     31:   PetscInt :: i, jstart, jend

     33: #ifdef PETSC_USE_OPENMP_KERNELS
     34:   !omp parallel do private(jstart,jend)
     35: #endif
     36:   do i = 1, n
     37:     jstart = ii(i)
     38:     jend = ii(i + 1)
     39:     y(i) = sum(a(jstart:jend - 1)*x(jj(jstart:jend - 1)))
     40:   end do
     41: end subroutine FortranMultAIJ