Skip to content

Commit 7fbebff

Browse files
dlfivefiftydlfivefifty
authored
Support AbstractQLayout (#70)
* Revert "Revert "abstractq"" This reverts commit 31b1fce. * Revert "Revert "Merge branch 'master' of https://github.com/JuliaMatrices/BlockBandedMatrices.jl"" This reverts commit 08d6dd4. * Update BlockBandedMatrices.jl * Use MatrixFactorizations.QR * update qr * Update Project.toml
1 parent 31b1fce commit 7fbebff

11 files changed

+147
-73
lines changed

‎Project.toml

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "BlockBandedMatrices"
22
uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
3-
version = "0.8.1"
3+
version = "0.8.2"
44

55
[deps]
66
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -18,9 +18,9 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
1818
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
1919

2020
[compat]
21-
ArrayLayouts = "0.2.3"
22-
BandedMatrices = "0.15.2"
23-
BlockArrays = "0.12"
24-
FillArrays = "0.8"
25-
MatrixFactorizations = "0.3.1"
21+
ArrayLayouts = "0.2.6"
22+
BandedMatrices = "0.15.6"
23+
BlockArrays = "0.12.4"
24+
FillArrays = "0.8.8"
25+
MatrixFactorizations = "0.4"
2626
julia = "1.2"

‎src/AbstractBlockBandedMatrix.jl

Lines changed: 23 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,24 @@
1-
####
2-
# Matrix memory layout traits
3-
#
4-
# if MemoryLayout(A) returns BandedColumnMajor, you must override
5-
# pointer and leadingdimension
6-
# in addition to the banded matrix interface
7-
####
81

2+
"""
3+
AbstractBlockBandedLayout
4+
5+
isa a `MemoryLayout` that indicates that the array implements the block-banded
6+
interface.
7+
"""
98
abstract type AbstractBlockBandedLayout <: AbstractBlockLayout end
9+
10+
"""
11+
AbstractBandedBlockBandedLayout
12+
13+
isa a `MemoryLayout` that indicates that the array implements the banded-block-banded
14+
interface.
15+
"""
1016
abstract type AbstractBandedBlockBandedLayout <: AbstractBlockBandedLayout end
1117

18+
19+
struct BandedBlockBandedLayout <: AbstractBandedBlockBandedLayout end
20+
struct BlockBandedLayout <: AbstractBlockBandedLayout end
21+
1222
struct BandedBlockBandedColumns{LAY} <: AbstractBandedBlockBandedLayout end
1323
struct BandedBlockBandedRows{LAY} <: AbstractBandedBlockBandedLayout end
1424
struct BlockBandedColumns{LAY} <: AbstractBlockBandedLayout end
@@ -64,6 +74,12 @@ for Func in (:blockbanded_blockcolstart, :blockbanded_blockcolstop, :blockbanded
6474
@eval $Func(A, i::Block{1}) = $Func(A, Int(i))
6575
end
6676

77+
blockbanded_blockcolstart(A, i::BlockRange) = blockbanded_blockcolstart(A, minimum(i))
78+
blockbanded_blocktowstart(A, i::BlockRange) = blockbanded_blockrowstart(A, minimum(i))
79+
blockbanded_blockcolstop(A, i::BlockRange) = blockbanded_blockcolstop(A, maximum(i))
80+
blockbanded_blocktowstop(A, i::BlockRange) = blockbanded_blockrowstop(A, maximum(i))
81+
82+
6783
@inline blockcolsupport(::AbstractBlockBandedLayout, A, i) = blockbanded_blockcolstart(A,i):blockbanded_blockcolstop(A,i)
6884
@inline blockrowsupport(::AbstractBlockBandedLayout, A, i) = blockbanded_blockrowstart(A,i):blockbanded_blockrowstop(A,i)
6985

‎src/BandedBlockBandedMatrix.jl

Lines changed: 32 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -175,8 +175,10 @@ function BandedBlockBandedMatrix{T,Blocks,RR}(A::AbstractMatrix, axes::NTuple{2,
175175
end
176176

177177

178-
BandedBlockBandedMatrix(A::AbstractMatrix) =
179-
BandedBlockBandedMatrix(A, axes(A), blockbandwidths(A), subblockbandwidths(A))
178+
BandedBlockBandedMatrix{T}(A::AbstractMatrix) where T =
179+
copyto!(BandedBlockBandedMatrix{T}(undef, axes(A), blockbandwidths(A), subblockbandwidths(A)), A)
180+
181+
BandedBlockBandedMatrix(A::AbstractMatrix{T}) where T = BandedBlockBandedMatrix{T}(A)
180182

181183
BandedBlockBandedMatrix(A::AbstractMatrix, rdims::AbstractVector{Int}, cdims::AbstractVector{Int}, lu::NTuple{2,Int}, λμ::NTuple{2,Int}) =
182184
BandedBlockBandedMatrix(A, (blockedrange(rdims), blockedrange(cdims)), lu, λμ)
@@ -316,13 +318,18 @@ end
316318
end
317319

318320
## structured matrix methods ##
319-
function Base.replace_in_print_matrix(A::BandedBlockBandedMatrix, i::Integer, j::Integer, s::AbstractString)
321+
function _bandedblockbanded_replace_in_print_matrix(A, i, j, s)
320322
bi = findblockindex.(axes(A), (i,j))
321323
I,J = block.(bi)
322324
i,j = blockindex.(bi)
323-
-A.l ≤ Int(J-I) ≤ A.u && -A.λ ≤ j-i ≤ A.μ ? s : Base.replace_with_centered_mark(s)
325+
l,u = blockbandwidths(A)
326+
λ,μ = subblockbandwidths(A)
327+
-l ≤ Int(J-I) ≤ u && -λ ≤ j-i ≤ μ ? s : Base.replace_with_centered_mark(s)
324328
end
325329

330+
Base.replace_in_print_matrix(A::BandedBlockBandedMatrix, i::Integer, j::Integer, s::AbstractString) =
331+
_bandedblockbanded_replace_in_print_matrix(A, i, j, s)
332+
326333

327334
############
328335
# Indexing #
@@ -411,37 +418,48 @@ sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},Bloc
411418
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{Block1}}}) = BandedBlockBandedColumnMajor()
412419
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()
413420

