! matvec_ex.f90
include "mpiSim.f90"
program main
include "useSim.f90"

! This is a sample MPI application that illustrates how to use the
! MPI library simulation package. To compile this example with
! Lahey(with run-time debugging):
! 1) put USE and INCLUDE statements in the main program, matvec_jb.f90
! 2) issue(from Hammer for example):
!    lf95 matvec_jb.f90 --chkglobal -g --co --nsav --pca --warn --staticlink

!*********************************************************************72
!  Code from: http://people.scs.fsu.edu/~burkardt/f_src/mpi/mpi.html
!  MAIN is the main program for MATVEC.
!
!  Discussion:
!
!    MATVEC uses MPI to compute a matrix-vector product b = A * x.
!
!    This is the simple self-scheduling version.  Each worker is given a copy
!    of x, and then is fed one row of A.  As soon as it computes
!    B(I) = A(I,1:N)*x(1:N), it is given another column of A, unless there are
!    no more, in which case it is sent a "terminate" message.  Thus, a faster
!    process will be given more work to do.
!
!    By using allocatable arrays, the amount of memory used has been controlled.
!    The master process allocates A and x, but the worker processes only
!    allocate enough memory for one row of A, and x.
!
!  Modified:
!
!    10 October 2002
!
!  Reference:
!
!    William Gropp, Ewing Lusk, Anthony Skjellum,
!    Using MPI: Portable Parallel Programming with the
!    Message-Passing Interface,
!    Second Edition,
!    MIT Press, 1999,
!    ISBN: 0262571323.
!
!    Snir, Otto, Huss-Lederman, Walker, Dongarra,
!    MPI - The Complete Reference,
!    Volume 1, The MPI Core,
!    second edition,
!    MIT Press, 1998.
!

!
!  Fortran77 include file:
!
!  include 'mpif.h'
!
!  Fortran90 module:
!
!  use mpi
!
! implicit none
!
  double precision, allocatable, dimension ( :, : ) :: a
  double precision, allocatable, dimension ( : ) :: a_row
  double precision ans
  double precision, allocatable, dimension ( : ) :: b
  integer dest
  integer dummy
  integer i
  integer ierr
  integer j
  integer j_one
  integer m
  integer, parameter :: master = 0
  integer my_id
  integer n
  integer num_procs
  integer num_rows
  integer num_workers
  double precision, parameter :: pi = 3.141592653589793D+00
  integer status(MPI_STATUS_SIZE)
  integer tag
  integer, parameter :: tag_done = -1
  double precision, allocatable, dimension ( : ) :: x
!
!  Initialize MPI.
!
  call MPI_Init ( ierr )
!
!  Get this processor's ID.
!
  call MPI_Comm_rank ( MPI_COMM_WORLD, my_id, ierr )
!
!  Get the number of processors.
!
  call MPI_Comm_size ( MPI_COMM_WORLD, num_procs, ierr )

  if ( my_id == 0 ) then

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MATVEC - Master process:'
    write ( *, '(a)' ) '  FORTRAN90 version'
    write ( *, '(a)' ) '  An MPI example program to compute'
    write ( *, '(a)' ) '  a matrix-vector product b = A * x.'
    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  The number of processes is ', num_procs

  end if
  write ( *, *) "my_id=", my_id
  write ( *, '(a)' ) ' '
  write ( *, '(a,i8,a)' ) '  Process ', my_id, ' is active.'
!
!  Let the size of the array be determined by process 0.
!  The other processes need to know N, but not M!
!
  if ( my_id == 0 ) then
    m = 100
    n = 50
    write ( *, '(a)' ) ' '
    write ( *, '(a,i8)' ) '  The number of rows is    ', m
    write ( *, '(a,i8)' ) '  The number of columns is ', n
  end if

  call MPI_Bcast ( n, 1, MPI_INTEGER, master, MPI_COMM_WORLD, ierr )
!
!  The master process allocates and initializes A and X.
!
  write(*,*) "my_id, master=",my_id,master
  if ( my_id == master ) then

    allocate ( a(1:m,1:n) )
    allocate ( x(1:n) )
    allocate ( b(1:m) )
    b=0
    do i = 1, m
      do j = 1, n
        a(i,j) = sqrt ( 2.0D+00 / dble ( n + 1 ) ) &
          * sin ( dble ( i * j ) * pi / dble ( n + 1 ) )
      end do
    end do
