BlockSkylineMatrix–Diagonal arithmetic, fix BlockDiagonal

[jagot]

I implemented a fix for #57 and made sure BlockDiagonal does not use BlockSizes anymore. I also exported BlockDiagonal.

The MulAdds still need to be fixed though:
https://github.com/JuliaMatrices/BlockBandedMatrices.jl/blob/657b48059ab3d5f720fdf7efb2cd9eecb1dfbdd3/src/linalg.jl#L113-L148

[jagot]

My fix for #57 basically mirrors LinearAlgebra's implementation for UniformScaling:

https://github.com/JuliaLang/julia/blob/18783434e9c0e48152b99d5e3717deb8711f492b/stdlib/LinearAlgebra/src/uniformscaling.jl#L163-L179

I don't know if that is a good idea, since LinearAlgebra.copy_oftype is not exported.

[dlfivefifty]

I use internal functionality in Base all the time, for some things there's no other way, for others its easier to just fix bugs introduced by this in future versions when they arise.

[jagot]

Fair enough. I should mention that I consciously did not try to support arithmetic with a BlockSkylineMatrix which does not have blocks on the diagonal. My rationale is that you want to be sure you get the same structure back (and figuring out which blocks to add is maybe ill-defined?, beyond growing the bandwidths to at least zero). What do you think?

[dlfivefifty]

What BlockArrays.jl does is to take the union of all blocks, that is, if you add a vector with block sizes [1,2,2] to one with block sizes [3,1,1] you get one with block sizes [1,2,1,1].

I feel like block skyline structure would be preserved under this, though perhaps I haven't thought it through enough

[jagot]

What I mean is this slightly odd case:

4×4-blocked 10×10 BlockSkylineMatrix{Float64,Array{Float64,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},BandedMatrices.BandedMatrix{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}}:
  ⋅   │   ⋅    ⋅   │   ⋅    ⋅    ⋅   │   ⋅    ⋅    ⋅    ⋅
 ─────┼────────────┼─────────────────┼────────────────────
 1.0  │   ⋅    ⋅   │   ⋅    ⋅    ⋅   │   ⋅    ⋅    ⋅    ⋅
 1.0  │   ⋅    ⋅   │   ⋅    ⋅    ⋅   │   ⋅    ⋅    ⋅    ⋅
 ─────┼────────────┼─────────────────┼────────────────────
 1.0  │  1.0  1.0  │   ⋅    ⋅    ⋅   │   ⋅    ⋅    ⋅    ⋅
 1.0  │  1.0  1.0  │   ⋅    ⋅    ⋅   │   ⋅    ⋅    ⋅    ⋅
 1.0  │  1.0  1.0  │   ⋅    ⋅    ⋅   │   ⋅    ⋅    ⋅    ⋅
 ─────┼────────────┼─────────────────┼────────────────────
  ⋅   │  1.0  1.0  │  1.0  1.0  1.0  │   ⋅    ⋅    ⋅    ⋅
  ⋅   │  1.0  1.0  │  1.0  1.0  1.0  │   ⋅    ⋅    ⋅    ⋅
  ⋅   │  1.0  1.0  │  1.0  1.0  1.0  │   ⋅    ⋅    ⋅    ⋅
  ⋅   │  1.0  1.0  │  1.0  1.0  1.0  │   ⋅    ⋅    ⋅    ⋅
10×10 Diagonal{Float64,Array{Float64,1}}:
 2.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅   2.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅   2.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅   2.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅   2.0   ⋅    ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅   2.0   ⋅    ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   2.0   ⋅    ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   2.0   ⋅    ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   2.0   ⋅
  ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   2.0

I think the easiest is to just grow the block-bandwidths to be at least 0.

[jagot]

I finally got around to fixing the general case as well, so now we can do:

julia> M
4×4-blocked 10×10 BlockSkylineMatrix{Int64,Array{Int64,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}},Array{Int64,1},Array{Int64,1},BandedMatrices.BandedMatrix{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}}:
 ⋅  │  ⋅  ⋅  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ───┼────────┼───────────┼────────────
 ⋅  │  ⋅  ⋅  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ⋅  │  ⋅  ⋅  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ───┼────────┼───────────┼────────────
 1  │  ⋅  ⋅  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 1  │  ⋅  ⋅  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 1  │  ⋅  ⋅  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ───┼────────┼───────────┼────────────
 ⋅  │  1  1  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ⋅  │  1  1  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ⋅  │  1  1  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ⋅  │  1  1  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅

julia> T
10×10 Tridiagonal{Int64,Array{Int64,1}}:
 2  3  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 1  2  3  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  1  2  3  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  1  2  3  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  1  2  3  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  1  2  3  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  1  2  3  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  2  3  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  2  3
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  2

julia> M+T
4×4-blocked 10×10 BlockSkylineMatrix{Int64,Array{Int64,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}},Array{Int64,1},Array{Int64,1},BandedMatrices.BandedMatrix{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}}:
 2  │  3  0  │  ⋅  ⋅  ⋅  │  ⋅  ⋅  ⋅  ⋅
 ───┼────────┼───────────┼────────────
 1  │  2  3  │  0  0  0  │  ⋅  ⋅  ⋅  ⋅
 0  │  1  2  │  3  0  0  │  ⋅  ⋅  ⋅  ⋅
 ───┼────────┼───────────┼────────────
 1  │  0  1  │  2  3  0  │  0  0  0  0
 1  │  0  0  │  1  2  3  │  0  0  0  0
 1  │  0  0  │  0  1  2  │  3  0  0  0
 ───┼────────┼───────────┼────────────
 ⋅  │  1  1  │  0  0  1  │  2  3  0  0
 ⋅  │  1  1  │  0  0  0  │  1  2  3  0
 ⋅  │  1  1  │  0  0  0  │  0  1  2  3
 ⋅  │  1  1  │  0  0  0  │  0  0  1  2

i.e. the block-bandwidths are grown as necessary. I've implemented + and - for combinations of BlockSkylineMatrix and Diagonal, Bidiagonal (upper and lower), Tridiagonal, and SymTridiagonal. It is trivial to add support different BandedMatrices as well, should we want to.

[codecov]

Codecov Report

Merging #58 into master will increase coverage by 2.30%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master      #58      +/-   ##
==========================================
+ Coverage   79.54%   81.84%   +2.30%     
==========================================
  Files          10       10              
  Lines         787      865      +78     
==========================================
+ Hits          626      708      +82     
+ Misses        161      157       -4     

Impacted Files	Coverage Δ
src/BlockSkylineMatrix.jl	94.93% <100.00%> (+3.32%)	⬆️
src/interfaceimpl.jl	94.28% <100.00%> (+7.18%)	⬆️
src/AbstractBlockBandedMatrix.jl	58.92% <0.00%> (+1.78%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update ff8d25a...f3d3d59. Read the comment docs.

[jagot]

I also added UniformScaling.

[jagot]

BTW, nightly is failing: https://travis-ci.org/github/JuliaMatrices/BlockBandedMatrices.jl/jobs/704067097#L326

[dlfivefifty]

Don’t worry about nightly. I’ll get 1.5 working once it’s rc2 as I’m waiting for a backport

[jagot]

Ok, then this PR is good to go from my point of view.