Differentiate between copy_oftype(A, T) and convert(AbstractArray{T}, A) for structured arrays

[daanhb]

This pull request redefines the copy_oftype function in LinearAlgebra and it implements convert(AbstractArray{T}, A) differently for a small number of structured arrays. The motivation for this pull request has been described in #34995. It also fixes #35071.

The request is split into two commits.

The first commit has the following two major changes:

• copy_oftype(A, T) is now implemented in terms of either copy(A) or copyto!(similar(A, T), A). It was already like that implicitly. The function used to invoke convert(AbstractArray{T}, A), which in turn called AbstractArray{T}(A), which in most cases would call copyto!(similar(A, T), A). The change makes the fallback explicit.
• AbstractMatrix{T}(A) is now implemented for all structured matrices in LinearAlgebra by invoking the AbstractVector{T} or AbstractMatrix{T} constructor on its underlying data. This was already the case for several arrays. After the first change, AbstractMatrix{T}(A) is no longer assumed to return something mutable.

It is difficult to illustrate with an example in Base, because the change is helpful for immutable arrays that support conversion to a different eltype. With the latest version of StaticArrays we have:

julia> using LinearAlgebra, StaticArrays

julia> convert(AbstractMatrix{Float64}, SA[1 2; 0 3])
2×2 SArray{Tuple{2,2},Float64,2,4} with indices SOneTo(2)×SOneTo(2):
 1.0  2.0
 0.0  3.0

julia> tri = UpperTriangular(SA[1 2; 0 3]);

julia> convert(AbstractMatrix{Float64}, tri)
2×2 UpperTriangular{Float64,SArray{Tuple{2,2},Float64,2,4}}:
 1.0  2.0
  ⋅   3.0

The point is that an immutable static matrix remains immutable after conversion to AbstractMatrix{T} (this is implemented in StaticArrays itself), and that the upper triangular matrix also remains immutable after conversion to AbstractMatrix{T} (this is implemented in the pull request). Previously, the conversion of the triangular matrix would surprisingly yield a mutable and less efficient MArray as the underlying data.

Design choices:

• I think it makes sense that copyto!(similar(A, T), A) should always work. This was not the case for the SymTridiagonal type (due to another issue which I'll come back to), so I've added a definition.
• AbstractMatrix{T}(A) mimicks what similar(A, T) does: it will preserve structure if similar preserves structure, and vice-versa. Thus, its application to a transposed vector is recursive, but its application to a transposed matrix is not.
• Should AbstractMatrix{T}(A) make a copy of A or not? It is clear that convert(AbstractArray{T}, A) is not expected to make a copy. The recursive implementation of AbstractMatrix{T}(A) in this pull request invokes the constructor on the underlying data, not the conversion, so the decision lies with the underlying data array and this pull request remains agnostic. I don't think the behaviour is uniform, but perhaps that is not such a big problem.

The second commit makes some more aggressive changes elsewhere in LinearAlgebra. It defines the new function copy_similar right after copy_oftype:

copy_similar(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = copyto!(similar(A, S, size(A)), A)

Typically, the three-argument call to similar discards the structure of A. The third argument is size(A) and it is subtly absent in the definition of copy_oftype. The result of copy_similar is used in many places to create a writable matrix without the structure of A which is suitable for an in-place algorithm. I've looked for this use-case and replaced all occurrences with a call to copy_similar. This hopefully improves the clarity of this subtle difference in the code, and also shortens some files (most changes are in triangular.jl). It is otherwise independent from the first commit.

On my machine the entire test suite passes with these changes, that is, it does not seem to affect any of the existing tests. I have yet to add tests for the specific things I've changed.

The code is actually shorter than all my explanations here and in #34995 :-)
Do these changes look acceptable?

[daanhb]

The issue with SymTridiagonal is the following:

julia> using LinearAlgebra

julia> A = SymTridiagonal(rand(4), rand(3))
4×4 SymTridiagonal{Float64,Array{Float64,1}}:
 0.75844    0.0577674   ⋅          ⋅
 0.0577674  0.0622232  0.909377    ⋅
  ⋅         0.909377   0.55185    0.0100493
  ⋅          ⋅         0.0100493  0.64809

julia> A[1,2] = 1
ERROR: ArgumentError: cannot set off-diagonal entry (1, 2)
Stacktrace:
 [1] setindex!(::SymTridiagonal{Float64,Array{Float64,1}}, ::Int64, ::Int64, ::Int64) at /home/daan/code/julia/julialang/julia/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/tridiag.jl:454
 [2] top-level scope at REPL[3]:1

julia> A[1,3] = 0
ERROR: ArgumentError: cannot set off-diagonal entry (1, 3)

I can not set the off-diagonal entry to another value, even though it seems valid to me. Also, the type does not allow setting other elements to zero. Both of these things make the default copyto! implementation fail, so I've added a definition of copyto! in the pull request. But perhaps setindex! itself should be changed?

[daanhb]

On second thought, not being able to set the off-diagonal entries of a SymTridiagonal matrix is correct, because it would have the non-generic side-effect of changing another entry as well (the symmetric one). The behaviour was clearly intended and is explicitly tested for.

[daanhb]

Adding tests has the benefit of going over the changes again, I fixed one error on my side.

The tests:

• check that AbstractArray{T}(A) (and hence convert(AbstractArray{T}, A)) preserves the underlying storage type, just like similar does, for the structured arrays in LinearAlgebra
• checks for Vector accidentally overwritten in one specialization of / #35071 (would this be better in a separate pull request?)
• check copyto! for SymTridiagonal, which was disabled before because it was not implemented

For testing purposes I needed an array for which similar(A,T) returns something different than convert(AbstractArray{T}, A). I did not find such a thing in Base. I then used a basic ImmutableArray type, which wraps another array and implements getindex but not setindex!.
However, I found no logical place to put it in the LinearAlgebra tests, because all the files can be run independently and ideally the definition of ImmutableArray is shared. For now, it is in src/LinearAlgebra.jl itself, to check that the tests work, but it obviously shouldn't stay there. Where else can it go, without having to copy it in all test files?

[dkarrasch]

@daanhb This is impressive work! We have a folder test/testhelpers, you could put your ImmutableArrays there, and then include the file, just like we include in several different tests Furlongs or Quaternions. I quickly scanned over the code, not yet in detail, and it looks good to me, though it would be good to have it in more final shape, possibly rebased and with all tests passing. 😄

[daanhb]

Thanks for the feedback @dkarrasch - I will update!

[daanhb]

Rebased and tests passing @dkarrasch. The ImmutableArrays are now in testhelpers, and I've included the file in the tests in the same way as was done for Furlongs.
Is there anything more that I can do here? My next steps would be to follow up on these changes in StaticArrays.jl.

[dkarrasch]

This got a bit lost. @daanhb Do you remember what was the last status? Ready for review and merge, or was there anything left for discussion/implementation? Maybe we should revive this by rebasing and fixing merge conflicts?

[daanhb]

Thanks, just leaving a comment here to say I still intend to rebase this pull request. At that point I'll reassess whether it is still useful and, if so, I'll bump here.

[daanhb]

I'm closing this in favour of #40831