I think the source of your issue isn’t matrix multiplication, it’s the fact that indexing with a scalar drops that dimension. This is consistent with the idea that in indexing the dimension of the result should match the dimensionality of the index. Scalars are 0-dimensional, so that dimension is dropped.
julia> A = rand(2,3,4);
julia> size(A[1, 1, 1])
()
julia> size(A[1:1, 1, 1])
(1,)
julia> size(A[1:1, 1:1, 1])
(1, 1)
julia> size(A[1:1, [1, 2], 1])
(1, 2)
julia> size(A[1:1, [1 2; 2 3], 1])
(1, 2, 2)
I think to get the behavior you want you’d want to overload getindex(::Array), which is the function that’s called for indexing, e.g. the method called for A[1, 2] is getindex(::Array, ::Int64, ::Int64). You can see this with @which:
julia> @which A[1,1]
getindex(A::Array, i1::Int64, i2::Int64, I::Int64...) in Base at array.jl:745
If you overloaded getindex so that indexing with a scalar kept that dimension, then M[:, 1] would give you a Nx1 Matrix, and M[1,:] would give you a 1xN Matrix.
You can poke around at that code to see how the dispatch happens if you want to tweak things. What other folks are warning about is that this will likely break other people’s code if you use any other packages (or if you publish a package that someone else uses). That said, if you’re interested in this, there’s nothing wrong with playing around with it, and likely learning more about Julia internals along the way. Just know that you’re modifying the global method table, so you may end up getting into a state where you want to restart your interpreter to get things back to the default state.
Happy hacking. 