MatrixMultiplicationAndTranspose

dimension A(10,10),B(10,10),C(10,10),Btranspose(10,10),AB(10,10),ABBtranspose(10,10),ABBtransposeB(10,10)
    arow=3
    acol=3
    brow=3
    bcol=2
    open(5,file='A.txt')
    do i=1,arow
        read (5,*)(A(i,j),j=1,acol)
    end do
    open(5,file='B.txt')
    do i=1,brow
        read (5,*)(B(i,j),j=1,bcol)
    end do

    !Do [A][B][B]'[B]
    call MatMul(A,B,AB,arow,acol,brow,bcol)
    call Transpose(B,Btranspose,brow,bcol)
    call MatMul(AB,Btranspose,ABBtranspose,arow,bcol,bcol,brow)
    call MatMul(ABBtranspose,B,ABBtransposeB,arow,brow,brow,bcol)

    do i=1,arow
        write (*,*)(ABBtransposeB(i,j),j=1,bcol)
    end do
    end program


subroutine Transpose(A,R,rows,cols)
dimension A(10,10),R(10,10) !CAUTION: This declaration is crucial
    do i=1,cols !i goes from 1 to cols
        do j=1,rows !and j from 1 to rows because transpose should have the dimension: cols x rows
            R(i,j)=A(j,i) !fills from column of A to row of R
        end do
    end do
    return
end subroutine Transpose


subroutine MatMul(A,B,C,arow,acol,brow,bcol)
!To see how this works, try writing the elements of the product matrix in sigma notation
dimension A(10,10),B(10,10),C(10,10) !CAUTION: This declaration is crucial
do k=1,arow !k runs from 1 to arow
    do l=1,bcol !and l from 1 to bcol because that's the output matrix dimension and we _have_ to populate each element of it
        sum=0
        do i=1,acol !this loop is the essence of the aforementioned sigma notation for the (k,l) element of the product matrix
            sum=sum+A(k,i)*B(i,l)
        end do
        C(k,l)=sum
    end do
end do
return
end subroutine MatMul