! MPI Tutorial
! Dr. Andrew C. Pineda, HPCERC/AHPCC
! Dr. Brian Smith, HPCERC/AHPCC
! The University of New Mexico
! November 17, 1997
! Last Revised: September 18, 1998

program matvec2
  ! Perform matrix vector product -- Y = AX
  ! This is method two -- distribute A by block columns
  ! and X in blocks (of rows) and the partial vector sum of Y is on
  ! each processor.
  include 'mpif.h'
  integer, parameter :: dim1 = 80, dim2 = 10, dim3 = dim1*dim2
  integer ierr, rank, size, root, i, j
  integer sec_start, nano_start
  integer sec_curr, nano_curr
  integer sec_startup, nano_startup
  integer sec_comp, nano_comp
  integer sec_cleanup, nano_cleanup
  real, dimension(dim1,dim1) :: a
  real, dimension(dim1,dim2) :: apart
  real, dimension(dim1) :: x, y, ypart
  real, dimension(dim2) :: xpart
  interface
     subroutine posix_timer(job_sec, job_nanosec)
       integer job_sec, job_nanosec
     end subroutine posix_timer
  end interface
  root = 0
  call MPI_INIT( ierr )
  call MPI_COMM_RANK( MPI_COMM_WORLD, rank, ierr )
  call MPI_COMM_SIZE( MPI_COMM_WORLD, size, ierr )
  print *, 'START process on processor ', rank
  if( rank == root ) then
     call posix_timer(sec_start, nano_start)
     do i = 1, dim1
        x(i) = 1.0
        do j = 1, dim1
           a(j,i) = i + j
        enddo
     enddo
  endif
  !  Distribute the 80x80 array A by columns as
  ! 80x10 blocks stored in APART.
  ! Distribute the 80 dimensional
  ! array x in blocks of length 10
  call MPI_SCATTER( a, dim3, MPI_REAL, apart, dim3, MPI_REAL, root,&
       & MPI_COMM_WORLD, ierr )
  call MPI_SCATTER( x, dim2, MPI_REAL, xpart, dim2, MPI_REAL, root,&
       & MPI_COMM_WORLD, ierr )
  if( rank == root ) then
     call posix_timer(sec_curr, nano_curr)
     sec_startup = sec_curr - sec_start
     nano_startup = nano_curr - nano_start
     sec_start = sec_curr
     nano_start = nano_curr
  endif
  do j = 1, dim1
     ypart(j) = 0.0
  enddo
  do i = 1, dim2
     do j = 1, dim1
        ypart(j) = ypart(j) + xpart(i)*apart(j,i)
     enddo
  enddo
  if( rank == root ) then
     call posix_timer(sec_curr, nano_curr)
     sec_comp = sec_curr - sec_start
     nano_comp = nano_curr - nano_start
     sec_start = sec_curr
     nano_start = nano_curr
  endif
  call MPI_REDUCE( ypart, y, dim1, MPI_REAL, MPI_SUM, root, &
       MPI_COMM_WORLD, ierr )
  if( rank == root ) then
     call posix_timer(sec_curr, nano_curr)
     sec_cleanup = sec_curr - sec_start
     nano_cleanup = nano_curr - nano_start
  endif
  print *, 'Finish processor ', rank
  if( rank == root ) then
     print *, 'Matrix vector product, elements 10 and 60, are: ',&
          & y(10), y(60)
     print *, 'Startup execution times (sec, nano): ',&
          & sec_startup, nano_startup
     print *, 'Computation execution times (sec, nano): ',&
          & sec_comp, nano_comp
     print *, 'Cleanup execution times (sec, nano): ',&
          & sec_cleanup, nano_cleanup
  endif
  call MPI_FINALIZE( ierr )
end program matvec2
