program matvec
!
!*******************************************************************************
!
!! MATVEC uses MPI to multiply a matrix times a vector.
!
!
!  Discussion:
!
!    This is the simple self-scheduling version.
!
!  Reference:
!
!    Gropp, Lusk, Skjellum,
!    Using MPI,
!    Portable Parallel Programming with the Message-Passing Interface,
!    MIT Press, 1997.
!
!    Snir, Otto, Huss-Lederman, Walker, Dongarra,
!    MPI - The Complete Reference,
!    Volume 1, The MPI Core,
!    second edition,
!    MIT Press, 1998.
!
!  Modified:
!
!    02 April 2002
!
  implicit none
!
  include 'mpif.h'

  integer, parameter :: max_cols = 1000
  integer, parameter :: max_rows = 1000
!
  double precision a(max_rows,max_cols)
  double precision ans
  integer anstype
  double precision b(max_cols)
  double precision buffer(max_cols)
  double precision c(max_cols)
  integer cols
  integer donetype
  integer i
  integer ierr
  integer j
  integer job(max_rows)
  integer master
  integer myid
  integer numprocs
  integer numrcvd
  integer numsent
  integer rows
  integer rowtype
  integer sender
  integer status(mpi_status_size)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'MATVEC'
  write ( *, '(a)' ) '  A program using MPI,'
  write ( *, '(a)' ) '  to multiply a matrix times a vector.'
  write ( *, '(a)' ) ' '

  call mpi_init ( ierr )

  call mpi_comm_rank ( mpi_comm_world, myid, ierr )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6,a)' ) 'Process ', myid, ' is alive.'

  call mpi_comm_size ( mpi_comm_world, numprocs, ierr )

  if ( numprocs < 2 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MATVEC - Fatal error!'
    write ( *, '(a)' ) '  Must have at least 2 processes!'
    call mpi_abort ( mpi_comm_world, 1 )
    stop
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) 'Number of processes is ', numprocs

  rowtype = 1
  anstype = 2
  donetype = 3

  master = 0
  rows = 100
  cols = 100
!
!  The master process initializes and then dispatches A and B.
!
  if ( myid == master ) then

    do i = 1, cols
      b(i) = 1
    end do

    do i = 1, cols
      a(i,1:rows) = i
    end do

    numsent = 0
    numrcvd = 0
!
!  Broadcast B to each other process.
!
    call mpi_bcast ( b, cols, mpi_double_precision, master, mpi_comm_world, &
      ierr )
!
!  Send a row of A to each other process.
!
    do i = 1, numprocs-1

      buffer(1:cols) = a(i,1:cols)

      call mpi_send ( buffer, cols, mpi_double_precision, i, rowtype, &
        mpi_comm_world, ierr )

      job(i) = i
      numsent = numsent + 1

    end do
     
    do i = 1, rows

      call mpi_recv ( ans, 1, mpi_double_precision, mpi_any_source, &
        anstype, mpi_comm_world, status, ierr )

      sender = status(mpi_source)
      c(job(sender)) = ans

      if ( numsent < rows ) then

        buffer(1:cols) = a(numsent+1,1:cols)

        call mpi_send ( buffer, cols, mpi_double_precision, sender, &
          rowtype, mpi_comm_world, ierr )

        job(sender) = numsent + 1
        numsent = numsent + 1

      else

        call mpi_send ( 1, 1, mpi_integer, sender, donetype, &
          mpi_comm_world, ierr )

      end if

    end do
! 
!  Print out the answer
!
    do i = 1, cols
      write ( *, '(a,i,a,g14.6)' )  'c(', i, ') = ', c(i)
    end do
!
!  Each slave process receives B, then computes dot products until receiving
!  the "DONE" message.
!
  else

    call mpi_bcast ( b, cols, mpi_double_precision, master, mpi_comm_world, &
      ierr )

    do

      call mpi_recv ( buffer, cols, mpi_double_precision, master, &
        mpi_any_tag, mpi_comm_world, status, ierr )

      if ( status(mpi_tag) == donetype ) then
        exit
      end if 

      ans = dot_product ( buffer(1:cols), b(1:cols) )

      call mpi_send ( ans, 1, mpi_double_precision, master, anstype, &
        mpi_comm_world, ierr )

    end do

  end if
!
!  End of execution.
!
  call mpi_finalize ( ierr )

  stop
end