Make GivensQ subtype AbstractQ, use API

Project.toml

@@ -1,7 +1,7 @@
 name = "UpdatableQRFactorizations"
 uuid = "8b110c62-ada2-4dd7-b53a-29f29fe8f7f4"
 authors = ["Sebastian Ament <sebastianeament@gmail.com>"]
-version = "1.0.0"
+version = "1.0.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/UpdatableQRFactorizations.jl

@@ -1,12 +1,13 @@
 module UpdatableQRFactorizations
 
 using LinearAlgebra
-using LinearAlgebra: Givens
+using LinearAlgebra: Givens, AbstractQ
 
 export AbstractQR, GivensQ, GivensQR, UpdatableGivensQR, UpdatableQR, UQR,
         add_column!, remove_column!
 
 const AbstractMatOrFac{T} = Union{AbstractMatrix{T}, Factorization{T}}
+const AdjointQ = isdefined(LinearAlgebra, :AdjointQ) ? LinearAlgebra.AdjointQ : Adjoint
 
 # IDEA: could have efficient version for sparse matrices
 include("abstract.jl")

src/abstract.jl

@@ -3,15 +3,18 @@ abstract type AbstractQR{T} <: Factorization{T} end
 # AbstractQR assumes existence of n, m, Q, R fields
 Base.size(F::AbstractQR) = (F.n, F.m)
 Base.size(F::AbstractQR, i::Int) = i > 2 ? 1 : size(F)[i]
-Base.eltype(F::AbstractQR{T}) where T = T
+Base.eltype(::AbstractQR{T}) where T = T
 
 Base.:\(F::AbstractQR, x::AbstractVecOrMat) = ldiv!(F, copy(x))
+# disambiguation
+Base.:\(F::AbstractQR{T}, B::VecOrMat{Complex{T}}) where {T<:LinearAlgebra.BlasReal} =
+    complex.(F \ real(B), F \ imag(B))
 
 # some helpers
 """
-```
+
     number_of_rotations(n::Int, m::Int)
-```
+
 Computes the number of Givens rotations that are necessary to compute the QR
 factorization of a general matrix of size n by m.
 """
