Skip to content

Commit 78c72fa

Browse files
dlfivefiftydlfivefifty
authored
update broadcasting (#44)
* update broadcasting * triangular ambiguity * Use MulAdd * Update Project.toml
1 parent ee31e15 commit 78c72fa

File tree

7 files changed

+80
-62
lines changed

7 files changed

+80
-62
lines changed

‎Project.toml‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,9 +19,9 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
1919
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
2020

2121
[compat]
22-
BandedMatrices = "0.11"
22+
BandedMatrices = "0.12"
2323
BlockArrays = "0.10"
2424
FillArrays = "0.7"
25-
LazyArrays = "0.11"
25+
LazyArrays = "0.12"
2626
MatrixFactorizations = "0.1,0.2"
2727
julia = "1"

‎src/BandedBlockBandedMatrix.jl‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -468,7 +468,10 @@ end
468468
######################################
469469
# BandedMatrix interface for Blocks #
470470
######################################
471-
@inline bandwidths(V::BandedBlockBandedBlock) = subblockbandwidths(parent(V))
471+
@inline function bandwidths(V::BandedBlockBandedBlock)
472+
inblockbands(V) && return subblockbandwidths(parent(V))
473+
(-720,-720)
474+
end
472475

473476

474477

@@ -477,7 +480,8 @@ blocks(V::BandedBlockBandedBlock)::Tuple{Int,Int} = Int(first(parentindices(V)).
477480
Int(last(parentindices(V)).block)
478481

479482

480-
function bandeddata(V::BandedBlockBandedBlock)
483+
function bandeddata(V::BandedBlockBandedBlock{T}) where T
484+
inblockbands(V) || return Array{T}(undef, 0, size(V,2))
481485
A = parent(V)
482486
u = A.u
483487
K_sl, J_sl = parentindices(V)

‎src/BlockBandedMatrices.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ import LazyArrays: AbstractStridedLayout, ColumnMajor, @lazymul, MatMulMatAdd, M
2626
triangularlayout, UpperTriangularLayout, TriangularLayout, MatLdivVec,
2727
triangulardata, subarraylayout, @lazyldiv, @lazylmul,
2828
AbstractColumnMajor, DenseColumnMajor, ColumnMajor,
29-
DiagonalLayout, apply!, materialize!, MulAdd
29+
DiagonalLayout, apply!, materialize!, MulAdd, mulapplystyle, MulAddStyle
3030

3131
import BlockArrays: BlockSizes, nblocks, blocksize, blockcheckbounds, global2blockindex,
3232
Block, BlockSlice, getblock, unblock, setblock!, globalrange,

‎src/broadcast.jl‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -357,3 +357,9 @@ function similar(bc::Broadcasted{BandedBlockBandedStyle, <:Any, typeof(+),
357357

358358
BandedBlockBandedMatrix{T}(undef, _combine_blocksizes(A,B), (max(Al,Bl), max(Au,Bu)), (max(Aλ,Bλ), max(Aμ,Bμ)))
359359
end
360+
361+
362+
####
363+
# Special case for Diagonal multiplicartion
364+
####
365+

‎src/linalg.jl‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -165,7 +165,8 @@ end
165165
similar(M::MulAdd{<:DiagonalLayout,<:AbstractBlockBandedLayout}, ::Type{T}) where T = similar(M.B,T)
166166
similar(M::MulAdd{<:AbstractBlockBandedLayout,<:DiagonalLayout}, ::Type{T}) where T = similar(M.A,T)
167167

168-
168+
mulapplystyle(::DiagonalLayout, ::AbstractBlockBandedLayout) = MulAddStyle()
169+
mulapplystyle(::AbstractBlockBandedLayout, ::DiagonalLayout) = MulAddStyle()
169170

170171

171172

‎src/triblockbanded.jl‎

Lines changed: 57 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -98,68 +98,71 @@ end
9898
end
9999

100100

101-
102-
@inline function materialize!(M::MatLdivVec{<:TriangularLayout{'U',UNIT,<:AbstractBlockBandedLayout},
103-
<:AbstractStridedLayout}) where UNIT
104-
U,dest = M.A,M.B
105-
T = eltype(dest)
106-
107-
A = triangulardata(U)
108-
@assert hasmatchingblocks(A)
109-
110-
@boundscheck size(A,1) == size(dest,1) || throw(BoundsError())
111-
112-
# impose block structure
113-
b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))
114-
115-
Bs = blocksizes(A)
116-
N = nblocks(Bs,1)
117-
118-
for K = N:-1:1
119-
b_2 = view(b, Block(K))
120-
Ũ = _triangular_matrix(Val('U'), Val(UNIT), view(A, Block(K,K)))
121-
apply!(\, Ũ, b_2)
122-
123-
if K ≥ 2
124-
KR = blockcolstart(A, K):Block(K-1)
125-
V_12 = view(A, KR, Block(K))
126-
b̃_1 = view(b, KR)
127-
b̃_1 .= (-one(T)).*Mul(V_12, b_2) .+ b̃_1
101+
for UNIT in ('U', 'N')
102+
@eval begin
103+
@inline function materialize!(M::MatLdivVec{<:TriangularLayout{'U',$UNIT,<:AbstractBlockBandedLayout},
104+
<:AbstractStridedLayout})
105+
U,dest = M.A,M.B
106+
T = eltype(dest)
107+
108+
A = triangulardata(U)
109+
@assert hasmatchingblocks(A)
110+
111+
@boundscheck size(A,1) == size(dest,1) || throw(BoundsError())
112+
113+
# impose block structure
114+
b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))
115+
116+
Bs = blocksizes(A)
117+
N = nblocks(Bs,1)
118+
119+
for K = N:-1:1
120+
b_2 = view(b, Block(K))
121+
Ũ = _triangular_matrix(Val('U'), Val($UNIT), view(A, Block(K,K)))
122+
apply!(\, Ũ, b_2)
123+
124+
if K ≥ 2
125+
KR = blockcolstart(A, K):Block(K-1)
126+
V_12 = view(A, KR, Block(K))
127+
b̃_1 = view(b, KR)
128+
materialize!(MulAdd(-one(T), V_12, b_2, one(T), b̃_1))
129+
end
130+
end
131+
132+
dest
128133
end
129-
end
130134