!
!  X is specially chosen so that b = A * x is known in advance.
!  The value of B will be zero, except that entry J_ONE will be 1.
!  Pick any value of J_ONE between 1 and M.
!
    j_one = 17

    do i = 1, n
      x(i) = sqrt ( 2.0D+00 / dble ( n + 1 ) ) &
        * sin ( dble ( i * j_one ) * pi / dble ( n + 1 ) )
    end do
!
!  Worker processes set aside room for one row of A, and for the vector X.
!
  else

    allocate ( a_row(1:n) )
    allocate ( x(1:n) )

  end if
!
!  Process 0 broadcasts the vector X to the other processes.
!
  call MPI_Bcast ( x, n, MPI_DOUBLE_PRECISION, master, MPI_COMM_WORLD, ierr )

  if ( my_id == master ) then
!
!  Process 0 sends one row of A to all the other processes.
!
!  Note that the call to MPI_Send uses a FORTRAN90 array section.  Even
!  though the elements of a 2D array row are not contiguous as stored in memory,
!  FORTRAN90 interprets the expression "A(I,1:N)" as requiring it to make
!  a temporary, and contiguous, copy of the indicated elements.
!
    num_rows = 0

    do i = 1, num_procs-1

      num_rows = num_rows + 1
      dest = i
      tag = num_rows

      call MPI_Send ( a(num_rows,1:n), n, MPI_DOUBLE_PRECISION, dest, tag, &
        MPI_COMM_WORLD, ierr )

    end do

    num_workers = num_procs - 1

    do

      call MPI_Recv ( ans, 1, MPI_DOUBLE_PRECISION, MPI_ANY_SOURCE, &
        MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr )

      tag = status(MPI_TAG)
      b(tag) = ans

      if ( num_rows < m ) then

        num_rows = num_rows + 1
        dest = status(MPI_SOURCE)
        tag = num_rows

        call MPI_Send ( a(num_rows,1:n), n, MPI_DOUBLE_PRECISION, dest, &
          tag, MPI_COMM_WORLD, ierr )

      else

        num_workers = num_workers - 1
        dummy = 0
        dest = status(MPI_SOURCE)
        tag = tag_done

        call MPI_Send ( dummy, 1, MPI_INTEGER, dest, tag, MPI_COMM_WORLD, ierr )

        if ( num_workers <= 0 ) then
          exit
        end if

      end if

    end do

    deallocate ( a )
    deallocate ( x )
!
!  Each worker process repeatedly receives rows of A (with TAG indicating
!  which row it is), computes dot products A(I,1:N) * X(1:N) and returns
!  the result (and TAG), until receiving the "DONE" message.
!
  else

    do

      call MPI_Recv ( a_row, n, MPI_DOUBLE_PRECISION, master, &
        MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr )

      tag = status(MPI_TAG)

      if ( tag == tag_done ) then
        exit
      end if

      ans = dot_product ( a_row(1:n), x(1:n) )

      call MPI_Send ( ans, 1, MPI_DOUBLE_PRECISION, master, tag, &
        MPI_COMM_WORLD, ierr )

    end do

    deallocate ( a_row )
    deallocate ( x )

  end if
!
!  Print out the answer.
!
  if ( my_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MATVEC - Master process:'
    write ( *, '(a)' ) '  Product vector b = A * x'
    write ( *, '(a,i8)' ) '  (Should be zero, except for a 1 in entry ', j_one
    write ( *, '(a)' ) ' '
    do i = 1, m
      write ( *, '(i8,g14.6)' )  i, b(i)
    end do

    deallocate ( b )

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

  if ( my_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MATVEC - Master process:'
    write ( *, '(a)' ) '  Normal end of execution.'
  end if

  stop
end

   deallocate ( b )

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

  if ( my_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'MATVEC - Master process:'
    write ( *, '(a)' ) '  Normal end of execution.'
  end if

  stop
end