Program Matrices_Multiplication
 

C... Example : matrices multiplication.
C... Each CPU calculates separate columns of matrix
      implicit real*8 (a-h,o-z)
      parameter (n=400)  ! matrix dimension
      common/matrix/ a(n,n), b(n,n), c(n,n)
C...   MPI initialization.  nProc=<Total Number of Processes>,  
C...                        myProcess=<Current Process Number>
      include 'mpif.h'
      integer status(MPI_STATUS_SIZE)

      call MPI_Init(ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD,myProcess,ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD,nProcs,ierr)
      nc=n/nProcs            ! columns distribution (we consider n>nProcs)
      nrest=mod(n,nProcs)
      if(myProcess.lt.nrest) then
        nc1=1+(nc+1)*myProcess
        nc2=nc1+nc
      else
        nc1=1+(nc+1)*nrest+nc*(myProcess-nrest)
        nc2=nc1+nc-1
      endif
      write(*,*) ' Process',myProcess,'  of',nProcs,
     -           '  started for columns',nc1,nc2
C...
      if (myProcess.eq.0) then ! process 0 sends matrix to others
        do i=1,nProcs-1
          call MPI_Send(a,2*n*n,MPI_DOUBLE_PRECISION,i,itag,
     -                  MPI_COMM_WORLD,ierr)   ! Both A and B matrices !!!
        enddo
      else    ! other processes receive matrix from process 0
        call MPI_RECV(a,2*n*n,MPI_DOUBLE_PRECISION,0,itag,
     -                  MPI_COMM_WORLD,status,ierr)
      endif
C--- start multiplication
        do i=1,n
          do j=nc1,nc2      ! only own columns !!! 
            s=0
            do k=1,n
              s=s+a(i,k)*b(k,j)
            enddo
            c(i,j)=s
          enddo
        enddo
C--- finish multiplication
      if(myProcess.eq.0) then
        do i=1,nProcs-1
          if(i.lt.nrest) then
            nc1=1+(nc+1)*i
            k=nc+1
          else
            nc1=1+(nc+1)*nrest+nc*(i-nrest)
            k=nc
          endif
          call MPI_RECV(c(1,nc1),n*k,MPI_DOUBLE_Precision,i,itag,
     -                  MPI_COMM_WORLD,status,ierr)
        enddo
      else
          call MPI_SEND(c(1,nc1),n*(nc2-nc1+1),MPI_DOUBLE_Precision,
     -                  0,itag,MPI_COMM_WORLD,ierr)
      endif
      call MPI_Finalize(ierr)
      end