From 4f74de81822707c7aace96346415cbb35c30e4fb Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Fri, 10 Feb 2023 18:41:54 +0100 Subject: [PATCH 1/3] Adjust to `AdjointFactorization` --- Project.toml | 2 +- src/BandedMatrices.jl | 6 ++++++ src/banded/BandedLU.jl | 6 ++++-- src/banded/linalg.jl | 14 +++++++------- src/symbanded/SplitCholesky.jl | 12 +++++++++--- 5 files changed, 27 insertions(+), 13 deletions(-) diff --git a/Project.toml b/Project.toml index d0f37ba5..7a15d2ed 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "BandedMatrices" uuid = "aae01518-5342-5314-be14-df237901396f" -version = "0.17.11" +version = "0.17.12" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" diff --git a/src/BandedMatrices.jl b/src/BandedMatrices.jl index a9521c5d..769ab1ff 100644 --- a/src/BandedMatrices.jl +++ b/src/BandedMatrices.jl @@ -40,6 +40,12 @@ import FillArrays: AbstractFill, getindex_value, _broadcasted_zeros, unique_valu const libblas = Base.libblas_name const liblapack = Base.liblapack_name +const AdjointFact = isdefined(LinearAlgebra, :AdjointFactorization) ? + LinearAlgebra.AdjointFactorization : + Adjoint +const TransposeFact = isdefined(LinearAlgebra, :TransposeFactorization) ? + LinearAlgebra.TransposeFactorization : + Transpose export BandedMatrix, bandrange, diff --git a/src/banded/BandedLU.jl b/src/banded/BandedLU.jl index 9f8553be..97be3473 100644 --- a/src/banded/BandedLU.jl +++ b/src/banded/BandedLU.jl @@ -32,8 +32,10 @@ Base.iterate(S::BandedLU, ::Val{:U}) = (S.U, Val(:p)) Base.iterate(S::BandedLU, ::Val{:p}) = (S.p, Val(:done)) Base.iterate(S::BandedLU, ::Val{:done}) = nothing -adjoint(F::BandedLU) = Adjoint(F) -transpose(F::BandedLU) = Transpose(F) +if !(isdefined(LinearAlgebra, :AdjointFactorization)) # VERSION < v"1.10-" + adjoint(F::BandedLU) = Adjoint(F) + transpose(F::BandedLU) = Transpose(F) +end lu(S::BandedLU) = S diff --git a/src/banded/linalg.jl b/src/banded/linalg.jl index 36daeb39..1e154392 100644 --- a/src/banded/linalg.jl +++ b/src/banded/linalg.jl @@ -38,7 +38,7 @@ function ldiv!(A::BandedLU{T}, B::AbstractVecOrMat{Complex{T}}) where T<:Real B .= a .+ im.*b end -function ldiv!(transA::Transpose{T,<:BandedLU{T,<:BandedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasFloat} +function ldiv!(transA::TransposeFact{T,<:BandedLU{T,<:BandedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasFloat} A = transA.parent m = size(A.factors,1) l,u = bandwidths(A.factors) @@ -46,15 +46,15 @@ function ldiv!(transA::Transpose{T,<:BandedLU{T,<:BandedMatrix}}, B::StridedVecO LAPACK.gbtrs!('T', l, u-l, m, data, A.ipiv, B) end -function ldiv!(transA::Transpose{<:Any,<:BandedLU{<:Any,<:BandedMatrix}}, B::AbstractVecOrMat) +function ldiv!(transA::TransposeFact{<:Any,<:BandedLU{<:Any,<:BandedMatrix}}, B::AbstractVecOrMat) A = transA.parent ldiv!(transpose(UnitLowerTriangular(A.factors)), ldiv!(transpose(UpperTriangular(A.factors)), B)) _apply_inverse_ipiv_rows!(A, B) end -ldiv!(adjF::Adjoint{T,<:BandedLU{T,<:BandedMatrix}}, B::AbstractVecOrMat{T}) where {T<:Real} = +ldiv!(adjF::AdjointFact{T,<:BandedLU{T,<:BandedMatrix}}, B::AbstractVecOrMat{T}) where {T<:Real} = (F = adjF.parent; ldiv!(transpose(F), B)) -function ldiv!(adjA::Adjoint{T,<:BandedLU{T,<:BandedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} +function ldiv!(adjA::AdjointFact{T,<:BandedLU{T,<:BandedMatrix}}, B::StridedVecOrMat{T}) where {T<:BlasComplex} A = adjA.parent m = size(A.factors,1) l,u = bandwidths(A.factors) @@ -62,14 +62,14 @@ function ldiv!(adjA::Adjoint{T,<:BandedLU{T,<:BandedMatrix}}, B::StridedVecOrMat LAPACK.gbtrs!('C', l, u-l, m, data, A.ipiv, B) end -function ldiv!(adjA::Adjoint{<:Any,<:BandedLU{<:Any,<:BandedMatrix}}, B::AbstractVecOrMat) +function ldiv!(adjA::AdjointFact{<:Any,<:BandedLU{<:Any,<:BandedMatrix}}, B::AbstractVecOrMat) error("Implement") A = adjA.parent ldiv!(adjoint(UnitLowerTriangular(A.factors)), ldiv!(adjoint(UpperTriangular(A.factors)), B)) _apply_inverse_ipiv!(A, B) end -\(A::Adjoint{<:Any,<:BandedLU}, B::Adjoint{<:Any,<:AbstractVecOrMat}) = A \ copy(B) -\(A::Transpose{<:Any,<:BandedLU}, B::Transpose{<:Any,<:AbstractVecOrMat}) = A \ copy(B) +\(A::AdjointFact{<:Any,<:BandedLU}, B::Adjoint{<:Any,<:AbstractVecOrMat}) = A \ copy(B) +\(A::TransposeFact{<:Any,<:BandedLU}, B::Transpose{<:Any,<:AbstractVecOrMat}) = A \ copy(B) _factorize(::AbstractBandedLayout, _, A) = size(A,1) == size(A,2) ? lu(A) : qr(A) diff --git a/src/symbanded/SplitCholesky.jl b/src/symbanded/SplitCholesky.jl index 80b40342..4147d7e0 100644 --- a/src/symbanded/SplitCholesky.jl +++ b/src/symbanded/SplitCholesky.jl @@ -45,7 +45,13 @@ function splitcholesky!(::SymmetricLayout{<:BandedColumnMajor}, SplitCholesky(A, A.uplo) end -adjoint(S::SplitCholesky) = Adjoint(S) +if !(isdefined(LinearAlgebra, :AdjointFactorization)) # VERSION < v"1.10-" + adjoint(S::SplitCholesky) = Adjoint(S) +else + transpose(S::SplitCholesky{<:Real}) = S' + transpose(::SplitCholesky) = + throw(ArgumentError("transpose of SplitCholesky decomposition is not supported, consider using adjoint")) +end function lmul!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}} require_one_based_indexing(B) @@ -103,7 +109,7 @@ function ldiv!(S::SplitCholesky{T,Symmetric{T,M}}, B::AbstractVecOrMat{T}) where B end -function lmul!(S::Adjoint{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}} +function lmul!(S::AdjointFact{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}} require_one_based_indexing(B) n, nrhs = size(B, 1), size(B, 2) if size(S, 1) != n @@ -138,7 +144,7 @@ function lmul!(S::Adjoint{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVecOrMa B end -function ldiv!(S::Adjoint{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}} +function ldiv!(S::AdjointFact{T,SplitCholesky{T,Symmetric{T,M}}}, B::AbstractVecOrMat{T}) where {T,M<:BandedMatrix{T}} require_one_based_indexing(B) n, nrhs = size(B, 1), size(B, 2) if size(S, 1) != n From 1d938931221061c6be0971e0dbab18d984063c59 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 11 Feb 2023 17:58:27 +0100 Subject: [PATCH 2/3] make sure transpose of LU is Transpose* --- src/banded/BandedLU.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/banded/BandedLU.jl b/src/banded/BandedLU.jl index 97be3473..7e2c8fb3 100644 --- a/src/banded/BandedLU.jl +++ b/src/banded/BandedLU.jl @@ -35,6 +35,8 @@ Base.iterate(S::BandedLU, ::Val{:done}) = nothing if !(isdefined(LinearAlgebra, :AdjointFactorization)) # VERSION < v"1.10-" adjoint(F::BandedLU) = Adjoint(F) transpose(F::BandedLU) = Transpose(F) +else # overwrite fallback for transpose, use fallback for adjoint + transpose(F::BandedLU) = LinearAlgebra.TransposeFactorization(F) end lu(S::BandedLU) = S From a8098b4db34d576d98f40ec11fc5ceeb8984e80c Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sat, 11 Feb 2023 18:07:41 +0100 Subject: [PATCH 3/3] simplify --- src/banded/BandedLU.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/banded/BandedLU.jl b/src/banded/BandedLU.jl index 7e2c8fb3..e1fb991a 100644 --- a/src/banded/BandedLU.jl +++ b/src/banded/BandedLU.jl @@ -34,10 +34,8 @@ Base.iterate(S::BandedLU, ::Val{:done}) = nothing if !(isdefined(LinearAlgebra, :AdjointFactorization)) # VERSION < v"1.10-" adjoint(F::BandedLU) = Adjoint(F) - transpose(F::BandedLU) = Transpose(F) -else # overwrite fallback for transpose, use fallback for adjoint - transpose(F::BandedLU) = LinearAlgebra.TransposeFactorization(F) end +transpose(F::BandedLU) = TransposeFact(F) lu(S::BandedLU) = S