program main
     implicit none

     ! Wymiary macierzy: A:x,y   B:y,z   C:x,z <- macierz wynikowa
     !Jak macierze nie są kwadratowe, to 3 sposób daje złe wyniki :(, ale luz bo i tak pomiary czasu są na kwadratowych macierzach
     integer, parameter :: x = 4
     integer, parameter :: y = 4
     integer, parameter :: z = 4

     ! Zmienne pomocnicze do iteracji
     integer :: i,j,k

     ! Zmienna pomocnicza tmp, przy generowaniu wartości losowych
     real(kind=4) :: tmp

     ! Macierze na których bede pracować:
     real :: A,B,C
     dimension A(x,y), B(y,z), C(x,z)

     ! Wektory pomocnicze, przy iloczynie skalarnym (3 sposób mnożenia)
     real:: vec_A(y)
     real:: vec_B(y)

     ! Wypełniam macierz A i B losowymi wartościami, macierz C zerami.
     do i=1,x
        do j=1,y
            CALL RANDOM_NUMBER(tmp)
            A(i,j)=tmp
        end do
     end do

     do i=1,y
        do j=1,z
            CALL RANDOM_NUMBER(tmp)
            B(i,j)=tmp
        end do
     end do

     do i=1,x
        do j=1,z
            C(i,j)=0
        end do
     end do
     write(*,*) "Macierz A:"
     print *,A
     write(*,*) "Macierz B:"
     print *,B
     write(*,*) "Macierz C:"
     print *,C

     ! Obliczam macierz C pierwszym sposobem: (naive math)
     do i=1,x
        do j=1,y
            do k=1,z
                C(i,j)=C(i,j)+A(i,k)*B(k,j)
            end do
        end do
     end do

     write(*,*) "Macierz C - wynik A x B (naive math):"
     print *,C
     do i=1,x
        do j=1,z
            C(i,j)=0
        end do
     end do

     ! Obliczam macierz C drugim sposobem: (better math)

     do j=1,y
        do k=1,z
            do i=1,x
                C(i,j)=C(i,j)+A(i,k)*B(k,j)
            end do
        end do
     end do

     write(*,*) "Macierz C - wynik A x B (better math):"
     print *,C
     do i=1,x
        do j=1,z
            C(i,j)=0
        end do
     end do

     ! Obliczam macierz C trzecim sposobem: (better math + dot_product)

    do k=1,z
        do  i=1,x
            vec_A= A(i,:)
            vec_B=B(:,k)
            C(i,k)=dot_product(vec_A,vec_B)
        end do
    end do

     write(*,*) "Macierz C - wynik A x B (better math + dot_product):"
     print *,C
     do i=1,x
        do j=1,z
            C(i,j)=0
        end do
     end do


     ! Obliczam macierz C czwartym sposobem: (MATMUL)
     C = matmul(A,B)
     write(*,*) "Macierz C - wynik A x B (matmul):"
     print *,C



end program main