program Matmul

  implicit none
  save
  include 'mpif.h'

  character(LEN=50) :: stringa
  
  integer, parameter :: N = 8
  integer :: nprocs, ierr, proc_me, status(MPI_STATUS_SIZE)
  integer :: ncount, i, j, iglob, k, z
  real, allocatable, dimension(:,:) :: A, B, C, buffer
  real, allocatable, dimension(:) :: sup
  integer :: Nrow, Ncol
  
  call MPI_INIT(ierr)
  call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr)
  call MPI_COMM_RANK(MPI_COMM_WORLD,proc_me,ierr)

  Nrow = N / nprocs
  NCol = N

  allocate(A(Nrow,Ncol))
  allocate(B(Nrow,Ncol))
  allocate(C(Nrow,Ncol))
  allocate(buffer(Nrow,Ncol))
  allocate(sup(Ncol))

  !---------- inizializzazione delle matrici A e B 
  do j=1,Ncol
     do i=1,Nrow
        A(i,j) = 1.0/j
! ((Nrow * proc_me)+i) * j
        B(i,j) = 2.0/ A(i,j)
        C = 0.0
     enddo
  enddo
  
  write(stringa, *) Ncol



  do ncount=0, nprocs-1
     call MPI_Allgather(B(:, (ncount*Nrow)+1:ncount*Nrow+Nrow), NRow*Nrow, MPI_REAL, buffer, &
                                                     Nrow*Nrow, MPI_REAL, MPI_COMM_WORLD, ierr)
     do k = (ncount*Nrow) + 1, (ncount+1) * Nrow
        do j = 1, Ncol 
           do i = 1, Nrow
              C(i,k) = C(i,k) + A(i,j) * buffer(mod(j-1,NRow)+1,((j/(Nrow+1)*Nrow))+k-(ncount*Nrow))
           end do
        end do
    end do
  end do
  
  if (proc_me == 0) then

! print matrix C
    print*, "Matrix C"
      
    do i=1,Nrow
       print '('//trim(stringa)//'(F8.2))', C(i,:)
    enddo
     
     do k=1, nprocs-1
        call MPI_RECV(C, Nrow*Ncol, MPI_REAL, k, 0, MPI_COMM_WORLD, status, ierr)
        do i=1,Nrow
           print '('//trim(stringa)//'(F8.2))', C(i,:)
        enddo
     enddo
  else
     call MPI_SEND(C, Nrow*Ncol, MPI_REAL, 0, 0, MPI_COMM_WORLD, ierr)
  end if

  call MPI_FINALIZE(ierr)
  
end program Matmul
