c     examples/mult/mult_mpi
c     multiplying two matrices, a and b
c
      include "mpif.h"
      parameter (IDIM=20)
      real*8 a(IDIM,IDIM),b(IDIM,IDIM),c(IDIM,IDIM)
      integer my_rank,nproc,ierr,rank
c     
c     INITIALIZE MPI AND DETERMINE BOTH INDIVIDUAL PROCESSOR #
C     AND THE TOTAL NUMBER OF PROCESSORS
c     
      call init(nproc,my_rank)
      if(my_rank.eq.0) then
         call input(n,a,b,IDIM)
      endif
      call MPI_Bcast(n,1,MPI_INTEGER,iroot,MPI_COMM_WORLD,ierr)
      call MPI_Bcast(a,IDIM*n,MPI_DOUBLE_PRECISION,iroot,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Bcast(b,IDIM*n,MPI_DOUBLE_PRECISION,iroot,
     *     MPI_COMM_WORLD,ierr)
      if(my_rank.eq.0) then
         call master(nproc,n,c,IDIM)
      else
         call slave(my_rank,n,a,b,IDIM)
      endif
      if(my_rank.eq.0) then
         do i=1,n
            write(*,*) (c(i,j),j=1,n)
         enddo
      endif
      call MPI_Finalize(ierr)
      stop
      end
c
      subroutine init(nproc,my_rank)
      include "mpif.h"
c
c     initialize the MPI
c
      integer my_rank,nproc
      character*30 nodename
      integer nchar
      call MPI_init(ierr)
      call MPI_Comm_rank(MPI_COMM_WORLD,my_rank,ierr)
      call MPI_Comm_size(MPI_COMM_WORLD,nproc,ierr)
      call MPI_Get_processor_name(nodename,nchar,ierr)
      write(*,*) 'my_rank=',my_rank,' processor=',nodename
      return
      end

      subroutine input(n,a,b,IDIM)
      real*8 a(IDIM,IDIM),b(IDIM,IDIM)
      open(8,file='in.dat',status='unknown')
      read(8,*) n
      do i=1,n
         read(8,*) (a(i,j),j=1,n)
      enddo
      do i=1,n
         read(8,*) (b(i,j),j=1,n)
      enddo
      return
      end
c
      subroutine master(nproc,n,c,IDIM)
      include "mpif.h"
      integer nproc,n,ijob,irecv
      integer status(MPI_STATUS_SIZE)
      real*8 c(IDIM,IDIM),cc
      n_tot=n*n
      nproc1=nproc-1
      ijob=0
      irecv=0
      do while (irecv.lt.n_tot)
        do i=1,nproc1
           if(ijob.lt.n_tot) then
              idest=i
              itag=10
              ijob=ijob+1
              write(*,*) 'master sending ijob=',ijob
              call MPI_Send(ijob,1,MPI_INTEGER,idest,itag,
     *             MPI_COMM_WORLD,ierr)
           endif
        enddo
        do i=1,nproc1
           if(irecv.lt.n_tot) then
              isrc=i
              itag=20
              irecv=irecv+1
              call MPI_Recv(i1,1,MPI_INTEGER,isrc,itag,
     1             MPI_COMM_WORLD,status,ierr)
              call MPI_Recv(j1,1,MPI_INTEGER,isrc,itag,
     1             MPI_COMM_WORLD,status,ierr)
              call MPI_Recv(cc,1,MPI_DOUBLE_PRECISION,isrc,itag,
     1             MPI_COMM_WORLD,status,ierr)
              c(i1,j1)=cc
           endif
        enddo
      enddo
      ijob=-1
      do i=1,nproc1
         idest=i
         itag=10
         call MPI_Send(ijob,1,MPI_INTEGER,idest,itag,
     *        MPI_COMM_WORLD,ierr)
      enddo
      return
      end
c
      subroutine slave(my_rank,n,a,b,IDIM)
      include "mpif.h"
      real*8 a(IDIM,IDIM),b(IDIM,IDIM),cc
      integer ijob,isrc,itag,my_rank
      integer status(MPI_STATUS_SIZE)
      ijob=10
      do while(ijob.gt.0)
         isrc=0
         itag=10
         call MPI_Recv(ijob,1,MPI_INTEGER,isrc,itag,
     1        MPI_COMM_WORLD,status,ierr)
         j=mod(ijob,n)
         if(j.eq.0) then
            i=ijob/n
            j=n
         else
            i=ijob/n+1
         endif
         cc=0.0
         do k=1,n
            cc=cc+a(i,k)*b(k,j)
         enddo
         itag=20
         idest=0
         call MPI_Send(i,1,MPI_INTEGER,idest,itag,
     *        MPI_COMM_WORLD,ierr)
         call MPI_Send(j,1,MPI_INTEGER,idest,itag,
     *        MPI_COMM_WORLD,ierr)
         call MPI_Send(cc,1,MPI_DOUBLE_PRECISION,idest,itag,
     *        MPI_COMM_WORLD,ierr)
      enddo
      return
      end