421+
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{Block1}}}) = BandedLayout()
422+
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) = BandedLayout()
423+
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) = BandedLayout()
424+
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) = BandedLayout()
425+
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockRange1}}}) = BandedBlockBandedLayout()
426+
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockRange1}}}) = BandedBlockBandedLayout()
427+
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{Block1}}}) = BandedBlockBandedLayout()
428+
sublayout(::AbstractBandedBlockBandedLayout, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedLayout()
429+
430+
431+
sub_materialize(::AbstractBandedBlockBandedLayout, V, _) = BandedBlockBandedMatrix(V)
432+
414433
isbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedColumnMajor()
415434
isbandedblockbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedBlockBandedColumnMajor()
416435

417436

418-
subblockbandwidths(V::SubBandedBlockBandedMatrix) = subblockbandwidths(parent(V))
437+
subblockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1},<:BlockSlice{BlockRange1}}}) =
438+
subblockbandwidths(parent(V))
419439

420-
function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,Block1})
440+
function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1},<:BlockSlice{Block1}}})
421441
A = parent(V)
422-
423442
KR = parentindices(V)[1].block.indices[1]
424443
J = parentindices(V)[2].block
425444
shift = Int(KR[1])-Int(J)
426445
blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
427446
end
428447

429-
function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,Block1,BlockRange1})
448+
function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{<:Block1},<:BlockSlice{<:BlockRange1}}})
430449
A = parent(V)
431-
432450
K = parentindices(V)[1].block
433451
JR = parentindices(V)[2].block.indices[1]
434452
shift = Int(K)-Int(JR[1])
435-
436-
blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
453+
l,u = blockbandwidths(A)
454+
l - shift, u + shift
437455
end
438456

439-
function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,BlockRange1})
457+
function blockbandwidths(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{<:BlockRange1},<:BlockSlice{<:BlockRange1}}})
440458
A = parent(V)
441459