131-
dest
132-
end
135+
@inline function materialize!(M::MatLdivVec{<:TriangularLayout{'L',$UNIT,<:AbstractBlockBandedLayout},
136+
<:AbstractStridedLayout})
137+
L,dest = M.A, M.B
138+
T = eltype(dest)
139+
A = triangulardata(L)
140+
@assert hasmatchingblocks(A)
133141

134-
@inline function materialize!(M::MatLdivVec{<:TriangularLayout{'L',UNIT,<:AbstractBlockBandedLayout},
135-
<:AbstractStridedLayout}) where UNIT
136-
L,dest = M.A, M.B
137-
T = eltype(dest)
138-
A = triangulardata(L)
139-
@assert hasmatchingblocks(A)
142+
@boundscheck size(A,1) == size(dest,1) || throw(BoundsError())
140143

141-
@boundscheck size(A,1) == size(dest,1) || throw(BoundsError())
144+
# impose block structure
145+
b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))
142146

143-
# impose block structure
144-
b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))
147+
Bs = blocksizes(A)
148+
N = nblocks(Bs,1)
145149

146-
Bs = blocksizes(A)
147-
N = nblocks(Bs,1)
150+
for K = 1:N
151+
b_2 = view(b, Block(K))
152+
L̃ = _triangular_matrix(Val('L'), Val($UNIT), view(A, Block(K,K)))
153+
b_2 .= Ldiv(L̃, b_2)
148154

149-
for K = 1:N
150-
b_2 = view(b, Block(K))
151-
L̃ = _triangular_matrix(Val('L'), Val(UNIT), view(A, Block(K,K)))
152-
b_2 .= Ldiv(L̃, b_2)
153-
154-
if K < N
155-
KR = Block(K+1):blockcolstop(A, K)
156-
V_12 = view(A, KR, Block(K))
157-
b̃_1 = view(b, KR)
158-
b̃_1 .= (-one(T)).*Mul(V_12, b_2) .+ b̃_1
159-
end
160-
end
155+
if K < N
156+
KR = Block(K+1):blockcolstop(A, K)
157+
V_12 = view(A, KR, Block(K))
158+
b̃_1 = view(b, KR)
159+
materialize!(MulAdd(-one(T), V_12, b_2, one(T), b̃_1))
160+
end
161+
end
161162

162-
dest
163+
dest
164+
end
165+
end
163166
end
164167

165168

‎test/test_bandedblockbanded.jl‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
using BlockArrays, BandedMatrices, BlockBandedMatrices, FillArrays, SparseArrays, Test, LazyArrays , LinearAlgebra
2-
import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolrange, blockrowrange, colrange, rowrange, isbandedblockbanded
2+
import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolrange, blockrowrange, colrange, rowrange, isbandedblockbanded, bandeddata
33

44
@testset "BandedBlockBandedMatrix" begin
55
@testset "constructors" begin
@@ -190,7 +190,11 @@ import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolrange, blockrowran
190190
@test A[1,2] == view(A,Block(1,2))[1,1] == 11
191191
@test A[1,3] == view(A,Block(1,2))[1,2] == view(A,Block(1,2))[2] == 19
192192

193-
@test view(A, Block(3),Block(1)) ≈ [0,0,0]
193+
@test bandwidths(view(A, Block(3),Block(1))) == (-720,-720)
194+
@test isempty(bandeddata(view(A, Block(3),Block(1))))
195+
196+
@test A[Block(3,1)] == view(A, Block(3),Block(1)) == zeros(3,1)
197+
@test A[Block(3,1)] ≈ view(A, Block(3),Block(1)) ≈ zeros(3,1)
194198
@test_throws BandError view(A, Block(3),Block(1))[1,1] = 4
195199
@test_throws BoundsError view(A, Block(5,1))
196200
end

0 commit comments

Comments
 (0)