Skip to content

Commit 2424ea1

Browse files
dlfivefiftydlfivefifty
authored
Support inv via QR (#113)
* Revert "Revert "Support inv via QR"" This reverts commit b1add34. * Fix reshape ambiguities * Update test_blocklinalg.jl * tests for reshape overrides
1 parent b1add34 commit 2424ea1

File tree

7 files changed

+81
-34
lines changed

7 files changed

+81
-34
lines changed

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
77
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
88

99
[compat]
10-
ArrayLayouts = "0.2.4"
10+
ArrayLayouts = "0.3"
1111
julia = "1.1"
1212

1313
[extras]

‎src/BlockArrays.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@ import LinearAlgebra: lmul!, rmul!, AbstractTriangular, HermOrSym, AdjOrTrans,
3535
StructuredMatrixStyle
3636
import ArrayLayouts: _fill_lmul!, MatMulVecAdd, MatMulMatAdd, MatLmulVec, MatLdivVec,
3737
materialize!, MemoryLayout, sublayout, transposelayout, conjlayout,
38-
triangularlayout, triangulardata
38+
triangularlayout, triangulardata, _inv
3939

4040
include("blockindices.jl")
4141
include("blockaxis.jl")

‎src/blockarray.jl

Lines changed: 21 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -74,11 +74,11 @@ end
7474
_BlockArray(blocks::R, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}} =
7575
_BlockArray(blocks, map(blockedrange, block_sizes))
7676

77-
# support non-concrete eltypes in blocks
77+
# support non-concrete eltypes in blocks
7878
_BlockArray(blocks::R, block_axes::BS) where {T, N, R<:AbstractArray{<:AbstractArray{V,N} where V,N}, BS<:NTuple{N,AbstractUnitRange{Int}}} =
7979
_BlockArray(convert(AbstractArray{AbstractArray{mapreduce(eltype,promote_type,blocks),N},N}, blocks), block_axes)
8080
_BlockArray(blocks::R, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R<:AbstractArray{<:AbstractArray{V,N} where V,N}} =
81-
_BlockArray(convert(AbstractArray{AbstractArray{mapreduce(eltype,promote_type,blocks),N},N}, blocks), block_sizes...)
81+
_BlockArray(convert(AbstractArray{AbstractArray{mapreduce(eltype,promote_type,blocks),N},N}, blocks), block_sizes...)
8282