@@ -25,9 +28,9 @@ number_of_rotations_to_append_column(n::Int, m::Int, k::Int = 1) = k*(n-m) - (k*
 number_of_rotations_to_remove_column(m::Int, k::Int) = m-k
 
 """
-```
+
     allocate_rotations(T::DataType, n::Int, m::Int)
-```
+
 Allocates a vector of Givens rotations of a length that is necessary to compute
 the QR factorization of a general matrix of size n by m.
 """

src/givensQ.jl

@@ -1,5 +1,5 @@
 # Givens representation of Q matrix, memory efficient compared to dense representation
-struct GivensQ{T, QT<:AbstractVector{<:Givens{T}}} <: Factorization{T}
+struct GivensQ{T, QT<:AbstractVector{<:Givens{T}}} <: AbstractQ{T}
     rotations::QT
     n::Int
     m::Int
@@ -15,67 +15,104 @@ end
 Base.size(F::GivensQ) = (F.n, F.n) # could be square but also (n, m)?
 Base.size(F::GivensQ, i::Int) = i > 2 ? 1 : size(F)[i]
 Base.copy(F::GivensQ) = GivensQ(copy(F.rotations), F.n, F.m)
-Base.adjoint(F::GivensQ) = Adjoint(F)
+if !isdefined(LinearAlgebra, :AdjointQ) # VERSION < v"1.10-"
+    Base.adjoint(F::GivensQ) = Adjoint(F)
+end
 # NOTE: this also works if rotations is empty (is identity operator)
 function LinearAlgebra.lmul!(F::GivensQ, X::AbstractVecOrMat)
-    m = size(X, 1)
     for G in Iterators.reverse(F.rotations) # the adjoint of a GivensQ reverses the order of the rotations
         lmul!(G', X) # ... and takes their adjoint
     end
     return X
 end
-function LinearAlgebra.lmul!(A::Adjoint{<:Any, <:GivensQ}, X::AbstractVecOrMat)
-    F = A.parent
+function LinearAlgebra.lmul!(A::AdjointQ{<:Any, <:GivensQ}, X::AbstractVecOrMat)
+    F = parent(A)
     for G in F.rotations # the adjoint of a GivensQ reverses the order of the rotations
         lmul!(G, X) # ... and takes their adjoint
     end
     return X
 end
-LinearAlgebra.rmul!(X::AbstractVecOrMat, F::GivensQ) = lmul!(F', X')'
-LinearAlgebra.rmul!(X::AbstractVecOrMat, F::Adjoint{<:Any, <:GivensQ}) = lmul!(F', X')'
-
-lmul(F, X) = lmul!(F, copy(X))
-rmul(X, F) = rmul!(copy(X), F)
+# disambiguation
+LinearAlgebra.lmul!(F::GivensQ, X::LinearAlgebra.AbstractTriangular) = lmul!(F, full!(X))
+LinearAlgebra.lmul!(A::AdjointQ{<:Any, <:GivensQ}, X::LinearAlgebra.AbstractTriangular) =
+    lmul!(F, full!(X))
 
-Base.:*(F::GivensQ, X::AbstractVector) = lmul(F, X)
-Base.:*(F::GivensQ, X::AbstractMatrix) = lmul(F, X)
-Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::AbstractVector) = lmul(F, X)
-Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = lmul(F, X)
-
-Base.:*(X::AbstractMatrix, F::GivensQ) = rmul(X, F)
-Base.:*(X::AbstractMatrix, F::Adjoint{<:Any, <:GivensQ}) = rmul(X, F)
+LinearAlgebra.rmul!(X::AbstractVecOrMat, F::GivensQ) = lmul!(F', X')'
+LinearAlgebra.rmul!(X::AbstractVecOrMat, F::AdjointQ{<:Any, <:GivensQ}) = lmul!(F', X')'
+# disambiguation
+LinearAlgebra.rmul!(X::LinearAlgebra.AbstractTriangular, F::GivensQ) = lmul!(F', X')'
+LinearAlgebra.rmul!(X::LinearAlgebra.AbstractTriangular, F::AdjointQ{<:Any, <:GivensQ}) =
+    lmul!(F', X')'
 
-function Base.:*(F::GivensQ, G::Adjoint{<:Any, <:GivensQ})
+function Base.:*(F::GivensQ, G::AdjointQ{<:Any, <:GivensQ})
     T = promote_type(eltype(F), eltype(G))
     F === G' ? (one(T) * I)(F.n) : F * Matrix(G)
 end
-function Base.:*(F::Adjoint{<:Any, <:GivensQ}, G::GivensQ)
+function Base.:*(F::AdjointQ{<:Any, <:GivensQ}, G::GivensQ)
     T = promote_type(eltype(F), eltype(G))
     F' === G ? (one(T) * I)(G.n) : F * Matrix(G)
 end
+
 function Base.:(==)(F::GivensQ, G::GivensQ)
     F.rotations == G.rotations && F.n == G.n && F.m == G.m
 end
 
-LinearAlgebra.ldiv!(F::GivensQ, x::AbstractVector) = lmul!(F', x)
-LinearAlgebra.ldiv!(F::Adjoint{<:Any, <:GivensQ}, x::AbstractVector) = lmul!(F', x)
-LinearAlgebra.ldiv!(F::GivensQ, X::AbstractMatrix) = lmul!(F', X)
-LinearAlgebra.ldiv!(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = lmul!(F', X)
+if VERSION < v"1.10-"
+    lmul(F, X) = lmul!(F, copy(X))
+    rmul(X, F) = rmul!(copy(X), F)
 
+    Base.:*(F::GivensQ, X::AbstractVector) = lmul(F, X)
+    Base.:*(F::GivensQ, X::AbstractMatrix) = lmul(F, X)
+    Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::AbstractVector) = lmul(F, X)
+    Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = lmul(F, X)
 
-function LinearAlgebra.mul!(y::AbstractVector, F::GivensQ, x::AbstractVector)
-    @. y = x
-    lmul!(F, y)
-end
-function LinearAlgebra.mul!(Y::AbstractMatrix, F::GivensQ, X::AbstractMatrix)
-    @. Y = X
-    lmul!(F, Y)
-end
-function LinearAlgebra.mul!(y::AbstractVector, F::Adjoint{<:Any, <:GivensQ}, x::AbstractVector)
-    @. y = x
-    lmul!(F, y)
-end
-function LinearAlgebra.mul!(Y::AbstractMatrix, F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix)
-    @. Y = X
-    lmul!(F, Y)
+    Base.:*(X::AbstractMatrix, F::GivensQ) = rmul(X, F)
+    Base.:*(X::AbstractMatrix, F::AdjointQ{<:Any, <:GivensQ}) = rmul(X, F)
+
+    Base.:*(F::GivensQ, X::StridedVector) = lmul(F, X)
+    Base.:*(F::GivensQ, X::StridedMatrix) = lmul(F, X)
+    Base.:*(F::GivensQ, X::Adjoint{<:Any, <:StridedVecOrMat}) = lmul(F, X)
+    Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::StridedVector) = lmul(F, X)
+    Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::StridedMatrix) = lmul(F, X)
+    Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::Adjoint{<:Any, <:StridedVecOrMat}) = lmul(F, X)
+
+    Base.:*(X::StridedMatrix, F::GivensQ) = rmul(X, F)
+    Base.:*(X::StridedMatrix, F::AdjointQ{<:Any, <:GivensQ}) = rmul(X, F)
+
+    LinearAlgebra.mul!(Y::StridedVector, F::GivensQ, X::StridedVector) =
+        invoke(LinearAlgebra.mul!, Tuple{AbstractVector, GivensQ, AbstractVector}, Y, F, X)
+    LinearAlgebra.mul!(Y::StridedMatrix, F::GivensQ, X::StridedMatrix) =
+        invoke(LinearAlgebra.mul!, Tuple{AbstractMatrix, GivensQ, AbstractMatrix}, Y, F, X)
+    LinearAlgebra.mul!(Y::StridedVector, F::AdjointQ{<:Any, <:GivensQ}, X::StridedVector) =
+        invoke(LinearAlgebra.mul!, Tuple{AbstractVector, AdjointQ{<:Any, <:GivensQ}, AbstractVector}, Y, F, X)
+    LinearAlgebra.mul!(Y::StridedMatrix, F::AdjointQ{<:Any, <:GivensQ}, X::StridedMatrix) =
+        invoke(LinearAlgebra.mul!, Tuple{AbstractMatrix, AdjointQ{<:Any, <:GivensQ}, AbstractMatrix}, Y, F, X)
+
+    if VERSION < v"1.8-"
+        Base.:\(F::GivensQ, X::AbstractVector) = F'X
+        Base.:\(F::GivensQ, X::AbstractMatrix) = F'X
+        Base.:\(F::Adjoint{<:Any, <:GivensQ}, X::AbstractVector) = F'X
+        Base.:\(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = F'X
+    end
+    LinearAlgebra.ldiv!(F::GivensQ, x::AbstractVector) = lmul!(F', x)
+    LinearAlgebra.ldiv!(F::Adjoint{<:Any, <:GivensQ}, x::AbstractVector) = lmul!(F', x)
+    LinearAlgebra.ldiv!(F::GivensQ, X::AbstractMatrix) = lmul!(F', X)
+    LinearAlgebra.ldiv!(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = lmul!(F', X)
+
+    function LinearAlgebra.mul!(y::AbstractVector, F::GivensQ, x::AbstractVector)
+        @. y = x
+        lmul!(F, y)
+    end
+    function LinearAlgebra.mul!(Y::AbstractMatrix, F::GivensQ, X::AbstractMatrix)
+        @. Y = X
+        lmul!(F, Y)
+    end
+    function LinearAlgebra.mul!(y::AbstractVector, F::AdjointQ{<:Any, <:GivensQ}, x::AbstractVector)
+        @. y = x
+        lmul!(F, y)
+    end
+    function LinearAlgebra.mul!(Y::AbstractMatrix, F::AdjointQ{<:Any, <:GivensQ}, X::AbstractMatrix)
+        @. Y = X
+        lmul!(F, Y)
+    end    
 end

src/updatableQR.jl

@@ -1,7 +1,7 @@
 """
-```
+
     UpdatableGivensQR <: AbstractQR <: Factorization
-```
+
 Holds storage for Givens rotations and the R matrix for a QR factorization.
 Allows for efficient addition and deletion of columns.
 See `add_column!` and `remove_column!`.
@@ -102,9 +102,9 @@ function LinearAlgebra.ldiv!(F::UQR, x::AbstractVecOrMat)
 end
 
 """
-```
+
     add_column!(F::UpdatableGivensQR, x::AbstractVector, k::Int = size(F, 2) + 1)
-```
+
 Given an existing QR factorization `F` of a matrix, computes the factorization of
 the same matrix after a new column `x` has been added as the `k`ᵗʰ column.
 Computational complexity: O(nm) where size(F) = (n, m).
@@ -195,9 +195,9 @@ function ensure_space_to_append_column!(rotations::AbstractVector{<:Givens},
 end
 
 """
-```
+
 remove_column!(F::UpdatableGivensQR, k::Int = size(F, 2))
-```
+
 Updates the existing QR factorization `F` of a matrix to the factorization of
 the same matrix after its kᵗʰ column has been deleted.
 Computational complexity: O(m²) where size(F) = (n, m).

test/UpdatableQRFactorizations.jl

@@ -74,9 +74,11 @@ using UpdatableQRFactorizations
         x = randn(m)
         Ax = A*x
         @test F \ (A*x) ≈ x
+        @test F \ complex(A*x) ≈ x
         r = 4
         X = randn(m, r)
         @test F \ (A*X) ≈ X
+        @test F \ complex(A*X) ≈ X
     end
 
     @testset "UpdatableQR" begin