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 implicit none (type, external) 9 PetscInt, intent(in) :: bs,ncols 10 MatScalar, intent(in) :: A(bs,ncols) 11 PetscScalar, intent(in) :: x(ncols) 12 PetscScalar, intent(out) :: y(bs) 13 14 PetscInt :: i 15 16 y(1:bs) = 0.0 17 do i=1,ncols 18 y(1:bs) = y(1:bs) + A(1:bs,i)*x(i) 19 end do 20end subroutine MSGemv 21 22pure subroutine MSGemvp(bs,ncols,A,x,y) 23 implicit none (type, external) 24 PetscInt, intent(in) :: bs,ncols 25 MatScalar, intent(in) :: A(bs,ncols) 26 PetscScalar, intent(in) :: x(ncols) 27 PetscScalar, intent(inout) :: y(bs) 28 29 PetscInt :: i 30 31 do i=1,ncols 32 y(1:bs) = y(1:bs) + A(1:bs,i)*x(i) 33 end do 34end subroutine MSGemvp 35 36pure subroutine MSGemvm(bs,ncols,A,x,y) 37 implicit none (type, external) 38 PetscInt, intent(in) :: bs,ncols 39 MatScalar, intent(in) :: A(bs,ncols) 40 PetscScalar, intent(in) :: x(ncols) 41 PetscScalar, intent(inout) :: y(bs) 42 43 PetscInt :: i 44 45 do i=1,ncols 46 y(1:bs) = y(1:bs) - A(1:bs,i)*x(i) 47 end do 48end subroutine MSGemvm 49 50pure subroutine MSGemvt(bs,ncols,A,x,y) 51 implicit none (type, external) 52 PetscInt, intent(in) :: bs,ncols 53 MatScalar, intent(in) :: A(bs,ncols) 54 PetscScalar, intent(in) :: x(bs) 55 PetscScalar, intent(inout) :: y(ncols) 56 57 PetscInt :: i 58 59 do i=1,ncols 60 y(i) = y(i) + sum(A(1:bs,i)*x(1:bs)) 61 end do 62end subroutine MSGemvt 63 64pure subroutine MSGemm(bs,A,B,C) 65 implicit none (type, external) 66 PetscInt, intent(in) :: bs 67 MatScalar, intent(in) :: B(bs,bs),C(bs,bs) 68 MatScalar, intent(inout) :: A(bs,bs) 69 70 PetscInt :: i,j 71 72 do i=1,bs 73 do j=1,bs 74 A(i,j) = A(i,j) - sum(B(i,1:bs)*C(1:bs,j)) 75 end do 76 end do 77end subroutine MSGemm 78 79pure subroutine MSGemmi(bs,A,C,B) 80 implicit none (type, external) 81 PetscInt, intent(in) :: bs 82 MatScalar, intent(in) :: B(bs,bs),C(bs,bs) 83 MatScalar, intent(out) :: A(bs,bs) 84 85 PetscInt :: i,j 86 87 do i=1,bs 88 do j=1,bs 89 A(i,j) = sum(B(i,1:bs)*C(1:bs,j)) 90 end do 91 end do 92end subroutine MSGemmi 93