subroutine matvec( n, m, lmatrix, lx, ly, counts, comm, work ) 
      include 'mpif.h' 
      integer n, m, comm, counts(*) 
      real lmatrix(n,m), lx(m), ly(m) 
      integer i, j 
      real sum 
      real work(n) 
C 
C Perform the local matrix-vector multiply 
C Should use the level-2 BLAS routine SGEMV 
      do i=1,n 
         sum = 0 
         do j=1,m 
            sum = sum + lmatrix(i,j)*lx(j) 
         enddo 
         work(i) = sum 
      enddo 
C 
C Perform the local matrix-vector product 
      call MPI_REDUCE_SCATTER( work, ly, counts,  
     *                         MPI_REAL, MPI_SUM, comm, ierr) 
C 
C We're done! 
      end