      program macwek
c
c Prosty program sprawdzajacy dzialanie rozproszonej procedury mnozenia
c macierzy A przez wektor V z wykorzystaniem procedury MPI_REDUCE_SCATTER.
c Kazdy procesor wykorzystuje "swoj" pasek macierzy A (macierz A_LOC) 
c oraz "swoja czesc wektora V (wektor V_LOC) a wyniki sa przypisywane do 
c wektora W, ktory jest rownomiernie rozproszony pomiedzy procesory (kazdy 
c posiada czesc W_LOC). Poniewaz jest to program przykladowy, kazdy procesor i
c sam tworzy cala macierz A i wektor V, potem uzywajac tylko ich czesci;
c w rzeczywistych wykorzystaniach mnozenia bardzo duzych tablic kazdy procesor
c bedzie faktycznie posiadal tylko swoje czesci tablic.
c
      implicit none
      include 'mpif.h'
      integer max_n,max_m,max_procs
      parameter (max_n=100,max_m=100,max_procs=16)
      integer n,m,n_loc,m_loc
      integer i,j,ii
      integer me,nprocs,ierror
      real a(max_m,max_n),v(max_n),w(max_m),wi
      real a_loc(max_m,max_n),v_loc(max_n),w_loc(max_m)
      integer counts(0:Max_Procs-1)
c
      call MPI_Init( IERROR )
      call MPI_Comm_Rank( MPI_COMM_WORLD, me, IERROR )
      call MPI_Comm_Size( MPI_COMM_WORLD, nprocs, IERROR )
c
c Obliczanie elementow A i V; tylko tak, aby wypelnic te tablice bez wiekszego
c zastanawiania sie a jednoczesnie w sposob niezbyt trywialny zeby mozna bylo
c program przetestowac dla przypadku jak najogolniejszego.
c
      n=2*nprocs
      m=3*nprocs
      write (*,*) "Jestem procesorem nr",me," z",nprocs 
      do i=1,m
        do j=1,n
          a(i,j)=cos(float(i+j))
        enddo
      enddo
      do i=1,n
        v(i)=sin(float(i))
      enddo

      if (me .eq. 0) then
c
c Kontrolne obliczanie A*V w trybie szeregowym
c
      do i=1,m
        wi=0.0
        do j=1,n
          wi=wi+a(i,j)*v(j)
        enddo
        w(i)=wi
      enddo
      write (*,*) "Wektor w - obliczenie szeregowe" 
      write (*,'(8f8.3)') (w(i),i=1,m)

      endif ! me .eq. 0
c
c Wypelnianie a_temp i v_temp dla kazdego procesora
c
      n_loc=n/nprocs
      m_loc=m/nprocs
      do i=me*n_loc+1,(me+1)*n_loc
        ii=i-me*n_loc
        v_loc(ii)=v(i)
        do j=1,m
          a_loc(j,ii)=a(j,i)
        enddo
      enddo
      do i=0,nprocs-1
        counts(i)=m_loc
      enddo 
c
c Drukowanie kawalkow macierzy A i wektora V odpowiadajacych poszczegolnym
c procesorom
c
      write (*,'(a,i3,a)') 
     &  "Processor",me,"> Macierz A_LOC i wektor V_LOC"
      write (*,'(a,i3,a,3x,20i8)') 
     &  "Processor",me,"> ",(j+me*n_loc,j=1,n_loc)
      do i=1,m
        write (*,'(a,i3,a,i3,20f8.3)') 
     &  "Processor",me,"> ",i,(a_loc(i,j),j=1,n_loc)
      enddo
      write (*,*)
      write (*,'(a,i3,a,3x,20f8.3)') 
     &  "Processor",me,"> ",(v_loc(j),j=1,n_loc)
c
c Rozproszone obliczanie iloczynu A*V. Nalezy zwrocic uwage, ze podprogram
c jest wolana z M a nie M_LOC (zob. komentarz w tekscie podprogramu).
c
      call matvec(n_loc, m, a_loc, v_loc, w_loc, counts, MPI_COMM_WORLD)
c
c Drukowanie kawalkow iloczyny A*V dla poszczegolnych procesorow
c
      write (*,'(a,i3,a)') 
     &  "Processor",me,"> Wektor w"
      do i=1,m_loc
        write (*,'(a,i3,a,i3,f10.5)') 
     &  "Processor",me,"> ",i+me*m_loc,w_loc(i)
      enddo
      call MPI_Finalize ( IERROR ) 
      stop
      end
c------------------------------------------------------------------------------
      subroutine matvec( n, m, lmatrix, lx, ly, counts, comm )
c
c Procedura rozproszonego mnozenia macierzy przez wektor z rozproszeniem
c elementow wynikowego wektora pomiedzy procesory. Znaczenie parametrow
c jest nastepujace:
c
c M - faktyczny wymiar wierszowy macierzy oraz wymiar wektora wynikowego,
c N - liczba kolumn macierzy oraz elementow wektora wejsciowego posiadanych
c     przez procesor,
c LMATRIX - "pasek" macierzy wejsciowej ("pocietej" wzdluz kolumn) posiadany
c     przez procesor,
c LX - czesc wektora wejsciowego posiadanego przez procesor,
c LY - czesc wektora wynikowego posiadanego przez procesor,
c COUNTS - COUNTS(ME) jest liczba elementow wektora wynikowego posiadanych
c      przez procesor,
c COMM - komunikator.
c 
      implicit none
      include 'mpif.h'
      integer comm,ierr
      integer max_n,max_m,max_procs
      parameter (max_n=100,max_m=100,max_procs=16)
      integer n,m,counts(0:max_procs-1)
      integer i,j
c
c Dla nie-Fortranowcow: faktyczny wymiar wierszowy tablicy LMATRIX wynosi M
c ale z uwagi na leksykograficzne przesylanie macierzy musimy dac maksymalny
c wymiar wierszowy.
c
      real lmatrix(max_m,n),lx(n),ly(m)
      real sum,tmp(max_m)
c
c Wykonywanie przez procesor przypisanej mu czesci obliczania wektora W. Wynik
c jest przechowywany w tablicy TMP. Nalezy zwrocic uwage, ze tablica TMP
c ma dla kazdego procesora M a nie M_LOC elementow, natomiast dla danego
c procesora jej elementy sa "czesciami" elementow wynikowego wektora W.
c
      do i=1,m
         sum = 0
         do j=1,n
            sum = sum + lmatrix(i,j)*lx(j)
         enddo
         tmp(i) = sum
      enddo
c
c Sumowanie skladnikow elementow wektora W obliczonych przez poszczegolne
c procesory i rozpraszanie wynikow pomiedzy procesory; kazdy z nich 
c dysponuje po tej operacji czescia wynikowego wektora zapamietana w 
c wektorze LY o wymiarze COUNTS(ME), gdzie ME jest rzedem danego procesora.
c
      call MPI_REDUCE_SCATTER( tmp, ly, counts, 
     &  MPI_REAL, MPI_SUM, comm, ierr)
      return
      end