8383
const BlockMatrix{T, R <: AbstractMatrix{<:AbstractMatrix{T}}} = BlockArray{T, 2, R}
8484
const BlockVector{T, R <: AbstractVector{<:AbstractVector{T}}} = BlockArray{T, 1, R}
@@ -165,7 +165,7 @@ initialized_blocks_BlockArray(::Type{R}, block_sizes::Vararg{AbstractVector{Int}
165165

166166

167167
@inline BlockArray{T,N,R,BS}(::UndefInitializer, sizes::NTuple{N,Int}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}, BS<:NTuple{N,AbstractUnitRange{Int}}} =
168-
BlockArray{T,N,R,BS}(undef, convert(BS, map(Base.OneTo, sizes)))
168+
BlockArray{T,N,R,BS}(undef, convert(BS, map(Base.OneTo, sizes)))
169169

170170
function BlockArray{T}(arr::AbstractArray{V, N}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T,V,N}
171171
for i in 1:N
@@ -176,7 +176,7 @@ function BlockArray{T}(arr::AbstractArray{V, N}, block_sizes::Vararg{AbstractVec
176176
BlockArray{T}(arr, map(blockedrange,block_sizes))
177177
end
178178

179-
BlockArray(arr::AbstractArray{T, N}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T,N} =
179+
BlockArray(arr::AbstractArray{T, N}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T,N} =
180180
BlockArray{T}(arr, block_sizes...)
181181

182182
@generated function BlockArray{T}(arr::AbstractArray{T, N}, baxes::NTuple{N,AbstractUnitRange{Int}}) where {T,N}
@@ -347,22 +347,22 @@ end
347347

348348

349349
@inline Base.similar(block_array::AbstractArray, ::Type{T}, axes::Tuple{BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
350-
BlockArray{T}(undef, axes)
350+
BlockArray{T}(undef, axes)
351351
@inline Base.similar(block_array::AbstractArray, ::Type{T}, axes::Tuple{BlockedUnitRange,BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
352-
BlockArray{T}(undef, axes)
352+
BlockArray{T}(undef, axes)
353353
@inline Base.similar(block_array::AbstractArray, ::Type{T}, axes::Tuple{AbstractUnitRange{Int},BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
354-
BlockArray{T}(undef, axes)
355-
354+
BlockArray{T}(undef, axes)
355+
356356
@inline Base.similar(block_array::Type{<:AbstractArray{T}}, axes::Tuple{BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
357-
BlockArray{T}(undef, axes)
357+
BlockArray{T}(undef, axes)
358358
@inline Base.similar(block_array::Type{<:AbstractArray{T}}, axes::Tuple{BlockedUnitRange,BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
359-
BlockArray{T}(undef, axes)
359+
BlockArray{T}(undef, axes)
360360
@inline Base.similar(block_array::Type{<:AbstractArray{T}}, axes::Tuple{AbstractUnitRange{Int},BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
361-
BlockArray{T}(undef, axes)
362-
361+
BlockArray{T}(undef, axes)
362+
363363
const OffsetAxis = Union{Integer, UnitRange, Base.OneTo, Base.IdentityUnitRange, Colon}
364364

365-
# avoid ambiguities
365+
# avoid ambiguities
366366
@inline Base.similar(block_array::BlockArray, ::Type{T}, dims::NTuple{N,Int}) where {T,N} =
367367
Array{T}(undef, dims)
368368
@inline Base.similar(block_array::BlockArray, ::Type{T}, axes::Tuple{OffsetAxis,Vararg{OffsetAxis}}) where T =
@@ -454,5 +454,11 @@ function rmul!(block_array::BlockArray, α::Number)
454454
end
455455

456456
# Temporary work around
457-
Base.reshape(block_array::BlockArray, axes::NTuple{N,AbstractUnitRange{Int}}) where N =
458-
reshape(PseudoBlockArray(block_array), axes)
457+
Base.reshape(block_array::BlockArray, axes::NTuple{N,AbstractUnitRange{Int}}) where N =
458+
reshape(PseudoBlockArray(block_array), axes)
459+
Base.reshape(block_array::BlockArray, dims::Tuple{Int,Vararg{Int}}) =
460+
reshape(PseudoBlockArray(block_array), dims)
461+
Base.reshape(block_array::BlockArray, axes::Tuple{Union{Integer,Base.OneTo}, Vararg{Union{Integer,Base.OneTo}}}) =
462+
reshape(PseudoBlockArray(block_array), axes)
463+
Base.reshape(block_array::BlockArray, dims::Tuple{Vararg{Union{Int,Colon}}}) =
464+
reshape(PseudoBlockArray(block_array), dims)

‎src/blocklinalg.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -255,3 +255,6 @@ for UNIT in ('U', 'N')
255255
end
256256
end
257257
end
258+
259+
# For now, use PseudoBlockArray
260+
_inv(::AbstractBlockLayout, axes, A) = BlockArray(inv(PseudoBlockArray(A)))

‎src/pseudo_blockarray.jl

Lines changed: 27 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -30,9 +30,9 @@ julia> A = PseudoBlockArray(rand(2,3), [1,1], [2,1])
3030
3131
julia> A = PseudoBlockArray(sprand(6, 0.5), [3,2,1])
3232
3-blocked 6-element PseudoBlockArray{Float64,1,SparseVector{Float64,Int64},Tuple{BlockedUnitRange{Array{Int64,1}}}}:
33-
0.0
33+
0.0
3434
0.5865981007905481
35-
0.0
35+
0.0
3636
───────────────────
3737
0.05016684053503706
3838
0.0
@@ -56,22 +56,22 @@ const PseudoBlockVecOrMat{T} = Union{PseudoBlockMatrix{T}, PseudoBlockVector{T}}
5656
PseudoBlockArray{T, N, R,BS}(blocks, baxes)
5757

5858
@inline PseudoBlockArray{T}(blocks::R, baxes::BS) where {T,N,R<:AbstractArray{T,N},BS<:NTuple{N,AbstractUnitRange{Int}}} =
59-
PseudoBlockArray{T, N, R,BS}(blocks, baxes)
60-
59+
PseudoBlockArray{T, N, R,BS}(blocks, baxes)
60+
6161
@inline PseudoBlockArray{T}(blocks::AbstractArray{<:Any,N}, baxes::NTuple{N,AbstractUnitRange{Int}}) where {T,N} =
62-
PseudoBlockArray{T}(convert(AbstractArray{T,N}, blocks), baxes)
62+
PseudoBlockArray{T}(convert(AbstractArray{T,N}, blocks), baxes)
6363

6464
@inline PseudoBlockArray(blocks::PseudoBlockArray, baxes::BS) where {N,BS<:NTuple{N,AbstractUnitRange{Int}}} =
6565
PseudoBlockArray(blocks.blocks, baxes)
6666

6767
@inline PseudoBlockArray{T}(blocks::PseudoBlockArray, baxes::BS) where {T,N,BS<:NTuple{N,AbstractUnitRange{Int}}} =
68-
PseudoBlockArray{T}(blocks.blocks, baxes)
68+
PseudoBlockArray{T}(blocks.blocks, baxes)
6969

7070
PseudoBlockArray(blocks::AbstractArray{T, N}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N} =
7171
PseudoBlockArray(blocks, map(blockedrange,block_sizes))
7272

7373
PseudoBlockArray{T}(blocks::AbstractArray{<:Any, N}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N} =
74-
PseudoBlockArray{T}(blocks, map(blockedrange,block_sizes))
74+
PseudoBlockArray{T}(blocks, map(blockedrange,block_sizes))
7575

7676
@inline PseudoBlockArray{T,N,R,BS}(::UndefInitializer, baxes::NTuple{N,AbstractUnitRange{Int}}) where {T,N,R,BS<:NTuple{N,AbstractUnitRange{Int}}} =
7777
PseudoBlockArray{T,N,R,BS}(R(undef, length.(baxes)), convert(BS, baxes))
@@ -117,7 +117,7 @@ convert(::Type{PseudoBlockArray{T,N}}, A::PseudoBlockArray{T,N}) where {T,N} = A
117117
convert(::Type{PseudoBlockArray{T}}, A::PseudoBlockArray{T}) where {T} = A
118118
convert(::Type{PseudoBlockArray}, A::PseudoBlockArray) = A
119119

120-
convert(::Type{PseudoBlockArray{T,N,R,BS}}, A::PseudoBlockArray) where {T,N,R,BS} =
120+
convert(::Type{PseudoBlockArray{T,N,R,BS}}, A::PseudoBlockArray) where {T,N,R,BS} =
121121
PseudoBlockArray{T,N,R,BS}(convert(R, A.blocks), convert(BS, A.axes))
122122

123123

@@ -144,18 +144,25 @@ function Base.similar(block_array::PseudoBlockArray{T,N}, ::Type{T2}) where {T,N
144144
end
145145

146146
@inline Base.similar(block_array::Type{<:Array{T}}, axes::Tuple{BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
147-
PseudoBlockArray{T}(undef, axes)
147+
PseudoBlockArray{T}(undef, axes)
148148
@inline Base.similar(block_array::Type{<:Array{T}}, axes::Tuple{BlockedUnitRange,BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
149-
PseudoBlockArray{T}(undef, axes)
149+
PseudoBlockArray{T}(undef, axes)
150150
@inline Base.similar(block_array::Type{<:Array{T}}, axes::Tuple{AbstractUnitRange{Int},BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
151-
PseudoBlockArray{T}(undef, axes)
151+
PseudoBlockArray{T}(undef, axes)
152152

153153
@inline Base.similar(block_array::Array, ::Type{T}, axes::Tuple{BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
154-
PseudoBlockArray{T}(undef, axes)
154+
PseudoBlockArray{T}(undef, axes)
155155
@inline Base.similar(block_array::Array, ::Type{T}, axes::Tuple{BlockedUnitRange,BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
156-
PseudoBlockArray{T}(undef, axes)
156+
PseudoBlockArray{T}(undef, axes)
157157
@inline Base.similar(block_array::Array, ::Type{T}, axes::Tuple{AbstractUnitRange{Int},BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
158-
PseudoBlockArray{T}(undef, axes)
158+
PseudoBlockArray{T}(undef, axes)
159+
160+
@inline Base.similar(block_array::PseudoBlockArray, ::Type{T}, axes::Tuple{BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
161+
PseudoBlockArray{T}(undef, axes)
162+
@inline Base.similar(block_array::PseudoBlockArray, ::Type{T}, axes::Tuple{BlockedUnitRange,BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
163+
PseudoBlockArray{T}(undef, axes)
164+
@inline Base.similar(block_array::PseudoBlockArray, ::Type{T}, axes::Tuple{AbstractUnitRange{Int},BlockedUnitRange,Vararg{AbstractUnitRange{Int}}}) where T =
165+
PseudoBlockArray{T}(undef, axes)
159166

160167
@inline function Base.getindex(block_arr::PseudoBlockArray{T, N}, i::Vararg{Integer, N}) where {T,N}
161168
@boundscheck checkbounds(block_arr, i...)
@@ -287,9 +294,13 @@ function rmul!(block_array::PseudoBlockArray, α::Number)
287294
block_array
288295
end
289296

297+
_pseudo_reshape(block_array, axes) = PseudoBlockArray(reshape(block_array.blocks,map(length,axes)),axes)
290298
Base.reshape(block_array::PseudoBlockArray, axes::NTuple{N,AbstractUnitRange{Int}}) where N =
291-
PseudoBlockArray(reshape(block_array.blocks,map(length,axes)),axes)
292-
299+
_pseudo_reshape(block_array, axes)
300+
Base.reshape(parent::PseudoBlockArray, shp::Tuple{Union{Integer,Base.OneTo}, Vararg{Union{Integer,Base.OneTo}}}) where N =
301+
reshape(parent, Base.to_shape(shp))
302+
Base.reshape(parent::PseudoBlockArray, dims::Tuple{Int,Vararg{Int}}) =
303+
Base._reshape(parent, dims)
293304

294305

295306
###########################

‎test/test_blockarrays.jl

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -105,6 +105,7 @@ end
105105
@test similar(Array{Float64}, axes(ret)) isa PseudoBlockArray
106106
@test similar(Vector{Float64}, axes(ret)) isa PseudoBlockArray
107107
@test similar(randn(5,5), Float64, axes(ret)) isa PseudoBlockArray
108+
@test similar(ret, Float64, (Base.IdentityUnitRange(1:3),)) isa BlockArray
108109

109110
ret = BlockArray{Float64}(undef, 1:3, 1:3)
110111
@test similar(typeof(ret), axes(ret)) isa BlockArray
@@ -371,6 +372,11 @@ end
371372
C = convert(PseudoBlockArray{Float32, 2}, A)
372373
@test C ≈ A ≈ PseudoBlockArray(A)
373374
@test eltype(C) == Float32
375+
376+
@test convert(BlockArray, A) === A
377+
@test convert(BlockArray{Float64}, A) === A
378+
@test convert(BlockMatrix{Float64}, A) === A
379+
@test convert(BlockMatrix{Float64,Matrix{Matrix{Float64}}}, A) === A
374380
end
375381

376382
@testset "string" begin
@@ -495,11 +501,15 @@ end
495501
@test reshape(A, Val(2)) isa PseudoBlockArray{Int64,2,Array{Int64,2},Tuple{BlockedUnitRange{Array{Int64,1}},Base.OneTo{Int64}}}
496502
@test reshape(A, Val(2)) == PseudoBlockArray(reshape(1:6,6,1), (blockedrange(1:3), Base.OneTo(1)))
497503
@test reshape(A, (blockedrange(Fill(2,3)),))[Block(1)] == 1:2
504+
@test reshape(A, 2, 3) == reshape(A, Base.OneTo(2), 3) == reshape(Vector(A), 2, 3)
505+
506+
@test_throws DimensionMismatch reshape(A,3)
498507

499508
A = PseudoBlockArray(1:6, 1:3)
500509
@test reshape(A, Val(2)) isa typeof(PseudoBlockArray(reshape(1:6,6,1), (blockedrange(1:3), Base.OneTo(1))))
501510
@test reshape(A, Val(2)) == PseudoBlockArray(reshape(1:6,6,1), (blockedrange(1:3), Base.OneTo(1)))
502511
@test reshape(A, (blockedrange(Fill(2,3)),))[Block(1)] == 1:2
512+
@test reshape(A, 2, 3) == reshape(A, Base.OneTo(2), 3) == reshape(Vector(A), 2, 3)
503513
end
504514

505515
@testset "*" begin

‎test/test_blocklinalg.jl

Lines changed: 18 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -172,6 +172,23 @@ import ArrayLayouts: DenseRowMajor
172172
@test UnitLowerTriangular(B) \ b ≈ UnitLowerTriangular(Matrix(B)) \ b
173173
end
174174
end
175-
end
176175

176+
@testset "inv" begin
177+
A = PseudoBlockArray{Float64}(randn(6,6), fill(2,3), 1:3)
178+
F = factorize(A)
179+
# Defaults to QR for generality
180+
@test F isa LinearAlgebra.QRCompactWY
181+
B = randn(6,6)
182+
@test ldiv!(F, copy(B)) ≈ Matrix(A) \ B
183+
B̃ = PseudoBlockArray(copy(B),1:3,fill(2,3))
184+
@test ldiv!(F, B̃) ≈ A\B ≈ Matrix(A) \ B
185+
186+
@test inv(A) isa PseudoBlockArray
187+
@test inv(A) ≈ inv(Matrix(A))
188+
@test inv(A)*A ≈ Matrix(I,6,6)
177189

190+
A = BlockArray{Float64}(randn(6,6), fill(2,3), 1:3)
191+
@test inv(A) isa BlockArray
192+
@test inv(A)*A ≈ Matrix(I,6,6)
193+
end
194+
end

0 commit comments

Comments
 (0)