442460
KR = parentindices(V)[1].block.indices[1]
443461
JR = parentindices(V)[2].block.indices[1]
444-
shift = Int(KR[1])-Int(JR[1])
462+
shift = Int(first(KR))-Int(first(JR))
445463

446464
blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
447465
end
@@ -487,7 +505,7 @@ parentblocks2Int(V::BandedBlockBandedBlock)::Tuple{Int,Int} = Int(first(parentin
487505
######################################
488506
# BandedMatrix interface for Blocks #
489507
######################################
490-
@inline function bandwidths(V::BandedBlockBandedBlock)
508+
@inline function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice1,BlockSlice1}}) where T
491509
inblockbands(V) && return subblockbandwidths(parent(V))
492510
(-720,-720)
493511
end

‎src/BlockBandedMatrices.jl

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -14,19 +14,19 @@ import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Br
1414

1515
import LinearAlgebra: UniformScaling, isdiag, rmul!, lmul!, ldiv!, rdiv!,
1616
AbstractTriangular, AdjOrTrans, HermOrSym, StructuredMatrixStyle,
17-
qr, qr!, QRPackedQ
17+
qr, qr!
1818
import LinearAlgebra.BLAS: BlasInt, BlasFloat, @blasfunc, libblas, BlasComplex, BlasReal
1919
import LinearAlgebra.LAPACK: chktrans, chkdiag, liblapack, chklapackerror, checksquare, chkstride1,
2020
chkuplo
21-
import MatrixFactorizations: ql, ql!, _ql, QLPackedQ
21+
import MatrixFactorizations: ql, ql!, _ql, QLPackedQ, AdjQRPackedQLayout, AdjQLPackedQLayout, QR, QRPackedQ
2222
import SparseArrays: sparse
2323

2424
import ArrayLayouts: BlasMatLmulVec,
2525
triangularlayout, UpperTriangularLayout, TriangularLayout, MatLdivVec,
26-
triangulardata, sublayout,
26+
triangulardata, sublayout, sub_materialize,
2727
AbstractColumnMajor, DenseColumnMajor, ColumnMajor,
2828
DiagonalLayout, MulAdd, mul, colsupport, rowsupport,
29-
_qr, _factorize, _copyto!
29+
_qr, _factorize, _copyto!, zero!
3030

3131
import BlockArrays: blocksize, blockcheckbounds, BlockedUnitRange, blockisequal, DefaultBlockAxis,
3232
Block, BlockSlice, getblock, unblock, setblock!, block, blockindex,
@@ -38,7 +38,7 @@ import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrang
3838
inbands_setindex!, inbands_getindex, banded_setindex!,
3939
banded_generic_axpy!,
4040
BlasFloat, banded_dense_axpy!, MemoryLayout,
41-
BandedColumnMajor,
41+
BandedLayout, BandedColumnMajor,
4242
BandedSubBandedMatrix, bandeddata, tribandeddata,
4343
_BandedMatrix, colstart, colstop, rowstart, rowstop,
4444
BandedStyle, _fill_lmul!, bandshift,
@@ -50,7 +50,7 @@ export BandedBlockBandedMatrix, BlockBandedMatrix, BlockSkylineMatrix, blockband
5050

5151
const Block1 = Block{1,Int}
5252
const BlockRange1 = BlockRange{1,Tuple{UnitRange{Int}}}
53-
const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}
53+
const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}
5454

5555
blockcolrange(A...) = blockcolsupport(A...)
5656
blockrowrange(A...) = blockrowsupport(A...)

