+ and - throw error with UpperTriangular matrices

[diegozea]

julia> m = [ 1 2 3
             4 5 6
             7 8 9 ]
3×3 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9

julia> UpperTriangular(m)
3×3 UpperTriangular{Int64,Array{Int64,2}}:
 1  2  3
 ⋅  5  6
 ⋅  ⋅  9

julia> UpperTriangular(m) + 2
ERROR: ArgumentError: cannot set index in the lower triangular part (2, 1) of an UpperTriangular matrix to a nonzero value (2)
 in setindex!(::UpperTriangular{Int64,Array{Int64,2}}, ::Int64, ::Int64, ::Int64) at ./linalg/triangular.jl:122
 in .+(::UpperTriangular{Int64,Array{Int64,2}}, ::Int64) at ./arraymath.jl:82
 in +(::UpperTriangular{Int64,Array{Int64,2}}, ::Int64) at ./arraymath.jl:94

julia> UpperTriangular(m) - 2
ERROR: ArgumentError: cannot set index in the lower triangular part (2, 1) of an UpperTriangular matrix to a nonzero value (-2)
 in setindex!(::UpperTriangular{Int64,Array{Int64,2}}, ::Int64, ::Int64, ::Int64) at ./linalg/triangular.jl:122
 in .-(::UpperTriangular{Int64,Array{Int64,2}}, ::Int64) at ./arraymath.jl:82
 in -(::UpperTriangular{Int64,Array{Int64,2}}, ::Int64) at ./arraymath.jl:96

julia> versioninfo()
Julia Version 0.5.0
Commit 3c9d753 (2016-09-19 18:14 UTC)
Platform Info:
  System: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5 CPU         750  @ 2.67GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Nehalem)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.7.1 (ORCJIT, nehalem)

[yuyichao]

So should these just return a full matrix?

[kshyatt]

@yuyichao yes, I think so. Probably something like +(a::Number, b::UpperTriangular) = a + Array(b) is the right way to go.

[tkelman]

not if the wrapped type is sparse.

[kshyatt]

Sorry - parent(A)? (This is why I said "something like" 😉 )

[andreasnoack]

@tkelman Is this the conclusion in one of the recent broadcast+sparse issues? (I have lost track of those discussions.) Right now, the behavior is mixed

julia> typeof(sprandn(10,10,0.1) .+ 1)
Array{Float64,2}

julia> typeof(broadcast(+, sprandn(10,10,0.1), 1))
SparseMatrixCSC{Float64,Int64}

but densifying seems more sensible for + and -.

[Sacha0]

#17623 fixed the inconsistency I believe

julia> typeof(sprandn(10,10,0.1) .+ 1)
SparseMatrixCSC{Float64,Int64}

Best!

[andreasnoack]

This works for the dot versions now that #17623 is merged. Anybody knows why these wouldn't just call broadcast?

[Sacha0]

Anybody knows why these wouldn't just call broadcast?

Perhaps for the special promote_op/promote_array_type behavior? Otherwise perhaps to avoid overhead (that may no longer exist)? Best!

[andreasnoack]

Maybe @stevengj would now if we can get rid of them.

[stevengj]

The promotion rules are slightly different for abstract types, so I didn't want to change them in that PR. It seemed best to do it in a separate PR.

[Sacha0]

Ref. related discussion https://github.com/JuliaLang/julia/issues/19595#issuecomment-268070182 and the forthcoming issue from @martinholters. Best!

[martinholters]

See #19669 re. the type behavior.
Some quick benchmarking shows that the comprehension is slightly slower than the explicit loop, which seems to be on par with broadcast, though. This is not a rigorous evaluation, just an indication that some more profound benchmarking might be in order.

[andreasnoack]

Fixed by #19746