getindex on WoodburyPDMat should return a WoodburyPDMat

[glennmoy]

getindex returns a generic matrix, which loses the useful decomposition contained in the WoodburyPDMat.

Maybe it would be a good idea to return another WoodburyPDMat? Although this might only make sense in certain instances e.g. extracting a square sub-matrix and not a slice.

MWE

julia> A = 0.5 * ones(4, 3);

julia> S = 0.5 * ones(4);

julia> D = ones(3);

julia> w = WoodburyPDMat(A, D, S);

julia> w[[1, 2], [1, 2]]
2×2 Matrix{Float64}:
 1.25  0.75
 0.75  1.25

[glennmoy]

although, this doesn't seem to be the case with other PDMats so maybe it's not a fair requirement

julia> m = PDMats.PDMat(Symmetric([1 0 0; 0 1 0; 0 0 1]))
3×3 PDMat{Float64, Matrix{Float64}}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> m[[1, 2], [1, 3]]
2×2 Matrix{Float64}:
 1.0  0.0
 0.0  0.0

[glennmoy]

An alternative proposal to the problem this would solve, might be to have a WoodburyPDMat constructor that takes a square matrix and computes the A, S, and D matrices, sort of how cholesky does. The you could just wrap your result from indexing to get back a Woodbury.