Skip to content

Commit ec5aac0

Browse files
dlfivefiftydlfivefifty
authored
fix inf-broadcasting with blocked fills (#52)
* fix inf-broadcasting with blocked fills * Update Project.toml * Support BlockInterlace * Update Project.toml * Update for RangeCumsum * work on TriangleRecurrence * move map for mortar to LazyArrays.jl * save * FillArrays v0.11 * Update runtests.jl * Update runtests.jl * Update blockbanded.jl * tests pass * Update runtests.jl * Update Project.toml * Update Project.toml * Update Project.toml * use _power_by_squaring * resave
1 parent 7b2e59f commit ec5aac0

File tree

6 files changed

+229
-86
lines changed

6 files changed

+229
-86
lines changed

‎Project.toml

Lines changed: 13 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "InfiniteLinearAlgebra"
22
uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
3-
version = "0.4.4"
3+
version = "0.4.5"
44

55
[deps]
66
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -16,22 +16,23 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
1616
SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
1717

1818
[compat]
19-
ArrayLayouts = "0.4.10"
20-
BandedMatrices = "0.15.23"
21-
BlockArrays = "0.12.12"
22-
BlockBandedMatrices = "0.9"
23-
FillArrays = "0.10"
24-
InfiniteArrays = "0.8.3"
25-
LazyArrays = "0.19"
26-
LazyBandedMatrices = "0.3.2"
27-
MatrixFactorizations = "0.6, 0.7, 0.8"
28-
SemiseparableMatrices = "0.2"
19+
ArrayLayouts = "0.5.1"
20+
BandedMatrices = "0.16"
21+
BlockArrays = "0.14"
22+
BlockBandedMatrices = "0.10"
23+
FillArrays = "0.11"
24+
InfiniteArrays = "0.9"
25+
LazyArrays = "0.20"
26+
LazyBandedMatrices = "0.4"
27+
MatrixFactorizations = "0.7.1, 0.8"
28+
SemiseparableMatrices = "0.2.2"
2929
julia = "1.5"
3030

3131
[extras]
3232
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
3333
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
34+
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
3435
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
3536

3637
[targets]
37-
test = ["Test", "Random", "SpecialFunctions"]
38+
test = ["Test", "Random", "SpecialFunctions", "StaticArrays"]

‎src/InfiniteLinearAlgebra.jl

Lines changed: 5 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -9,11 +9,11 @@ import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted
99

1010
import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, TridiagonalLayout,
1111
sublayout, _qr, __qr, MatLmulVec, MatLmulMat, AbstractQLayout, materialize!, diagonaldata, subdiagonaldata, supdiagonaldata,
12-
_bidiag_forwardsub!, mulreduce
12+
_bidiag_forwardsub!, mulreduce, RangeCumsum
1313
import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
1414
_default_banded_broadcast, banded_similar
1515
import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row
16-
import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange, AbstractInfUnitRange, InfAxes
16+
import InfiniteArrays: OneToInf, InfUnitRange, Infinity, InfStepRange, AbstractInfUnitRange, InfAxes, InfRanges
1717
import LinearAlgebra: lmul!, ldiv!, matprod, qr, AbstractTriangular, AbstractQ, adjoint, transpose
1818
import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector, DenseColumnMajor, FillLayout, ApplyMatrix, check_mul_axes, ApplyStyle, LazyArrayApplyStyle, LazyArrayStyle,
1919
resizedata!, MemoryLayout,
@@ -23,11 +23,11 @@ import LazyArrays: applybroadcaststyle, CachedArray, CachedMatrix, CachedVector,
2323
import MatrixFactorizations: ul, ul!, _ul, ql, ql!, _ql, QLPackedQ, getL, getR, getU, reflector!, reflectorApply!, QL, QR, QRPackedQ,
2424
QRPackedQLayout, AdjQRPackedQLayout, QLPackedQLayout, AdjQLPackedQLayout, LayoutQ
2525

26-
import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange, blockcolsupport, BlockLayout, AbstractBlockLayout
26+
import BlockArrays: AbstractBlockVecOrMat, sizes_from_blocks, _length, BlockedUnitRange, blockcolsupport, BlockLayout, AbstractBlockLayout, BlockSlice
2727

2828
import BandedMatrices: BandedMatrix, bandwidths, AbstractBandedLayout, _banded_qr!, _banded_qr, _BandedMatrix
2929

30-
import LazyBandedMatrices: ApplyBandedLayout, BroadcastBandedLayout, _krontrav_axes
30+
import LazyBandedMatrices: ApplyBandedLayout, BroadcastBandedLayout, _krontrav_axes, _block_interlace_axes
3131

3232
import BlockBandedMatrices: _BlockSkylineMatrix, _BandedMatrix, _BlockSkylineMatrix, blockstart, blockstride,
3333
BlockSkylineSizes, BlockSkylineMatrix, BlockBandedMatrix, _BlockBandedMatrix, BlockTridiagonal,
@@ -38,7 +38,7 @@ import SemiseparableMatrices: AbstractAlmostBandedLayout, _almostbanded_qr!
3838

3939
# BroadcastStyle(::Type{<:BandedMatrix{<:Any,<:Any,<:OneToInf}}) = LazyArrayStyle{2}()
4040

41-
function ^(A::BandedMatrix{T,<:Any,<:OneToInf}, p::Integer) where T
41+
function ArrayLayouts._power_by_squaring(_, ::NTuple{2,Infinity}, A::AbstractMatrix{T}, p::Integer) where T
4242
if p < 0
4343
inv(A)^(-p)
4444
elseif p == 0
@@ -61,38 +61,5 @@ include("infql.jl")
6161
include("infqr.jl")
6262
include("inful.jl")
6363

64-
#######
65-
# block broadcasted
66-
######
67-
68-
const CumsumOneToInf2 = Cumsum{<:Any,1,<:OneToInf}
69-
BlockArrays.sortedunion(a::CumsumOneToInf2, ::CumsumOneToInf2) = a
70-
71-
function BlockArrays.sortedunion(a::Vcat{Int,1,<:Tuple{<:AbstractVector{Int},InfStepRange{Int,Int}}},
72-
b::Vcat{Int,1,<:Tuple{<:AbstractVector{Int},InfStepRange{Int,Int}}})
73-
@assert a == b
74-
a
75-
end
76-
77-
78-
map(::typeof(length), A::BroadcastArray{OneTo{Int},1,Type{OneTo}}) = A.args[1]
79-
map(::typeof(length), A::BroadcastArray{<:Fill,1,Type{Fill}}) = A.args[2]
80-
map(::typeof(length), A::BroadcastArray{<:Zeros,1,Type{Zeros}}) = A.args[1]
81-
map(::typeof(length), A::BroadcastArray{<:Vcat,1,Type{Vcat}}) = broadcast(+,map.(length,A.args)...)
82-
broadcasted(::LazyArrayStyle{1}, ::typeof(length), A::BroadcastArray{OneTo{Int},1,Type{OneTo}}) =
83-
A.args[1]
84-
broadcasted(::LazyArrayStyle{1}, ::typeof(length), A::BroadcastArray{<:Fill,1,Type{Fill}}) =
85-
A.args[2]
86-
87-
BlockArrays._length(::BlockedUnitRange, ::OneToInf) = ∞
88-
BlockArrays._last(::BlockedUnitRange, ::OneToInf) = ∞
89-
90-
###
91-
# KronTrav
92-
###
93-
94-
_krontrav_axes(A::NTuple{N,OneToInf{Int}}, B::NTuple{N,OneToInf{Int}}) where N =
95-
@. blockedrange(OneTo(length(A)))
96-
9764

9865
end # module

‎src/banded/infbanded.jl

Lines changed: 13 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -373,14 +373,19 @@ _BandedMatrix(::PertToeplitzLayout, A::AbstractMatrix) =
373373
# UniformScaling
374374
##
375375

376-
for op in (:+, :-), Typ in (:(BandedMatrix{<:Any,<:Any,OneToInf{Int}}),
377-
:(Adjoint{<:Any,<:BandedMatrix{<:Any,<:Any,OneToInf{Int}}}),
378-
:(Transpose{<:Any,<:BandedMatrix{<:Any,<:Any,OneToInf{Int}}}))
379-
@eval begin
380-
$op(A::$Typ, λ::UniformScaling) = $op(A, Diagonal(Fill(λ.λ,∞)))
381-
$op(λ::UniformScaling, A::$Typ) = $op(Diagonal(Fill(λ.λ,∞)), A)
382-
end
383-
end
376+
# for op in (:+, :-), Typ in (:(BandedMatrix{<:Any,<:Any,OneToInf{Int}}),
377+
# :(Adjoint{<:Any,<:BandedMatrix{<:Any,<:Any,OneToInf{Int}}}),
378+
# :(Transpose{<:Any,<:BandedMatrix{<:Any,<:Any,OneToInf{Int}}}))
379+
# @eval begin
380+
# $op(A::$Typ, λ::UniformScaling) = $op(A, Diagonal(Fill(λ.λ,∞)))
381+
# $op(λ::UniformScaling, A::$Typ) = $op(Diagonal(Fill(λ.λ,∞)), A)
382+
# end
383+
# end
384+
385+
386+
387+
ArrayLayouts._apply(_, ::NTuple{2,Infinity}, op, Λ::UniformScaling, A::AbstractMatrix) = op(Diagonal(Fill(Λ.λ,∞)), A)
388+
ArrayLayouts._apply(_, ::NTuple{2,Infinity}, op, A::AbstractMatrix, Λ::UniformScaling) = op(A, Diagonal(Fill(Λ.λ,∞)))
384389

385390
_default_banded_broadcast(bc::Broadcasted, ::Tuple{<:OneToInf,<:Any}) = copy(Broadcasted{LazyArrayStyle{2}}(bc.f, bc.args))
386391

‎src/blockbanded/blockbanded.jl

Lines changed: 59 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,67 @@
1+
const OneToInfCumsum = InfiniteArrays.RangeCumsum{Int,OneToInf{Int}}
2+
const OneToCumsum = InfiniteArrays.RangeCumsum{Int,OneTo{Int}}
3+
4+
BlockArrays.sortedunion(a::OneToInfCumsum, ::OneToInfCumsum) = a
5+
BlockArrays.sortedunion(a::OneToCumsum, ::OneToCumsum) = a
6+
7+
function BlockArrays.sortedunion(a::Vcat{Int,1,<:Tuple{<:AbstractVector{Int},InfStepRange{Int,Int}}},
8+
b::Vcat{Int,1,<:Tuple{<:AbstractVector{Int},InfStepRange{Int,Int}}})
9+
@assert a == b
10+
a
11+
end
12+
113
sizes_from_blocks(A::AbstractVector, ::Tuple{OneToInf{Int}}) = (map(length,A),)
214

15+
const OneToInfBlocks = BlockedUnitRange{OneToInfCumsum}
16+
const OneToBlocks = BlockedUnitRange{OneToCumsum}
17+
18+
axes(a::OneToInfBlocks) = (a,)
19+
axes(a::OneToBlocks) = (a,)
20+
21+
22+
function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{Ones{T,1,Tuple{OneToInfBlocks}},AbstractArray{V,N}}}) where {N,T,V}
23+
a,b = bc.args
24+
@assert bc.axes == axes(b)
25+
convert(AbstractArray{promote_type(T,V),N}, b)
26+
end
27+
28+
function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractArray{T,N},Ones{V,1,Tuple{OneToInfBlocks}}}}) where {N,T,V}
29+
a,b = bc.args
30+
@assert bc.axes == axes(a)
31+
convert(AbstractArray{promote_type(T,V),N}, a)
32+
end
333

