subroutine matvec(n, m, lmatrix, lx, ly, counts, comm) use mpi integer n, m, comm, counts(*) real lmatrix(n,m), lx(m), ly(m) integer i, j real sum real, allocatable :: tmp(:) allocate (tmp(n)) ! Perform the local matrix-vector multiply ! 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 tmp(i) = sum enddo ! Perform the local matrix-vector product call MPI_REDUCE_SCATTER(tmp, ly, counts, & MPI_REAL, MPI_SUM, comm, ierr) deallocate (tmp) ! We're done! end