hvcat indexing issue with misaligned blocks

[JeffFessler]

While doing more thorough testing of LinearMapsAA I uncovered this issue with the new blockmap features (that are otherwise working great):

A = LinearMap(ones(3,2))
B = [I A; A I]
Matrix(B) # works fine
Matrix(B') # fails with error DimensionMismatch("mul!")

Now, one could argue that the way B is constructed in this MWE, with blocks that do not quite line up, is not the way we write block matrices in math, and I agree. But that type of hvcat works fine with ordinary arrays, because the two block rows have the same number of columns. I'm not sure it is worth trying to support multiplication with such constructs. Perhaps LinearMaps.hvcat could just throw an error saying that such constructions are unsupported, e.g., by checking that all blocks in a given column have the same number of columns.

[dkarrasch]

Hm, I’m confused. Shouldn’t calling B' already throw with an informative error message, like inconsistent number of rows or something like this? At least, this is what I would want to happen.

[dkarrasch]

Oh no, I forgot. At adjoining, we don’t construct a new BlockMap, but simply wrap it as an AdjointMap. The reason is that „transposing" the linearly ordered tuple of maps in a type-stable fashion was impossible for me. It was much easier and possible to write At_mul_B and Ac_mul_B multiplication code for virtually transposed/adjoint block maps. So, indeed, the issue with your example becomes apparent only after calling some sort of multiplication code like Matrix. And this makes it very hard to generate an error message: generally, we don’t insist on BlockMaps being „block-square“ in the first place (consistently with the hvcat code for UniformScalings and AbstractMatrixes), which in turn makes a column-wise check unclear. Also, the construction of B is supported, consistently with the matrix case in Julia Base : A = ones(3,2); B = [I A; A I].

[JeffFessler]

Ah, I see now that perhaps what would be needed is some modification of At_mul_B and Ac_mul_B to either handle such cases or throw errors there. I'm not sure if it is worth doing unless someone has a use case. I'm removing those cases from my test suite... In case someone comes looking later, the specific error message is ERROR: DimensionMismatch("mul!") from this line: https://github.com/Jutho/LinearMaps.jl/blob/master/src/blockmap.jl#L283 Thanks.

[dkarrasch]

I slept over it, and think we should consider including some size logic in the construction of transposes and adjoints. In your original example, the construction of B is fine, it's the construction of B' that yields a LinearMap that is no good. So the error message is then slightly misleading, because it makes it look like there is an issue with the multiplication, but the actual issue is with the adjoint/transpose itself, because it is not of block type (in that example, one would have to split and merge blocks).

I guess the question is whether we want to make the construction of AdjointMap/TransposeMap more expensive by introducing size-checking logic, or whether we leave things as they are at the expense of a not fully explicit error message.

[JeffFessler]

Actually, the existing indexing is fine. There is a small bug in the transpose multiplies. Will submit PR soon. See #59