‎src/BlockSkylineMatrix.jl

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -65,11 +65,14 @@ BlockSkylineSizes(rows::AbstractVector{Int}, cols::AbstractVector{Int}, l::Abstr
6565
BlockSkylineSizes((blockedrange(rows),blockedrange(cols)), l, u)
6666

6767
BlockBandedSizes(b_axes::NTuple{2,AbstractUnitRange{Int}}, l::Int, u::Int) =
68-
BlockSkylineSizes(b_axes, Fill(l, blocklength(b_axes[2])), Fill(u, blocklength(b_axes[2])))
68+
convert(BlockBandedSizes, BlockSkylineSizes(b_axes, Fill(l, blocklength(b_axes[2])), Fill(u, blocklength(b_axes[2]))))
6969

7070

7171
BlockBandedSizes(rows::AbstractVector{Int}, cols::AbstractVector{Int}, l::Int, u::Int) =
72-
BlockSkylineSizes(rows, cols, Fill(l, length(cols)), Fill(u, length(cols)))
72+
convert(BlockBandedSizes, BlockSkylineSizes(rows, cols, Fill(l, length(cols)), Fill(u, length(cols))))
73+
74+
convert(::Type{BlockSkylineSizes{BS,LL,UU,BStarts,BStrides}}, b::BlockSkylineSizes) where {BS,LL,UU,BStarts,BStrides} =
75+
BlockSkylineSizes(convert(BS, b.axes), convert(BStarts, b.block_starts), convert(BStrides,b.block_strides), convert(LL,b.l), convert(UU,b.u))
7376

7477
colblockbandwidths(bs::BlockSkylineSizes) = (bs.l, bs.u)
7578

@@ -227,7 +230,7 @@ BlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
227230
rdims::AbstractVector{Int}, cdims::AbstractVector{Int},
228231
lu::NTuple{2,Int}) = BlockBandedMatrix{eltype(A)}(A, rdims, cdims, lu)
229232

230-
BlockBandedMatrix(A::BlockBandedMatrix, lu::NTuple{2,Int}) = BlockBandedMatrix(A, BlockBandedSizes(axes(A), lu...))
233+
BlockBandedMatrix(A::AbstractMatrix, lu::NTuple{2,Int}) = BlockBandedMatrix(A, BlockBandedSizes(axes(A), lu...))
231234

232235
function convert(::Type{BlockSkylineMatrix}, A::AbstractMatrix)
233236
@assert isblockbanded(A)

‎src/blockskylineqr.jl

Lines changed: 33 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -1,36 +1,47 @@
1-
qrf!(A,τ) = _qrf!(MemoryLayout(typeof(A)),MemoryLayout(typeof(τ)),A,τ)
1+
qrf!(A,τ) = _qrf!(MemoryLayout(A),MemoryLayout(τ),A,τ)
22
_qrf!(::AbstractColumnMajor,::AbstractStridedLayout,A::AbstractMatrix{T},τ::AbstractVector{T}) where T<:BlasFloat =
33
LAPACK.geqrf!(A,τ)
44
_apply_qr!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasReal =
55
LAPACK.ormqr!('L','T',A,τ,B)
66
_apply_qr!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasComplex =
77
LAPACK.ormqr!('L','C',A,τ,B)
8-
apply_qr!(A, τ, B) = _apply_qr!(MemoryLayout(typeof(A)), MemoryLayout(typeof(τ)), MemoryLayout(typeof(B)), A, τ, B)
8+
apply_qr!(A, τ, B) = _apply_qr!(MemoryLayout(A), MemoryLayout(τ), MemoryLayout(B), A, τ, B)
99

10-
qlf!(A,τ) = _qlf!(MemoryLayout(typeof(A)),MemoryLayout(typeof(τ)),A,τ)
10+
qlf!(A,τ) = _qlf!(MemoryLayout(A),MemoryLayout(τ),A,τ)
1111
_qlf!(::AbstractColumnMajor,::AbstractStridedLayout,A::AbstractMatrix{T},τ::AbstractVector{T}) where T<:BlasFloat =
1212
LAPACK.geqlf!(A,τ)
1313
_apply_ql!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasReal =
1414
LAPACK.ormql!('L','T',A,τ,B)
1515
_apply_ql!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasComplex =
1616
LAPACK.ormql!('L','C',A,τ,B)
17-
apply_ql!(A, τ, B) = _apply_ql!(MemoryLayout(typeof(A)), MemoryLayout(typeof(τ)), MemoryLayout(typeof(B)), A, τ, B)
17+
apply_ql!(A, τ, B) = _apply_ql!(MemoryLayout(A), MemoryLayout(τ), MemoryLayout(B), A, τ, B)
1818

