      subroutine mxmy_mpi(n,a,b,c)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      include "mpif.h"

      dimension istat(mpi_status_size)

      dimension a(n,*),b(n,*),c(n,*)

      data zero/0.0d0/

      call mpi_comm_rank(mpi_comm_world,myid,ierr)
      call mpi_comm_size(mpi_comm_world,numprocs,ierr)

      nprocesses=numprocs-1

      if(myid.eq.0) then
                          isent=0

                          do 1001 islave=0,numprocs-1

                          ia=n-isent
                          ib=numprocs-islave

                          nloop=ia/ib
                          if(mod(ia,ib).ne.0) nloop=nloop+1

                          call mpi_bsend(nloop,1,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)
                          call mpi_bsend(isent,1,mpi_integer,islave,
     &                                   islave,mpi_comm_world,ierr)

                          do 1002 i=1,nloop
 1002                     call mpi_bsend(b(1,isent+i),n,
     &                                   mpi_real8,islave,islave,
     &                                   mpi_comm_world,ierr)

 1001                     isent=isent+nloop
                    endif

      call mpi_recv(nloop,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierr)
      call mpi_recv(isent,1,mpi_integer,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierr)

      do 1003 i=1,nloop
 1003 call mpi_recv(b(1,isent+i),n,mpi_real8,mpi_any_source,
     &              mpi_any_tag,mpi_comm_world,istat,ierr)

      do 1004 j=isent+1,isent+nloop
      do 1005 i=1,n
 1005 c(i,j)=zero
      do 1004 k=1,n
      do 1004 i=1,n
 1004 c(i,j)=c(i,j)+a(i,k)*b(k,j)

      call mpi_bsend(nloop,1,mpi_integer,0,01,mpi_comm_world,ierr)

      do 1006 i=1,nloop
 1006 call mpi_bsend(c(1,isent+i),n,mpi_real8,0,02,mpi_comm_world,ierr)

      if(myid.eq.0) then
                          ireceived=0

                          do 1007 islave=0,numprocs-1

                          call mpi_recv(nloop,1,mpi_integer,islave,
     &                                  01,mpi_comm_world,istat,ierr)

                          do 1008 i=1,nloop
 1008                     call mpi_recv(c(1,ireceived+i),n,
     &                                  mpi_real8,islave,02,
     &                                  mpi_comm_world,istat,ierr)

 1007                     ireceived=ireceived+nloop
                    endif

      call mpi_bcast(c,n*n,mpi_real8,0,mpi_comm_world,ierr)

      return
      end
