matrix multiplication not supported?

[benninkrs]

julia> A = OffsetArray([1 2; 3 4], 0:1, 0:1)
2×2 OffsetArray(::Array{Int64,2}, 0:1, 0:1) with eltype Int64 with indices 0:1×0:1:
 1  2
 3  4

julia> b = OffsetArray([5; 6], 0:1)
2-element OffsetArray(::Array{Int64,1}, 0:1) with eltype Int64 with indices 0:1:
 5
 6

julia> A*b
ERROR: ArgumentError: offset arrays are not supported but got an array with index other than 1
Stacktrace:
 [1] require_one_based_indexing at .\abstractarray.jl:89 [inlined]
 [2] generic_matvecmul!(::OffsetArray{Int64,1,Array{Int64,1}}, ::Char, ::OffsetArray{Int64,2,Array{Int64,2}}, ::OffsetArray{Int64,1,Array{Int64,1}}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\LinearAlgebra\src\matmul.jl:501
 [3] mul! at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\LinearAlgebra\src\matmul.jl:77 [inlined]
 [4] *(::OffsetArray{Int64,2,Array{Int64,2}}, ::OffsetArray{Int64,1,Array{Int64,1}}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.2\LinearAlgebra\src\matmul.jl:51
 [5] top-level scope at none:0

[timholy]

A pull request would be welcome. The easy version is to specialize it for two OffsetArray inputs with parents that have 1-based indexing, in which case you can strip the wrappers, call the ordinary routines, and then re-wrap.

The more ambitious version is to get this working generally in Julia itself. There are other packages that introduce arrays with unconventional axes, so it won't be a perfect solution unless we make it more general.

[benninkrs]

I'm not ready to tackle Base yet, but I'll look into making a PR for this here. A question about semantics arises, similar to #45. I see three options for A*B, from most conservative to most liberal:

1. Both A and B must be OffsetArrays, and axes(A,2) must "match" (be semantically equivalent to) axes(B,1).
2. Either A or B may be a non-OffsetArray, but the axes must match as in option 1. (This means the participating axis of the OffsetArray must start at 1.)
3. Either A or B may be a non-OffsetArray, and the contracted axes need only have the same length; the offset of the OffsetArray is irrelevant.

In all three cases, the result would be an OffsetArray with axes (axes(A,1), axes(B,2)).

Option 2 is probably just as safe as option 1; it effectively "promotes" the non-OffsetArray into an OffsetArray (Maybe this should be an actual promotion?)

The tempting motivation for option 3 is that sometimes I just want to apply a matrix A (say, some transform kernel) to a vector/matrix B and I don't know or don't care what the offset of B is. (Maybe B was provided by a function from some other package that doesn't use the same indexing convention the my code does.) This could perhaps be justified by the interpretation that an "ordinary" (non-offset) array is merely a data container whose indices simply provide a way of accessing different values but have no intrinsic meaning, whereas an OffsetArray expresses that index values are intrinsically meaningful. But even with this interpretation, I can see that option 3 has a problem: If A is a non-OffsetArray and B is an OffsetArray, then C=A*B is an OffsetArray (whose indices are purportedly meaningful), but the first axis of C is axes(A,1) whose indices were potentially not meaningful.

What do you think? After arguing with myself, I'm inclined to pursue option 2.

[timholy]

I like option 2 also. You can manually call no_view_offset and resolve the axes how you see fit when you want 3.