!
!    Fortran kernel for gemv() BLAS operation. This version supports
!  matrix array stored in single precision but vectors in double
!
#include <petsc/finclude/petscsys.h>

pure subroutine MSGemv(bs,ncols,A,x,y)
  use, intrinsic :: ISO_C_binding
  implicit none (type, external)
  PetscInt, intent(in) :: bs,ncols
  MatScalar, intent(in) :: A(bs,ncols)
  PetscScalar, intent(in) :: x(ncols)
  PetscScalar, intent(out) :: y(bs)

  PetscInt :: i

  y(1:bs) = 0.0
  do i=1,ncols
    y(1:bs) = y(1:bs) + A(1:bs,i)*x(i)
  end do
end subroutine MSGemv

pure subroutine MSGemvp(bs,ncols,A,x,y)
  use, intrinsic :: ISO_C_binding
  implicit none (type, external)
  PetscInt, intent(in) :: bs,ncols
  MatScalar, intent(in) :: A(bs,ncols)
  PetscScalar, intent(in) :: x(ncols)
  PetscScalar, intent(inout) :: y(bs)

  PetscInt :: i

  do i=1,ncols
    y(1:bs) = y(1:bs) + A(1:bs,i)*x(i)
  end do
end subroutine MSGemvp

pure subroutine MSGemvm(bs,ncols,A,x,y)
  use, intrinsic :: ISO_C_binding
  implicit none (type, external)
  PetscInt, intent(in) :: bs,ncols
  MatScalar, intent(in) :: A(bs,ncols)
  PetscScalar, intent(in) :: x(ncols)
  PetscScalar, intent(inout) :: y(bs)

  PetscInt :: i

  do i=1,ncols
    y(1:bs) = y(1:bs) - A(1:bs,i)*x(i)
  end do
end subroutine MSGemvm

pure subroutine MSGemvt(bs,ncols,A,x,y)
  use, intrinsic :: ISO_C_binding
  implicit none (type, external)
  PetscInt, intent(in) :: bs,ncols
  MatScalar, intent(in) :: A(bs,ncols)
  PetscScalar, intent(in) :: x(bs)
  PetscScalar, intent(inout) :: y(ncols)

  PetscInt :: i

  do  i=1,ncols
    y(i) = y(i) + sum(A(1:bs,i)*x(1:bs))
  end do
end subroutine MSGemvt

pure subroutine MSGemm(bs,A,B,C)
  use, intrinsic :: ISO_C_binding
  implicit none (type, external)
  PetscInt, intent(in) :: bs
  MatScalar, intent(in) :: B(bs,bs),C(bs,bs)
  MatScalar, intent(inout) :: A(bs,bs)

  PetscInt :: i,j

  do i=1,bs
    do j=1,bs
      A(i,j) = A(i,j) - sum(B(i,1:bs)*C(1:bs,j))
    end do
  end do
end subroutine MSGemm

pure subroutine MSGemmi(bs,A,C,B)
  use, intrinsic :: ISO_C_binding
  implicit none (type, external)
  PetscInt, intent(in) :: bs
  MatScalar, intent(in) :: B(bs,bs),C(bs,bs)
  MatScalar, intent(out) :: A(bs,bs)

  PetscInt :: i,j

  do i=1,bs
    do j=1,bs
      A(i,j) = sum(B(i,1:bs)*C(1:bs,j))
    end do
  end do
end subroutine MSGemmi