1! 2! Fortran kernel for gemv() BLAS operation. This version supports 3! matrix array stored in single precision but vectors in double 4! 5#include <petsc/finclude/petscsys.h> 6 7pure subroutine MSGemv(bs, ncols, A, x, y) 8 use, intrinsic :: ISO_C_binding 9 implicit none(type, external) 10 PetscInt, intent(in) :: bs, ncols 11 MatScalar, intent(in) :: A(bs, ncols) 12 PetscScalar, intent(in) :: x(ncols) 13 PetscScalar, intent(out) :: y(bs) 14 15 PetscInt :: i 16 17 y(1:bs) = 0.0 18 do i = 1, ncols 19 y(1:bs) = y(1:bs) + A(1:bs, i)*x(i) 20 end do 21end subroutine MSGemv 22 23pure subroutine MSGemvp(bs, ncols, A, x, y) 24 use, intrinsic :: ISO_C_binding 25 implicit none(type, external) 26 PetscInt, intent(in) :: bs, ncols 27 MatScalar, intent(in) :: A(bs, ncols) 28 PetscScalar, intent(in) :: x(ncols) 29 PetscScalar, intent(inout) :: y(bs) 30 31 PetscInt :: i 32 33 do i = 1, ncols 34 y(1:bs) = y(1:bs) + A(1:bs, i)*x(i) 35 end do 36end subroutine MSGemvp 37 38pure subroutine MSGemvm(bs, ncols, A, x, y) 39 use, intrinsic :: ISO_C_binding 40 implicit none(type, external) 41 PetscInt, intent(in) :: bs, ncols 42 MatScalar, intent(in) :: A(bs, ncols) 43 PetscScalar, intent(in) :: x(ncols) 44 PetscScalar, intent(inout) :: y(bs) 45 46 PetscInt :: i 47 48 do i = 1, ncols 49 y(1:bs) = y(1:bs) - A(1:bs, i)*x(i) 50 end do 51end subroutine MSGemvm 52 53pure subroutine MSGemvt(bs, ncols, A, x, y) 54 use, intrinsic :: ISO_C_binding 55 implicit none(type, external) 56 PetscInt, intent(in) :: bs, ncols 57 MatScalar, intent(in) :: A(bs, ncols) 58 PetscScalar, intent(in) :: x(bs) 59 PetscScalar, intent(inout) :: y(ncols) 60 61 PetscInt :: i 62 63 do i = 1, ncols 64 y(i) = y(i) + sum(A(1:bs, i)*x(1:bs)) 65 end do 66end subroutine MSGemvt 67 68pure subroutine MSGemm(bs, A, B, C) 69 use, intrinsic :: ISO_C_binding 70 implicit none(type, external) 71 PetscInt, intent(in) :: bs 72 MatScalar, intent(in) :: B(bs, bs), C(bs, bs) 73 MatScalar, intent(inout) :: A(bs, bs) 74 75 PetscInt :: i, j 76 77 do i = 1, bs 78 do j = 1, bs 79 A(i, j) = A(i, j) - sum(B(i, 1:bs)*C(1:bs, j)) 80 end do 81 end do 82end subroutine MSGemm 83 84pure subroutine MSGemmi(bs, A, C, B) 85 use, intrinsic :: ISO_C_binding 86 implicit none(type, external) 87 PetscInt, intent(in) :: bs 88 MatScalar, intent(in) :: B(bs, bs), C(bs, bs) 89 MatScalar, intent(out) :: A(bs, bs) 90 91 PetscInt :: i, j 92 93 do i = 1, bs 94 do j = 1, bs 95 A(i, j) = sum(B(i, 1:bs)*C(1:bs, j)) 96 end do 97 end do 98end subroutine MSGemmi 99