cccccccccccccccccccccccccccccccccccccccccccccccc
c      Matrix Multiplication MPI Program       c
c      For this simple version, # of procssors c
c      equals # of columns in matrix           c
c                                              c
c      To run, mpirun -np 4 a.out              c
cccccccccccccccccccccccccccccccccccccccccccccccc

      include 'mpif.h'

      parameter (ncols=4, nrows=4)
      integer a(ncols,nrows), b(ncols,nrows), c(ncols,nrows)
      integer buf(ncols),ans(nrows)
      integer myid, root, numprocs, ierr, status(MPI_STATUS_SIZE)
      integer sender, count

      call MPI_INIT(ierr)
      call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
      call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )

      if(numprocs.ne.4) then
        print *, "Please run this exercise on 4 processors"
        call MPI_FINALIZE(ierr)            
        stop
      endif

      root = 0
      tag = 100
      count = nrows*ncols

c     Master initializes and then dispatches to others
      IF ( myid .eq. root ) then

         do j=1,ncols
            do i=1,nrows
               a(i,j) = 1
               b(i,j) = j
            enddo
         enddo

c        send a to all other process
         call MPI_BCAST(a,count,MPI_INTEGER,root,MPI_COMM_WORLD,ierr)

c        send one column of b to each other process
         do j=1,numprocs-1
            do i = 1,nrows
               buf(i) = b(i,j+1)
            enddo
         call MPI_SEND(buf,nrows,MPI_INTEGER,j,tag,MPI_COMM_WORLD,ierr)
         enddo
         
c        Master does his own part here
         do i=1,nrows
            ans(i) = 0
            do j=1,ncols
               ans(i) = ans(i) + a(i,j) * b(i,1)
            enddo
         c(i,1) = ans(i)
         enddo

c        then receives answers from others

         do j=1,numprocs-1
            call MPI_RECV(ans, nrows, MPI_INTEGER, MPI_ANY_SOURCE,
     $           MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

            sender = status(MPI_SOURCE)
            do i=1,nrows
               c(i,sender+1) = ans(i)
            enddo

         enddo
         
         do i=1,nrows
            write(6,*)(c(i,j),j=1,ncols)
         enddo

      ELSE

c        slaves receive a, and one column of b, then compute dot product
         call MPI_BCAST(a,count,MPI_INTEGER,root,MPI_COMM_WORLD,ierr)

         call MPI_RECV(buf, nrows, MPI_INTEGER, root,
     $        MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

         do i=1,nrows
            ans(i) = 0
            do j=1,ncols
               ans(i) = ans(i) + a(i,j) * buf(j)
            enddo
         enddo

         call MPI_SEND(ans,nrows,MPI_INTEGER,root,0,MPI_COMM_WORLD,ierr)

      ENDIF

      call MPI_FINALIZE(ierr)

      stop
      end