mxdsmm Subroutine

      pure subroutine mxdsmm(n,a,x,y)

      integer,intent(in) :: n      !! order of the matrix a.
      real(wp),intent(in) :: a(*)  !! `a(n*(n+1)/2)`  dense symmetric matrix stored in the packed form.
      real(wp),intent(in) :: x(*)  !! x(n)  input vector.
      real(wp),intent(out) :: y(*)  !! y(n)  output vector equal to ` a*x`.

      real(wp) :: temp
      integer :: i , j , k , l

      k = 0
      do i = 1 , n
         temp = 0.0_wp
         l = k
         do j = 1 , i
            l = l + 1
            temp = temp + a(l)*x(j)
         end do
         do j = i + 1 , n
            l = l + j - 1
            temp = temp + a(l)*x(j)
         end do
         y(i) = temp
         k = k + i
      end do

      end subroutine mxdsmm