prod_elemental_x Subroutine

   subroutine prod_elemental_x(n, nz, nelt, nvar, x, y, a_elt, eltptr, eltvar)
   !! \([A] \{x\}\), elemental CC format
   implicit none
   integer(kind=I4), intent(in) :: n, nz, nelt, nvar
   real(kind=R8),    dimension(nz     ), intent(in) :: a_elt
   integer(kind=I4), dimension(nelt +1), intent(in) :: eltptr
   integer(kind=I4), dimension(nvar   ), intent(in) :: eltvar
   real(kind=R8), dimension(n), intent(in)      :: x
   real(kind=R8), dimension(n), intent(out)     :: y
      integer(kind=I4) :: i, j, k, kk, inc_nz, inc_nn, n_elem, i_elem, max_n_elem
      real(kind=R8), dimension(:), allocatable :: a_elt_tmp, x_tmp, y_tmp
      inc_nz = 0
      inc_nn = 0
      n_elem = 0
      y      = 0._R8

      max_n_elem = 0
      do i_elem = 1, nelt
         max_n_elem = max(max_n_elem, eltptr(i_elem +1) -eltptr(i_elem))
      enddo

      allocate( a_elt_tmp(1:max_n_elem**2) )
      allocate(     x_tmp(1:max_n_elem   ) )
      allocate(     y_tmp(1:max_n_elem   ) )

      do i_elem = 1, nelt                                               ! browse all elemental matrices
         inc_nn = inc_nn +n_elem                                        ! step in eltvar for matrix number i_elem
         inc_nz = inc_nz +n_elem**2                                     ! step in a_elt for matrix number i_elem
         n_elem = eltptr(i_elem +1) -eltptr(i_elem)                     ! elemental matrix size
         a_elt_tmp(1:n_elem**2) = a_elt(inc_nz +1 : inc_nz +n_elem**2)  ! elemental matrix coefficients
         ! --- elemental rhs
         kk = 0
         do k = inc_nn +1, inc_nn +n_elem
            kk = kk +1
            x_tmp(kk) = x(eltvar(k))
         enddo
         ! --- elemental product
         y_tmp(1:n_elem) = 0._R8
         do i = 1, n_elem
            do j = 1, n_elem
               y_tmp(i) = y_tmp(i) +a_elt_tmp(i +n_elem*(j-1))*x_tmp(j)
            enddo
         enddo
         ! --- elemental product in global vector y
         kk = 0
         do k = inc_nn +1, inc_nn +n_elem
            kk = kk +1
            y(eltvar(k)) = y(eltvar(k)) +y_tmp(kk)
         enddo
      enddo
      deallocate(a_elt_tmp, x_tmp, y_tmp)
   return
   endsubroutine prod_elemental_x