program matrix
	    implicit none
	    integer, parameter :: maxsize = 500
	    integer, parameter :: maxnum = 100
	    integer, parameter :: rows1 = 2, cols1 = 4
	    integer, parameter :: rows2 = 4, cols2 = 3
	    integer m(rows1, cols1)
	    integer n(rows2, cols2)
	    integer prod(rows1, cols2)
	    integer prod2(rows1, cols2)
	    real, dimension(5) :: randomarray
	    integer :: counter = 1

	    
	    call random_seed()
	
	    write(*,*) 'Checking matrix multiplication algorithm...'
	    call makeTestMatrix(m, rows1, cols1, 1)
	    call makeTestMatrix(n, rows2, cols2, 1)
	
	    write(*,*) 'Matrix m:'
	    call printMatrix(m, rows1, cols1)
	    write(*,*) 'Matrix n:'
	    call printMatrix(n, rows2, cols2)
	    write(*,*) 'Multiplying matrices'
	    call matMult(m,n,prod, rows1, cols1, rows2, cols2)
	
	    write(*,*) 'Product matrix:'
	    call printmatrix(prod, rows1, cols2)
	    
	    write(*,*) '2nd product matrix:'
	    prod2 = matmul(m,n)
	    call printmatrix(prod2, rows1, cols2)
	    
	contains 
	    subroutine makeVector (v)
	        implicit none
		integer, dimension (maxsize), intent(out) :: v
		integer i
		
		do i = 1, maxsize
		    call random_number(randomarray)
		    v(i) = int(randomarray(1)*maxnum)
		end do
	    end subroutine
	    
	    subroutine printVector(v, whichone)
	        implicit none
		integer, dimension(maxsize), intent(in) :: v
		integer, intent(in) :: whichone
		integer i
		write (*,*) 'V', whichone, ':', (v(i),i=1,maxsize)
	     end subroutine
	     
	     subroutine makeTestMatrix(m, rows, cols, start)
	        implicit none
		integer, dimension(rows,cols), intent(out) :: m
		integer, intent(in) :: rows, cols
		integer, intent(in) :: start
		integer i, j, counter
		
		counter = start
		do i=1, rows
		    do j = 1, cols
		        m(i,j) = counter
			counter = counter + 1
		    end do
		end do
	     end subroutine
	     
	     subroutine makeMatrix(m, rows, cols)
	        implicit none
		integer, dimension(rows,cols), intent(out) :: m
		integer, intent(in) :: rows, cols
		integer i, j
		
		do i=1,rows
		    do j = 1, cols
		        call random_number(randomarray)
			m(i,j) = int(randomarray(1)*maxnum)
		    end do
		end do
	     end subroutine
		
	     subroutine printMatrix(m, rows, cols)
	        implicit none
		integer, dimension(rows,cols), intent(out) :: m
		integer, intent(in) :: rows, cols
		integer i, j	
		
		do i=1,rows
		    write(*,*) (m(i,j), j = 1,cols)		
		end do
	     end subroutine	
	     
	     subroutine matMult(m, n, prod, rows1, cols1, rows2, cols2)
	        implicit none
		integer, dimension(rows1,cols1), intent(in) :: m
		integer, dimension(rows2,cols2), intent(in) :: n
		integer, dimension(rows1,cols2), intent(out) :: prod
		integer, intent(in) :: rows1, cols1, rows2, cols2
		integer i, j, k
	        do i=1, rows1 
		    do j=1, cols2
		        prod(i,j)=0
	                do k = 1, cols1 ! or rows2: cols1 = rows2 in matrix multiplication
		            prod(i,j) = prod(i,j) + m(i,k)*n(k,j)
		        end do
	            end do
	        end do
	     end subroutine
	
	END PROGRAM