4-
# for LazyLay in (:(BlockLayout{LazyLayout}), :(TriangularLayout{UPLO,UNIT,BlockLayout{LazyLayout}} where {UPLO,UNIT}))
5-
# @eval begin
6-
# combine_mul_styles(::$LazyLay) = LazyArrayApplyStyle()
7-
# mulapplystyle(::AbstractQLayout, ::$LazyLay) = LazyArrayApplyStyle()
8-
# end
9-
# end
34+
_block_interlace_axes(::Int, ax::Tuple{BlockedUnitRange{OneToInf{Int}}}...) = (blockedrange(Fill(length(ax), ∞)),)
1035

11-
# BlockArrays.blockbroadcaststyle(::LazyArrayStyle{N}) where N = LazyArrayStyle{N}()
36+
_block_interlace_axes(nbc::Int, ax::NTuple{2,BlockedUnitRange{OneToInf{Int}}}...) =
37+
(blockedrange(Fill(length(ax) ÷ nbc, ∞)),blockedrange(Fill(mod1(length(ax),nbc), ∞)))
1238

1339

1440
include("infblocktridiagonal.jl")
1541

42+
43+
#######
44+
# block broadcasted
45+
######
46+
47+
48+
BroadcastStyle(::Type{<:SubArray{T,N,Arr,<:NTuple{N,BlockSlice{BlockRange{1,Tuple{II}}}},false}}) where {T,N,Arr<:BlockArray,II<:InfRanges} =
49+
LazyArrayStyle{N}()
50+
51+
# TODO: generalise following
52+
for Ax in (:(RangeCumsum{Int,OneToInf{Int}}), :(OneToInf{Int}))
53+
@eval begin
54+
BroadcastStyle(::Type{BlockArray{T,N,Arr,NTuple{N,BlockedUnitRange{$Ax}}}}) where {T,N,Arr} = LazyArrayStyle{N}()
55+
BroadcastStyle(::Type{PseudoBlockArray{T,N,Arr,NTuple{N,BlockedUnitRange{$Ax}}}}) where {T,N,Arr} = LazyArrayStyle{N}()
56+
end
57+
end
58+
59+
BlockArrays._length(::BlockedUnitRange, ::OneToInf) = ∞
60+
BlockArrays._last(::BlockedUnitRange, ::OneToInf) = ∞
61+
62+
###
63+
# KronTrav
64+
###
65+
66+
_krontrav_axes(A::NTuple{N,OneToInf{Int}}, B::NTuple{N,OneToInf{Int}}) where N =
67+
@. blockedrange(OneTo(length(A)))

0 commit comments

Comments
 (0)