PROGRAM MatrMult 
   IMPLICIT NONE
   INCLUDE "mpif.h"

   REAL(8), DIMENSION(:,:), ALLOCATABLE :: A, B, C, tC
   INTEGER :: L, M, N, N0, N1
   REAL ::   tempo2, tempo1, tempro

   INTEGER :: i, j, k, t, st
   REAL(8) :: rmin, rmax

!  MPI variables
   CHARACTER(MPI_MAX_PROCESSOR_NAME) :: server
   INTEGER :: ierr, my_rank, numprocs, ls, dest, tag
   INTEGER status(MPI_STATUS_SIZE)
   LOGICAL :: LetsGo
!
   call MPI_INIT( ierr )
   call MPI_COMM_RANK( MPI_COMM_WORLD, my_rank, ierr )
   call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )
   call MPI_Get_processor_name(server,ls,ierr)
!
   LetsGo = .TRUE.
   IF ( my_rank == 0 ) THEN
      ! Master process looks for input data
      WRITE(*,"(A)") " Input Matrices dimensions L, M, N "
      READ*,L, M, N
      WRITE(*,"(A,3(1x,I))") " Matrices dimensions are: ",L,M,N
      IF ( LetsGo ) THEN
         !  Master allocates entire matrices
         ALLOCATE(A(L,M), B(M,N), C(L,N), STAT=st)
         IF ( st /= 0 ) THEN
   	        WRITE(*,*) "Error allocating matrices "
      	     LetsGo = .FALSE.
         END IF
      END IF
   END IF
   ! Every process synchronizes
   CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
   ! Send start flag
   call MPI_BCAST(LetsGo,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
   IF ( .NOT.LetsGo ) CALL Terminate(my_rank)
   ! Send data to processes
   CALL MPI_BCAST(L,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
   CALL MPI_BCAST(M,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
   CALL MPI_BCAST(N,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
   WRITE(*,*) "Initializing matrices  on server ",TRIM(server),". Input data are: ", &
   &           L, M, N
   IF ( my_rank /= 0 ) THEN
      ALLOCATE(A(L,M), B(M,N), C(L,N), STAT=st)
      IF ( st /= 0 ) THEN
   	 WRITE(*,*) my_rank,": error allocating matrices "
      END IF
   END IF

   ! Every process executes
   A = 0.0; B = 0.0; C = 0.0

   CALL ProcessBand(N,N0,N1,numprocs,my_rank)
   !  Initializations
   DO i = 1, L
      DO k = 1, M
         A(i,k) = REAL(i-k) / REAL(i+k)
      END DO
   END DO

   DO j = 1, N
      DO k = 1, M
         B(k,j) = 0.01 * REAL(k-j) / REAL(k+j)
      END DO
   END DO

   CALL CPU_TIME(tempo1)
   !  Compute multiplication
   print*," Process ",my_rank," computes C(:,",N0,",",N1,")"
   DO j = N0, N1
      DO i = 1, L
         C(i,j) = 0.0
         DO k = 1, M
            C(i,j) = C(i,j) + A(i,k) * B(k,j)
         END DO
      END DO
   END DO
    
   IF ( my_rank /= 0 ) THEN
   !   Send data to master
      dest = 0; tag = 100
      call MPI_Send(C(1,N0), (L*(N1-N0+1)), MPI_DOUBLE_PRECISION, dest, &
                       tag, MPI_COMM_WORLD, ierr)
   END IF
	
   IF ( my_rank == 0 ) THEN
      ALLOCATE(tC(L,N), STAT=st)
      IF ( st /= 0 ) THEN
   	 WRITE(*,*) my_rank,": error allocating matrix tC "
      END IF
      
      ! Master process gathers results
      DO t = 1, numprocs-1
         tag = 100
         CALL ProcessBand(N,N0,N1,numprocs,t)
         call MPI_Recv(C(1,N0), (L*(N1-N0+1)), MPI_DOUBLE_PRECISION, t, &
                       tag, MPI_COMM_WORLD, status, ierr)       
!!          PRINT*,"Process ",my_rank," after MPI_Recv ierr = ",ierr
      END DO
   
      CALL CPU_TIME(tempo2)
   
      tempro =  (tempo2 - tempo1)

      
      WRITE(*,*) "Elapsed time is ",tempro
      	   	   
      tC = MATMUL(A,B)

      rmin = MINVAL(tc-C); rmax = MAXVAL(tC-C)

      WRITE(*,*) " rmin, rmax = ",rmin,rmax
      
   END IF
   CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
   CALL Terminate(my_rank)
   STOP
END PROGRAM MatrMult  

   SUBROUTINE ProcessBand(length,start,end,intervals,n)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: length,intervals,n
      INTEGER, INTENT(OUT) :: start,end

      INTEGER :: delta

      delta = ( length + (intervals - 1) ) / intervals
      start = delta * n + 1; end = delta * ( n + 1)
      if ( end > length ) end = length
   
   END SUBROUTINE ProcessBand

 SUBROUTINE Terminate(p)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: P
    INCLUDE "mpif.h"
    INTEGER :: ierr
 
    WRITE(*,*) "Process ",p," terminates"
    CALL MPI_Finalize(ierr)
 
    STOP
 END SUBROUTINE Terminate