      program prog
!
      use numerics
      implicit none
!
      integer,  parameter ::  n = 3
      real(rp)  ::  A(n,n), x(n), y(n)
      real(rp)  ::  sum, xydot
      integer   ::  i, j
!
!     initialize arrays
      do i = 1, n
         x(i) = 2.0_rp
         y(i) = real(i, rp)
      end do
      do j = 1, n
         do i = 1, n
            A(i,j) = real(i*j - i, rp)
         enddo
      enddo

!     compute the dot product of x and y
      call dot(x,y,n,xydot)
      write (6, '(A,3F9.6)') 'dot product of x & y =  ', xydot

!     compute y = Ax
      do i = 1, n
         y(i) = 0.0_rp
      end do
      do j = 1, n
         do i = 1, n
            y(i) = y(i) + A(i,j)*x(j)
         enddo
      enddo
      write (6, '(A,3F10.6)')  'y =', y
!
      stop
      end program prog

