diff --git a/docs/make.jl b/docs/make.jl index 438b82e..90da4fc 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -15,15 +15,9 @@ makedocs(; linkcheck = false, checkdocs = :exports, format = Documenter.HTML(; - canonical = "https://docs.sciml.ai/FastAlmostBandedMatrices/stable/" + canonical = "https://docs.sciml.ai/FastAlmostBandedMatrices/stable/", ), - pages = [ - "Home" => "index.md", - "API" => "api.md", - ] + pages = ["Home" => "index.md", "API" => "api.md"], ) -deploydocs(; - repo = "github.com/SciML/FastAlmostBandedMatrices.jl.git", - push_preview = true -) +deploydocs(; repo = "github.com/SciML/FastAlmostBandedMatrices.jl.git", push_preview = true) diff --git a/src/FastAlmostBandedMatrices.jl b/src/FastAlmostBandedMatrices.jl index f01eb59..323033f 100644 --- a/src/FastAlmostBandedMatrices.jl +++ b/src/FastAlmostBandedMatrices.jl @@ -8,8 +8,18 @@ import ArrayLayouts: LayoutMatrix, LayoutVector, Ldiv, Lmul, TriangularLayout import ConcreteStructs: @concrete import LazyArrays: LazyArray, Mul import LinearAlgebra -import LinearAlgebra: LowerTriangular, NoPivot, UnitLowerTriangular, UnitUpperTriangular, - UpperTriangular, diagind, lmul!, lu, qr, rank, triu! +import LinearAlgebra: + LowerTriangular, + NoPivot, + UnitLowerTriangular, + UnitUpperTriangular, + UpperTriangular, + diagind, + lmul!, + lu, + qr, + rank, + triu! import MatrixFactorizations # The BandedMatrices.jl surface that FastAlmostBandedMatrices reexports (see the second @@ -17,8 +27,19 @@ import MatrixFactorizations # the `bands` argument of an `AlmostBandedMatrix`, populate it, and query its band # structure. The BandedMatrices.jl definitions remain canonical; the two docstrings below # document their public bindings in this module. -using BandedMatrices: Band, BandError, BandRange, BandedMatrix, band, bandrange, bandwidth, - bandwidths, brand, brandn, colrange, rowrange +using BandedMatrices: + Band, + BandError, + BandRange, + BandedMatrix, + band, + bandrange, + bandwidth, + bandwidths, + brand, + brandn, + colrange, + rowrange """ Band(i) @@ -81,9 +102,21 @@ true """ BandError -import ArrayLayouts: MemoryLayout, sublayout, MatLdivVec, materialize!, - triangularlayout, triangulardata, colsupport, - rowsupport, _qr, _qr!, _factorize, muladd!, QRPackedQLayout, AdjQRPackedQLayout +import ArrayLayouts: + MemoryLayout, + sublayout, + MatLdivVec, + materialize!, + triangularlayout, + triangulardata, + colsupport, + rowsupport, + _qr, + _qr!, + _factorize, + muladd!, + QRPackedQLayout, + AdjQRPackedQLayout import BandedMatrices: _banded_qr!, banded_qr_lmul! import LinearAlgebra: ldiv! import MatrixFactorizations: QR, QRPackedQ, getQ, getR @@ -98,8 +131,8 @@ import MatrixFactorizations: QR, QRPackedQ, getQ, getR A lazy representation of the union of two ranges, supporting iteration and indexing without heap allocation. """ -struct DisjointRange{T <: Integer, R1 <: AbstractUnitRange{T}, R2 <: AbstractUnitRange{T}} <: - AbstractVector{T} +struct DisjointRange{T<:Integer,R1<:AbstractUnitRange{T},R2<:AbstractUnitRange{T}} <: + AbstractVector{T} r1::R1 r2::R2 end @@ -113,7 +146,7 @@ Base.length(d::DisjointRange) = length(d.r1) + length(d.r2) if i <= n1 return @inbounds d.r1[i] else - return @inbounds d.r2[i - n1] + return @inbounds d.r2[i-n1] end end @@ -239,9 +272,12 @@ of columns, or if the lower bandwidth of `bands` is too small for the fill rank. end function AlmostBandedMatrix( - ::UndefInitializer, ::Type{T}, mn::NTuple{2, Integer}, - lu::NTuple{2, Integer}, rank::Integer - ) where {T} + ::UndefInitializer, + ::Type{T}, + mn::NTuple{2,Integer}, + lu::NTuple{2,Integer}, + rank::Integer, +) where {T} @assert lu[2] ≥ rank - 1 @assert rank ≥ 1 "Rank 0 fill array makes it a BandedMatrix." bands = BandedMatrix{T}(undef, mn, lu) @@ -250,17 +286,20 @@ function AlmostBandedMatrix( end function AlmostBandedMatrix{T}( - ::UndefInitializer, mn::NTuple{2, Integer}, - lu::NTuple{2, Integer}, rank::Integer - ) where {T} + ::UndefInitializer, + mn::NTuple{2,Integer}, + lu::NTuple{2,Integer}, + rank::Integer, +) where {T} return AlmostBandedMatrix(undef, T, mn, lu, rank) end function AlmostBandedMatrix( - ::UndefInitializer, mn::NTuple{2, Integer}, lu::NTuple{ - 2, Integer, - }, rank::Integer - ) + ::UndefInitializer, + mn::NTuple{2,Integer}, + lu::NTuple{2,Integer}, + rank::Integer, +) return AlmostBandedMatrix(undef, Float64, mn, lu, rank) end @@ -313,7 +352,7 @@ end @inline function finish_part_setindex!(bands, fill) # copy `fill` into `bands` in the correct locations l, u = bandwidths(bands) - for i in 1:size(fill, 1), j in max(1, i - l):min(size(bands, 2), i + u) + for i = 1:size(fill, 1), j = max(1, i-l):min(size(bands, 2), i+u) @inbounds bands[i, j] = fill[i, j] end @@ -389,7 +428,7 @@ E = exclusive_bandpart(A) # Returns a view of rows 3:10 of the banded part """ @inline function exclusive_bandpart(A) B, F = bandpart(A), fillpart(A) - return @view(B[(size(F, 1) + 1):end, :]) + return @view(B[(size(F, 1)+1):end, :]) end """ @@ -481,9 +520,9 @@ end @inline function rowsupport(::AbstractAlmostBandedLayout, A, k) l, _ = almostbandwidths(A) if maximum(k) ≤ almostbandedrank(A) - return max(1, minimum(k) - l):size(A, 2) + return max(1, minimum(k)-l):size(A, 2) else - return max(1, minimum(k) - l):min(maximum(k) + l, size(A, 2)) + return max(1, minimum(k)-l):min(maximum(k)+l, size(A, 2)) end end @@ -523,12 +562,9 @@ end # TODO: Support views properly function sublayout( - ::AlmostBandedLayout, ::Type{ - <:Tuple{ - AbstractUnitRange{Int}, AbstractUnitRange{Int}, - }, - } - ) + ::AlmostBandedLayout, + ::Type{<:Tuple{AbstractUnitRange{Int},AbstractUnitRange{Int}}}, +) return AlmostBandedLayout() end @@ -564,9 +600,10 @@ end # --------------- function _almost_banded_summary(io, B::AlmostBandedMatrix{T}, inds) where {T} return print( - io, Base.dims2string(length.(inds)), + io, + Base.dims2string(length.(inds)), " AlmostBandedMatrix{$T} with bandwidths $(almostbandwidths(B)) and fill \ - rank $(almostbandedrank(B))" + rank $(almostbandedrank(B))", ) end function Base.array_summary(io::IO, B::AlmostBandedMatrix, inds::Tuple{Vararg{Base.OneTo}}) @@ -591,7 +628,7 @@ end function ArrayInterface.fast_scalar_indexing(A::AlmostBandedMatrix) return ArrayInterface.fast_scalar_indexing(typeof(A.bands)) && - ArrayInterface.fast_scalar_indexing(typeof(A.fill)) + ArrayInterface.fast_scalar_indexing(typeof(A.fill)) end function ArrayInterface.qr_instance(A::AlmostBandedMatrix{T}, pivot = NoPivot()) where {T} @@ -613,7 +650,8 @@ function _almostbanded_qr(_, A) # Expand the bandsize for the QR factorization ## Bypass the safety checks in `AlmostBandedMatrix` return almostbanded_qr!( - AlmostBandedMatrix{eltype(A)}(BandedMatrix(copy(B), (l, l + u)), copy(L)), Val(true) + AlmostBandedMatrix{eltype(A)}(BandedMatrix(copy(B), (l, l + u)), copy(L)), + Val(true), ) end @@ -653,10 +691,10 @@ end k = 1 while k ≤ ncols - kr = k:min(k + l + u, m) - jr1 = k:min(k + u, n) - jr2 = (k + u + 1):min(last(kr) + u, n) - jr3 = k:min(k + u, n, ncols) + kr = k:min(k+l+u, m) + jr1 = k:min(k+u, n) + jr2 = (k+u+1):min(last(kr)+u, n) + jr3 = k:min(k+u, n, ncols) S = B[kr, jr1] τv = τ[jr3] R, _ = _banded_qr!(S, τv, length(jr3)) @@ -665,17 +703,13 @@ end B_right = B[kr, jr2] L_right = L[:, jr2] U′ = U[kr, :] - for j in 1:length(jr2) - muladd!( - -one(T), U′[(j + 1):end, :], L_right[:, j], one(T), B_right[(j + 1):end, j] - ) + for j = 1:length(jr2) + muladd!(-one(T), U′[(j+1):end, :], L_right[:, j], one(T), B_right[(j+1):end, j]) end banded_qr_lmul!(Q', B_right) banded_qr_lmul!(Q', U′) - for j in 1:length(jr2) - muladd!( - one(T), U′[(j + 1):end, :], L_right[:, j], one(T), B_right[(j + 1):end, j] - ) + for j = 1:length(jr2) + muladd!(one(T), U′[(j+1):end, :], L_right[:, j], one(T), B_right[(j+1):end, j]) end k = last(jr1) + 1 end @@ -683,10 +717,10 @@ end return AlmostBandedMatrix{eltype(A)}(B, Mul(U, L)), τ end -function getQ(F::QR{<:Any, <:AlmostBandedMatrix}) +function getQ(F::QR{<:Any,<:AlmostBandedMatrix}) return LinearAlgebra.QRPackedQ(bandpart(F.factors), F.τ) end -function getR(F::QR{<:Any, <:AlmostBandedMatrix}) +function getR(F::QR{<:Any,<:AlmostBandedMatrix}) n = min(size(F.factors, 1), size(F.factors, 2)) return UpperTriangular(view(F.factors, 1:n, 1:n)) end @@ -720,13 +754,13 @@ function _almostbanded_ldiv!(A::QR, B) end end -ldiv!(A::QR{T, <:AlmostBandedMatrix}, B::StridedVector{T}) where {T} = +ldiv!(A::QR{T,<:AlmostBandedMatrix}, B::StridedVector{T}) where {T} = _almostbanded_ldiv!(A, B) -ldiv!(A::QR{T, <:AlmostBandedMatrix}, B::StridedMatrix{T}) where {T} = +ldiv!(A::QR{T,<:AlmostBandedMatrix}, B::StridedMatrix{T}) where {T} = _almostbanded_ldiv!(A, B) -ldiv!(A::QR{T, <:AlmostBandedMatrix}, B::LayoutVector{T}) where {T} = +ldiv!(A::QR{T,<:AlmostBandedMatrix}, B::LayoutVector{T}) where {T} = _almostbanded_ldiv!(A, B) -ldiv!(A::QR{T, <:AlmostBandedMatrix}, B::LayoutMatrix{T}) where {T} = +ldiv!(A::QR{T,<:AlmostBandedMatrix}, B::LayoutMatrix{T}) where {T} = _almostbanded_ldiv!(A, B) # needed for adaptive QR @@ -738,7 +772,7 @@ function Base.materialize!(M::Lmul{<:AdjQRPackedQLayout{<:AlmostBandedLayout}}) return lmul!(QRPackedQ(bandpart(Q.factors), Q.τ)', M.B) end -triangularlayout(::Type{Tri}, ::ML) where {Tri, ML <: AlmostBandedLayout} = Tri{ML}() +triangularlayout(::Type{Tri}, ::ML) where {Tri,ML<:AlmostBandedLayout} = Tri{ML}() @inline function __arguments(x::LazyArray, ::AlmostBandedMatrix, ::Val) return LazyArrays.arguments(x) @@ -769,8 +803,11 @@ end @inline __original_almostbandedrank(A) = size(first(__lowrankfillpart(A)), 2) @views function _almostbanded_upper_ldiv!( - ::Type{Tri}, R::AbstractMatrix, b::AbstractVector{T}, buffer - ) where {T, Tri} + ::Type{Tri}, + R::AbstractMatrix, + b::AbstractVector{T}, + buffer, +) where {T,Tri} B = bandpart(R) U, V = __lowrankfillpart(R) fill!(buffer, zero(T)) @@ -779,9 +816,9 @@ end k = n = size(R, 2) while k > 0 - kr = max(1, k - u):k - jr1 = (k + 1):(k + u + 1) - jr2 = (k + u + 2):(k + 2u + 2) + kr = max(1, k-u):k + jr1 = (k+1):(k+u+1) + jr2 = (k+u+2):(k+2u+2) bv = b[kr] if jr2[1] < n muladd!(one(T), V[:, jr2], b[jr2], one(T), buffer) @@ -797,7 +834,7 @@ end return b end -function Base.materialize!(M::MatLdivVec{TriangularLayout{'U', 'N', AlmostBandedLayout}}) +function Base.materialize!(M::MatLdivVec{TriangularLayout{'U','N',AlmostBandedLayout}}) R, x = M.A, M.B A = triangulardata(R) r = __original_almostbandedrank(A) @@ -805,7 +842,7 @@ function Base.materialize!(M::MatLdivVec{TriangularLayout{'U', 'N', AlmostBanded return x end -function Base.materialize!(M::MatLdivVec{TriangularLayout{'U', 'U', AlmostBandedLayout}}) +function Base.materialize!(M::MatLdivVec{TriangularLayout{'U','U',AlmostBandedLayout}}) R, x = M.A, M.B A = triangulardata(R) r = __original_almostbandedrank(A) @@ -813,14 +850,14 @@ function Base.materialize!(M::MatLdivVec{TriangularLayout{'U', 'U', AlmostBanded return x end -function Base.materialize!(M::MatLdivVec{TriangularLayout{'L', 'N', AlmostBandedLayout}}) +function Base.materialize!(M::MatLdivVec{TriangularLayout{'L','N',AlmostBandedLayout}}) R, x = M.A, M.B A = triangulardata(R) materialize!(Ldiv(LowerTriangular(bandpart(A)), x)) return x end -function Base.materialize!(M::MatLdivVec{TriangularLayout{'L', 'U', AlmostBandedLayout}}) +function Base.materialize!(M::MatLdivVec{TriangularLayout{'L','U',AlmostBandedLayout}}) R, x = M.A, M.B A = triangulardata(R) materialize!(Ldiv(UnitLowerTriangular(bandpart(A)), x)) @@ -832,11 +869,15 @@ end # --------------- @views function muladd!( - α, A::AlmostBandedMatrix, B::AbstractVecOrMat, β, C::AbstractVecOrMat - ) + α, + A::AlmostBandedMatrix, + B::AbstractVecOrMat, + β, + C::AbstractVecOrMat, +) L = fillpart(A) muladd!(α, L, B, β, selectdim(C, 1, 1:size(L, 1))) - muladd!(α, exclusive_bandpart(A), B, β, selectdim(C, 1, (size(L, 1) + 1):size(C, 1))) + muladd!(α, exclusive_bandpart(A), B, β, selectdim(C, 1, (size(L, 1)+1):size(C, 1))) return C end @@ -863,14 +904,29 @@ end end end -export AlmostBandedMatrix, bandpart, fillpart, exclusive_bandpart, finish_part_setindex!, - almostbandwidths, almostbandedrank +export AlmostBandedMatrix, + bandpart, + fillpart, + exclusive_bandpart, + finish_part_setindex!, + almostbandwidths, + almostbandedrank # Reexported BandedMatrices.jl names; approved via `reexports_allow` in test/qa/qa.jl. # `AlmostBandedMatrix(bands::BandedMatrix, fill)` is the documented constructor and every # documented example builds `bands` with `brand`, so these must keep coming along with # `using FastAlmostBandedMatrices`. -export Band, BandError, BandRange, BandedMatrix, band, bandrange, bandwidth, bandwidths, - brand, brandn, colrange, rowrange +export Band, + BandError, + BandRange, + BandedMatrix, + band, + bandrange, + bandwidth, + bandwidths, + brand, + brandn, + colrange, + rowrange end diff --git a/test/core_tests.jl b/test/core_tests.jl index 8fc8035..e04f4b5 100644 --- a/test/core_tests.jl +++ b/test/core_tests.jl @@ -7,8 +7,18 @@ using SafeTestsets # Kept in sync with the reexport `export` block in src/FastAlmostBandedMatrices.jl, # `REEXPORTED_API` in test/qa/qa.jl, and the reexport section of docs/src/api.md. reexports = ( - :Band, :BandError, :BandRange, :BandedMatrix, :band, :bandrange, :bandwidth, - :bandwidths, :brand, :brandn, :colrange, :rowrange, + :Band, + :BandError, + :BandRange, + :BandedMatrix, + :band, + :bandrange, + :bandwidth, + :bandwidths, + :brand, + :brandn, + :colrange, + :rowrange, ) exported = Set(names(FastAlmostBandedMatrices)) @@ -32,7 +42,8 @@ using SafeTestsets # # Both files are read line-wise after normalising the line endings: git checks these # files out with CRLF on Windows, so nothing here may assume "\n". - readlines_lf(path) = split(replace(read(path, String), "\r\n" => "\n", "\r" => "\n"), '\n') + readlines_lf(path) = + split(replace(read(path, String), "\r\n" => "\n", "\r" => "\n"), '\n') qa_file = joinpath(@__DIR__, "qa", "qa.jl") if isfile(qa_file) @@ -45,7 +56,8 @@ using SafeTestsets declared = Set(Symbol(m[1]) for m in eachmatch(r":(\w+)", block[1])) declared == Set(reexports) || @error( "reexport sync: REEXPORTED_API disagrees with the exports", - file = qa_file, only_in_file = setdiff(declared, Set(reexports)), + file = qa_file, + only_in_file = setdiff(declared, Set(reexports)), only_in_exports = setdiff(Set(reexports), declared) ) @test declared == Set(reexports) @@ -65,7 +77,7 @@ using SafeTestsets @test start !== nothing if start !== nothing stop = findnext(line -> startswith(line, "## "), lines, start + 1) - section = lines[start:(stop === nothing ? lastindex(lines) : stop - 1)] + section = lines[start:(stop===nothing ? lastindex(lines) : stop-1)] # Only the first contiguous run of bullets under the heading is the list of # reexported names (" - Building the bands: `BandedMatrix`, ..."); the later # bullet list in the same section is the deliberate *exclusions*. @@ -77,15 +89,16 @@ using SafeTestsets end isempty(bullets) && @error( "reexport sync: no ` - ` role bullets under the heading", - file = docs_file, heading + file = docs_file, + heading ) @test !isempty(bullets) - documented = Set( - Symbol(m[1]) for line in bullets for m in eachmatch(r"`(\w+)`", line) - ) + documented = + Set(Symbol(m[1]) for line in bullets for m in eachmatch(r"`(\w+)`", line)) documented == Set(reexports) || @error( "reexport sync: the docs section disagrees with the exports", - file = docs_file, only_in_docs = setdiff(documented, Set(reexports)), + file = docs_file, + only_in_docs = setdiff(documented, Set(reexports)), only_in_exports = setdiff(Set(reexports), documented) ) @test documented == Set(reexports) @@ -162,8 +175,8 @@ end B, L = bandpart(A), fillpart(A) F = qr(A) - @test F.Q isa LinearAlgebra.QRPackedQ{Float64, <:BandedMatrix} - @test F.R isa UpperTriangular{Float64, <:SubArray{Float64, 2, <:AlmostBandedMatrix}} + @test F.Q isa LinearAlgebra.QRPackedQ{Float64,<:BandedMatrix} + @test F.R isa UpperTriangular{Float64,<:SubArray{Float64,2,<:AlmostBandedMatrix}} @test F.Q' * A ≈ F.R @test A == Ã @@ -188,8 +201,7 @@ end n = 80 A = AlmostBandedMatrix(BandedMatrix(fill(2.0, n, n), (1, 1)), fill(3.0, 1, n)) b = randn(n) - @test MemoryLayout(UpperTriangular(A)) == - TriangularLayout{'U', 'N', AlmostBandedLayout}() + @test MemoryLayout(UpperTriangular(A)) == TriangularLayout{'U','N',AlmostBandedLayout}() @test_broken UpperTriangular(Matrix(A)) \ b ≈ UpperTriangular(A) \ b @test_broken UnitUpperTriangular(Matrix(A)) \ b ≈ UnitUpperTriangular(A) \ b @test LowerTriangular(Matrix(A)) \ b ≈ LowerTriangular(A) \ b diff --git a/test/qa/qa.jl b/test/qa/qa.jl index 84a5f75..bcc9486 100644 --- a/test/qa/qa.jl +++ b/test/qa/qa.jl @@ -5,8 +5,18 @@ using SciMLTesting, FastAlmostBandedMatrices # the documented constructor and the documented examples build `bands` with `brand`. # Kept in sync with the reexport `export` block in src/FastAlmostBandedMatrices.jl. const REEXPORTED_API = ( - :Band, :BandError, :BandRange, :BandedMatrix, - :band, :bandrange, :bandwidth, :bandwidths, :brand, :brandn, :colrange, :rowrange, + :Band, + :BandError, + :BandRange, + :BandedMatrix, + :band, + :bandrange, + :bandwidth, + :bandwidths, + :brand, + :brandn, + :colrange, + :rowrange, ) run_qa( @@ -19,10 +29,21 @@ run_qa( # QR/QRPackedQ/getQ/getR, BandedMatrices _banded_qr!/banded_qr_lmul!. all_explicit_imports_are_public = (; ignore = ( - :MatLdivVec, :sublayout, :triangulardata, :triangularlayout, - :_qr, :_qr!, :_factorize, :QRPackedQLayout, :AdjQRPackedQLayout, - :QR, :QRPackedQ, :getQ, :getR, - :_banded_qr!, :banded_qr_lmul!, + :MatLdivVec, + :sublayout, + :triangulardata, + :triangularlayout, + :_qr, + :_qr!, + :_factorize, + :QRPackedQLayout, + :AdjQRPackedQLayout, + :QR, + :QRPackedQ, + :getQ, + :getR, + :_banded_qr!, + :banded_qr_lmul!, ), ), # Qualified accesses of non-public names: Base OneTo/array_summary/dims2string/ @@ -30,8 +51,15 @@ run_qa( # ArrayInterface fast_scalar_indexing/qr_instance. all_qualified_accesses_are_public = (; ignore = ( - :OneTo, :array_summary, :dims2string, :inds2string, :materialize!, - :QRPackedQ, :arguments, :fast_scalar_indexing, :qr_instance, + :OneTo, + :array_summary, + :dims2string, + :inds2string, + :materialize!, + :QRPackedQ, + :arguments, + :fast_scalar_indexing, + :qr_instance, ), ), ),