xref: /petsc/src/mat/impls/aij/seq/ftn-kernels/fmult.F90 (revision ccfb0f9f40a0131988d7995ed9679700dae2a75a)
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  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:*)
13
14  PetscInt :: i, jstart, jend
15
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
22end subroutine FortranMultTransposeAddAIJ
23
24pure 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:*)
30
31  PetscInt :: i, jstart, jend
32
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
41end subroutine FortranMultAIJ
42