19-
function qr!(A::BlockBandedMatrix{T}) where T
19+
function _blockbanded_qr!(A::AbstractMatrix, τ::AbstractVector)
20+
M,N = Block.(blocksize(A))
21+
(M < N ? axes(A,1) : axes(A,2)) == axes(τ,1) || throw(DimensionMismatch(""))
22+
_blockbanded_qr!(A, τ, min(N,M))
23+
QR(A,τ.blocks)
24+
end
25+
26+
function _blockbanded_qr!(A::AbstractMatrix, τ::AbstractVector, NCOLS::Block{1})
2027
l,u = blockbandwidths(A)
21-
M,N = blocksize(A)
22-
ax1 = M < N ? axes(A,1) : axes(A,2)
23-
τ = PseudoBlockVector{T}(undef, (ax1,))
24-
for K = 1:min(N,M)
25-
KR = Block.(K:min(K+l,M))
26-
V = view(A,KR,Block(K))
27-
t = view(τ,Block(K))
28+
M,N = Block.(blocksize(A))
29+
for K = Block(1):NCOLS
30+
KR = K:min(K+l,M)
31+
V = view(A,KR,K)
32+
t = view(τ,K)
2833
qrf!(V,t)
2934
for J = K+1:min(K+u,N)
30-
apply_qr!(V, t, view(A,KR,Block(J)))
35+
apply_qr!(V, t, view(A,KR,J))
3136
end
3237
end
33-
QR(A,τ.blocks)
38+
A,τ
39+
end
40+
41+
function qr!(A::BlockBandedMatrix{T}) where T
42+
M,N = blocksize(A)
43+
ax1 = M < N ? axes(A,1) : axes(A,2)
44+
_blockbanded_qr!(A, PseudoBlockVector{T}(undef, (ax1,)))
3445
end
3546

3647
function ql!(A::BlockBandedMatrix{T}) where T
@@ -57,11 +68,13 @@ function ql!(A::BlockBandedMatrix{T}) where T
5768
QL(A,τ.blocks)
5869
end
5970

60-
_qr(::AbstractBlockBandedLayout, _, A) = qr!(BlockBandedMatrix(A, (blockbandwidth(A,1), blockbandwidth(A,1)+blockbandwidth(A,2))))
71+
_qr(::AbstractBlockBandedLayout, ::Tuple{Integer,Integer}, A) = qr!(BlockBandedMatrix(A, (blockbandwidth(A,1), blockbandwidth(A,1)+blockbandwidth(A,2))))
72+
_qr(lay::AbstractBlockBandedLayout, ax::Tuple{AbstractUnitRange{Int},AbstractUnitRange{Int}}, A) = _qr(lay, map(length, ax), A)
6173
_ql(::AbstractBlockBandedLayout, _, A) = ql!(BlockBandedMatrix(A, (blockbandwidth(A,1)+blockbandwidth(A,2),blockbandwidth(A,2))))
6274
_factorize(::AbstractBlockBandedLayout, _, A) = qr(A)
6375

