Actual source code: sgemv.F90

      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>

      7: pure 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)

     15:   PetscInt :: i

     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
     21: end subroutine MSGemv

     23: pure 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)

     31:   PetscInt :: i

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

     38: pure 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)

     46:   PetscInt :: i

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

     53: pure 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)

     61:   PetscInt :: i

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

     68: pure 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)

     75:   PetscInt :: i, j

     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
     82: end subroutine MSGemm

     84: pure 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)

     91:   PetscInt :: i, j

     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
     98: end subroutine MSGemmi