Skip to content

Commit 546687f

Browse files
dlfivefiftydlfivefifty
committed
Support inv via QR
1 parent 089374f commit 546687f

File tree

6 files changed

+49
-24
lines changed

6 files changed

+49
-24
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.2.7"
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: 1 addition & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -451,8 +451,4 @@ function rmul!(block_array::BlockArray, α::Number)
451451
rmul!(block, α)
452452
end
453453
block_array
454-
end
455-
456-
# Temporary work around
457-
Base.reshape(block_array::BlockArray, axes::NTuple{N,AbstractUnitRange{Int}}) where N =
458-
reshape(PseudoBlockArray(block_array), axes)
454+
end

‎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: 25 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,11 @@ 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(block_array::PseudoBlockArray, axes::Tuple{Union{Integer,Base.OneTo}, Vararg{Union{Integer,Base.OneTo}}}) where N =
301+
_pseudo_reshape(block_array, axes)
293302

294303

295304
###########################

‎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 ≈ I(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 ≈ I(6)
193+
end
194+
end

0 commit comments

Comments
 (0)