64-
function lmul!(adjQ::Adjoint{<:Any,<:QRPackedQ{<:Any,<:BlockSkylineMatrix}}, Bin::AbstractVector)
76+
function materialize!(Mul::Lmul{<:AdjQRPackedQLayout{<:AbstractBlockBandedLayout}})
77+
adjQ,Bin = Mul.A,Mul.B
6578
Q = parent(adjQ)
6679
A = Q.factors
6780
l,u = blockbandwidths(A)
@@ -78,23 +91,22 @@ function lmul!(adjQ::Adjoint{<:Any,<:QRPackedQ{<:Any,<:BlockSkylineMatrix}}, Bin
7891
end
7992
Bin
8093
end
81-
82-
function lmul!(adjQ::Adjoint{<:Any,<:QLPackedQ{<:Any,<:BlockSkylineMatrix}}, Bin::AbstractVector)
83-
Q = parent(adjQ)
94+
function materialize!(Mul::Lmul{<:AdjQLPackedQLayout{<:AbstractBlockBandedLayout}})
95+
Q = parent(Mul.A)
8496
A = Q.factors
8597
l,u = blockbandwidths(A)
8698
N,M = blocksize(A)
8799
# impose block structure
88100
ax1 = (axes(A,1),)
89101
τ = PseudoBlockArray(Q.τ, ax1)
90-
B = PseudoBlockArray(Bin, ax1)
102+
B = PseudoBlockArray(Mul.B, ax1)
91103
for K = N:-1:1
92104
KR = Block.(max(1,K-u):K)
93105
V = view(A,KR,Block(K))
94106
t = view(τ,Block(K))
95107
apply_ql!(V, t, view(B,KR))
96108
end
97-
Bin
109+
Mul.B
98110
end
99111

100112
# avoid LinearALgebra Strided obsession

‎src/broadcast.jl

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,7 @@ similar(bc::Broadcasted{<:AbstractBlockBandedStyle}, ::Type{T}) where T =
4646

4747
checkblocks(A, B) = blockisequal(axes(A), axes(B)) || throw(DimensionMismatch("*"))
4848

49-
function _blockbanded_copyto!(dest::AbstractMatrix{T}, src::AbstractMatrix) where T
49+
function _blockbanded_copyto!(dest::AbstractMatrix, src::AbstractMatrix)
5050
@boundscheck checkblocks(dest, src)
5151

5252
dl, du = colblockbandwidths(dest)
@@ -57,13 +57,13 @@ function _blockbanded_copyto!(dest::AbstractMatrix{T}, src::AbstractMatrix) wher
5757

5858
for J = 1:N
5959
for K = max(1,J-du[J]):min(J-su[J]-1,M)
60-
view(dest,Block(K),Block(J)) .= zero(T)
60+
zero!(view(dest,Block(K),Block(J)))
6161
end
6262
for K = max(1,J-su[J]):min(J+sl[J],M)
63-
view(dest,Block(K),Block(J)) .= view(src,Block(K),Block(J))
63+
copyto!(view(dest,Block(K),Block(J)), view(src,Block(K),Block(J)))
6464
end
6565
for K = max(1,J+sl[J]+1):min(J+dl[J],M)
66-
view(dest,Block(K),Block(J)) .= zero(T)
66+
zero!(view(dest,Block(K),Block(J)))
6767
end
6868
end
6969
dest

‎src/linalg.jl

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -77,8 +77,9 @@ function blockbandwidths(V::SubBlockSkylineMatrix{<:Any,LL,UU,BlockRange1,BlockR
7777
KR = parentindices(V)[1].block.indices[1]
7878
JR = parentindices(V)[2].block.indices[1]
7979

80-
@assert KR[1] == JR[1] == 1
81-
blockbandwidths(A)
80+
shift = first(KR) - first(JR)
81+
l,u = blockbandwidths(A)
82+
l-shift, u+shift
8283
end
8384

8485

@@ -92,6 +93,10 @@ sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:
9293
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
9394
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockIndexRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
9495

96+
isblockbanded(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockRange1}}}) =
97+
isblockbanded(parent(V))
98+
99+
sub_materialize(::AbstractBlockBandedLayout, V, _) = BlockBandedMatrix(V)
95100

96101
strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,<:Union{BlockRange1,Block1},Block1}) where {LL,UU} =
97102
(1,parent(V).block_sizes.block_strides[Int(parentindices(V)[2].block)])
@@ -133,5 +138,3 @@ strides(V::SubBlockSkylineMatrix{T,LL,UU,BlockIndexRange1,BlockIndexRange1}) whe
133138
(1,parent(V).block_sizes.block_strides[Int(parentindices(V)[2].block.block)])
134139

135140
MemoryLayout(V::SubBlockSkylineMatrix{T,LL,UU,BlockIndexRange1,BlockIndexRange1}) where {T,LL,UU} = ColumnMajor()
136-
137-

0 commit comments

Comments
 (0)