Prep for separation of AbstractMatrix & AbstractQ

Author: dkarrasch

src/LinearMaps.jl

@@ -1,11 +1,10 @@
 module LinearMaps
 
-export LinearMap
+export LinearMap, FunctionMap, FillMap, InverseMap
 export ⊗, squarekron, kronsum, ⊕, sumkronsum, khatrirao, facesplitting
-export FillMap
-export InverseMap
 
 using LinearAlgebra
+using LinearAlgebra: AbstractQ
 import LinearAlgebra: mul!
 using SparseArrays
 
@@ -18,8 +17,9 @@ using Base: require_one_based_indexing
 
 abstract type LinearMap{T} end
 
-const MapOrVecOrMat{T} = Union{LinearMap{T}, AbstractVecOrMat{T}}
-const MapOrMatrix{T} = Union{LinearMap{T}, AbstractMatrix{T}}
+const AbstractVecOrMatOrQ{T} = Union{AbstractVecOrMat{T}, AbstractQ{T}}
+const MapOrVecOrMat{T} = Union{LinearMap{T}, AbstractVecOrMatOrQ{T}}
+const MapOrMatrix{T} = Union{LinearMap{T}, AbstractMatrix{T}, AbstractQ{T}}
 const TransposeAbsVecOrMat{T} = Transpose{T,<:AbstractVecOrMat}
 const RealOrComplex = Union{Real, Complex}
 
@@ -31,7 +31,7 @@ Base.eltype(::LinearMap{T}) where {T} = T
 
 # conversion to LinearMap
 Base.convert(::Type{LinearMap}, A::LinearMap) = A
-Base.convert(::Type{LinearMap}, A::AbstractVecOrMat) = LinearMap(A)
+Base.convert(::Type{LinearMap}, A::AbstractVecOrMatOrQ) = LinearMap(A)
 
 convert_to_lmaps() = ()
 convert_to_lmaps(A) = (convert(LinearMap, A),)
@@ -49,6 +49,7 @@ MulStyle(::FiveArg, ::ThreeArg) = ThreeArg()
 MulStyle(::ThreeArg, ::ThreeArg) = ThreeArg()
 MulStyle(::LinearMap) = ThreeArg() # default
 MulStyle(::AbstractVecOrMat) = FiveArg()
+MulStyle(::AbstractQ) = ThreeArg()
 MulStyle(A::LinearMap, As::LinearMap...) = MulStyle(MulStyle(A), MulStyle(As...))
 
 Base.isreal(A::LinearMap) = eltype(A) <: Real
@@ -337,11 +338,11 @@ end
 include("transpose.jl") # transposing linear maps
 include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties
 include("left.jl") # left multiplication by a matrix/transpose or adjoint vector
+include("functionmap.jl") # using a function as linear map
 include("uniformscalingmap.jl") # the uniform scaling map, to be able to make linear combinations of LinearMap objects and multiples of I
 include("linearcombination.jl") # defining linear combinations of linear maps
 include("scaledmap.jl") # multiply by a (real or complex) scalar
 include("composition.jl") # composition of linear maps
-include("functionmap.jl") # using a function as linear map
 include("blockmap.jl") # block linear maps
 include("kronecker.jl") # Kronecker product of linear maps
 include("khatrirao.jl") # Khatri-Rao and face-splitting products
@@ -355,34 +356,36 @@ include("chainrules.jl") # AD rules through ChainRulesCore
 
 """
     LinearMap(A::LinearMap; kwargs...)::WrappedMap
-    LinearMap(A::AbstractVecOrMat; kwargs...)::WrappedMap
+    LinearMap(A::AbstractVecOrMatOrQ; kwargs...)::WrappedMap
     LinearMap(J::UniformScaling, M::Int)::UniformScalingMap
     LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...)::FunctionMap
     LinearMap(A::MapOrVecOrMat, dims::Dims{2}, index::NTuple{2, AbstractVector{Int}})::EmbeddedMap
     LinearMap(A::MapOrVecOrMat, dims::Dims{2}; offset::Dims{2})::EmbeddedMap
 
 Construct a linear map object, either
 
-1. from an existing `LinearMap` or `AbstractVecOrMat` `A`, with the purpose of redefining
-  its properties via the keyword arguments `kwargs`;
+1. from an existing `LinearMap` or `AbstractVecOrMat`/`AbstractQ` `A`, with the purpose of
+  redefining its properties via the keyword arguments `kwargs`, see below;
 2. a `UniformScaling` object `J` with specified (square) dimension `M`;
 3. from a function or callable object `f`;
-4. from an existing `LinearMap` or `AbstractVecOrMat` `A`, embedded in a larger zero map.
+4. from an existing `LinearMap` or `AbstractVecOrMat`/`AbstractQ` `A`, embedded in a larger
+   zero map.
 
 In the case of item 3, one also needs to specify the size of the equivalent matrix
-representation `(M, N)`, i.e., for functions `f` acting
-on length `N` vectors and producing length `M` vectors (with default value `N=M`).
-Preferably, also the `eltype` `T` of the corresponding matrix representation needs to be
-specified, i.e., whether the action of `f` on a vector will be similar to, e.g., multiplying
-by numbers of type `T`. If not specified, the devault value `T=Float64` will be assumed.
-Optionally, a corresponding function `fc` can be specified that implements the adjoint
-(=transpose in the real case) of `f`.
-
-The keyword arguments and their default values for the function-based constructor are:
-*   `issymmetric::Bool = false` : whether `A` or `f` act as a symmetric matrix
-*   `ishermitian::Bool = issymmetric & T<:Real` : whether `A` or `f` act as a Hermitian
-    matrix
-*   `isposdef::Bool = false` : whether `A` or `f` act as a positive definite matrix.
+representation `(M, N)`, i.e., for functions `f` acting on length `N` vectors and producing
+length `M` vectors (with default value `N=M`). Preferably, also the `eltype` `T` of the
+corresponding matrix representation needs to be specified, i.e., whether the action of `f`
+on a vector will be similar to, e.g., multiplying by numbers of type `T`. If not specified,
+the devault value `T=Float64` will be assumed. Optionally, a corresponding function `fc`
+can be specified that implements the adjoint (or transpose in the real case) of `f`.
+
+The keyword arguments and their default values are:
+
+* `issymmetric::Bool = false` : whether `A` or `f` act as a symmetric matrix
+* `ishermitian::Bool = issymmetric & T<:Real` : whether `A` or `f` act as a Hermitian
+  matrix
+* `isposdef::Bool = false` : whether `A` or `f` act as a positive definite matrix.
+
 For existing linear maps or matrices `A`, the default values will be taken by calling
 internal functions `_issymmetric`, `_ishermitian` and `_isposdef` on the existing object `A`.
 These in turn dispatch to (overloads of) `LinearAlgebra`'s `issymmetric`, `ishermitian`,
@@ -391,11 +394,11 @@ known at compile time as for certain structured matrices, but return `false` for
 `AbstractMatrix` types.
 
 For the function-based constructor, there is one more keyword argument:
-*   `ismutating::Bool` : flags whether the function acts as a mutating matrix multiplication
-    `f(y,x)` where the result vector `y` is the first argument (in case of `true`),
-    or as a normal matrix multiplication that is called as `y=f(x)` (in case of `false`).
-    The default value is guessed by looking at the number of arguments of the first
-    occurrence of `f` in the method table.
+* `ismutating::Bool` : flags whether the function acts as a mutating matrix multiplication
+  `f(y,x)` where the result vector `y` is the first argument (in case of `true`),
+  or as a normal matrix multiplication that is called as `y=f(x)` (in case of `false`).
+  The default value is guessed by looking at the number of arguments of the first
+  occurrence of `f` in the method table.
 
 For the `EmbeddedMap` constructors, `dims` specifies the total dimensions of the map. The
 `index` argument specifies two collections of indices `inds1` and `inds2`, such that for

src/blockmap.jl

@@ -63,7 +63,7 @@ Base.size(A::BlockMap) = (last(last(A.rowranges)), last(last(A.colranges)))
 # hcat
 ############
 """
-    hcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...)::BlockMap
+    hcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMatOrQ}...)::BlockMap
 
 Construct a (lazy) representation of the horizontal concatenation of the arguments.
 All arguments are promoted to `LinearMap`s automatically.
@@ -81,7 +81,7 @@ julia> L * ones(Int, 6)
  6
 ```
 """
-function Base.hcat(As::Union{LinearMap, UniformScaling, AbstractVecOrMat}...)
+function Base.hcat(As::Union{LinearMap, UniformScaling, AbstractVecOrMatOrQ}...)
     T = promote_type(map(eltype, As)...)
     nbc = length(As)
 
@@ -98,7 +98,7 @@ end
 # vcat
 ############
 """
-    vcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...)::BlockMap
+    vcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMatOrQ}...)::BlockMap
 
 Construct a (lazy) representation of the vertical concatenation of the arguments.
 All arguments are promoted to `LinearMap`s automatically.
@@ -119,7 +119,7 @@ julia> L * ones(Int, 3)
  3
 ```
 """
-function Base.vcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...)
+function Base.vcat(As::Union{LinearMap,UniformScaling,AbstractVecOrMatOrQ}...)
     T = promote_type(map(eltype, As)...)
     nbr = length(As)
 
@@ -137,7 +137,7 @@ end
 # hvcat
 ############
 """
-    hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling,AbstractVecOrMat}...)::BlockMap
+    hvcat(rows::Tuple{Vararg{Int}}, As::Union{LinearMap,UniformScaling,AbstractVecOrMatOrQ}...)::BlockMap
 
 Construct a (lazy) representation of the horizontal-vertical concatenation of the arguments.
 The first argument specifies the number of arguments to concatenate in each block row.
@@ -165,7 +165,7 @@ julia> L * ones(Int, 6)
 Base.hvcat
 
 function Base.hvcat(rows::Tuple{Vararg{Int}},
-                    As::Union{LinearMap, UniformScaling, AbstractVecOrMat}...)
+                    As::Union{LinearMap, UniformScaling, AbstractVecOrMatOrQ}...)
     nr = length(rows)
     T = promote_type(map(eltype, As)...)
     sum(rows) == length(As) ||
@@ -322,7 +322,7 @@ end
 # provide one global intermediate storage vector if necessary
 __blockmul!(::FiveArg, y, A, x::AbstractVecOrMat, α, β)  = ___blockmul!(y, A, x, α, β, nothing)
 __blockmul!(::ThreeArg, y, A, x::AbstractVecOrMat, α, β) = ___blockmul!(y, A, x, α, β, similar(y))
-function ___blockmul!(y, A, x::AbstractVecOrMat, α, β, ::Nothing)
+function ___blockmul!(y, A, x, α, β, ::Nothing)
     maps, rows, yinds, xinds = A.maps, A.rows, A.rowranges, A.colranges
     mapind = 0
     for (row, yi) in zip(rows, yinds)
@@ -336,7 +336,7 @@ function ___blockmul!(y, A, x::AbstractVecOrMat, α, β, ::Nothing)
     end
     return y
 end
-function ___blockmul!(y, A, x::AbstractVecOrMat, α, β, z)
+function ___blockmul!(y, A, x, α, β, z)
     maps, rows, yinds, xinds = A.maps, A.rows, A.rowranges, A.colranges
     mapind = 0
     for (row, yi) in zip(rows, yinds)
@@ -497,7 +497,7 @@ BlockDiagonalMap(maps::LinearMap...) =
 # since the below methods are more specific than the Base method,
 # they would redefine Base/SparseArrays behavior
 for k in 1:8 # is 8 sufficient?
-    Is = ntuple(n->:($(Symbol(:A, n))::AbstractVecOrMat), Val(k-1))
+    Is = ntuple(n->:($(Symbol(:A, n))::AbstractVecOrMatOrQ), Val(k-1))
     # yields (:A1, :A2, :A3, ..., :A(k-1))
     L = :($(Symbol(:A, k))::LinearMap)
     # yields :Ak
@@ -524,7 +524,7 @@ for k in 1:8 # is 8 sufficient?
 end
 
 """
-    blockdiag(As::Union{LinearMap,AbstractVecOrMat}...)::BlockDiagonalMap
+    blockdiag(As::Union{LinearMap,AbstractVecOrMatOrQ}...)::BlockDiagonalMap
 
 Construct a (lazy) representation of the diagonal concatenation of the arguments.
 To avoid fallback to the generic `SparseArrays.blockdiag`, there must be a `LinearMap`
@@ -533,7 +533,7 @@ object among the first 8 arguments.
 SparseArrays.blockdiag
 
 """
-    cat(As::Union{LinearMap,AbstractVecOrMat}...; dims=(1,2))::BlockDiagonalMap
+    cat(As::Union{LinearMap,AbstractVecOrMatOrQ}...; dims=(1,2))::BlockDiagonalMap
 
 Construct a (lazy) representation of the diagonal concatenation of the arguments.
 To avoid fallback to the generic `Base.cat`, there must be a `LinearMap`

src/show.jl

@@ -8,7 +8,7 @@ end
 Base.show(io::IO, A::LinearMap) = print(io, map_show(io, A, 0))
 
 map_show(io::IO, A::LinearMap, i) = ' '^i * map_summary(A) * _show(io, A, i)
-map_show(io::IO, A::AbstractVecOrMat, i) = ' '^i * summary(A)
+map_show(io::IO, A::AbstractVecOrMatOrQ, i) = ' '^i * summary(A)
 _show(io::IO, ::LinearMap, _) = ""
 function _show(io::IO, A::FunctionMap{T,F,Nothing}, _) where {T,F}
     "($(A.f); ismutating=$(A._ismutating), issymmetric=$(A._issymmetric), ishermitian=$(A._ishermitian), isposdef=$(A._isposdef))"

src/wrappedmap.jl

@@ -23,21 +23,24 @@ WrappedMap(lmap::MapOrVecOrMat{T}; kwargs...) where {T} = WrappedMap{T}(lmap; kw
 
 # cheap property checks (usually by type)
 _issymmetric(A::AbstractMatrix) = false
+_issymmetric(A::AbstractQ) = false
 _issymmetric(A::AbstractSparseMatrix) = issymmetric(A)
 _issymmetric(A::LinearMap) = issymmetric(A)
 _issymmetric(A::LinearAlgebra.RealHermSymComplexSym) = issymmetric(A)
 _issymmetric(A::Union{Bidiagonal,Diagonal,SymTridiagonal,Tridiagonal}) = issymmetric(A)
 
 _ishermitian(A::AbstractMatrix) = false
+_ishermitian(A::AbstractQ) = false
 _ishermitian(A::AbstractSparseMatrix) = ishermitian(A)
 _ishermitian(A::LinearMap) = ishermitian(A)
 _ishermitian(A::LinearAlgebra.RealHermSymComplexHerm) = ishermitian(A)
 _ishermitian(A::Union{Bidiagonal,Diagonal,SymTridiagonal,Tridiagonal}) = ishermitian(A)
 
 _isposdef(A::AbstractMatrix) = false
+_isposdef(A::AbstractQ) = false
 _isposdef(A::LinearMap) = isposdef(A)
 
-const VecOrMatMap{T} = WrappedMap{T,<:AbstractVecOrMat}
+const VecOrMatMap{T} = WrappedMap{T,<:AbstractVecOrMatOrQ}
 
 MulStyle(A::VecOrMatMap) = MulStyle(A.lmap)
 

test/wrappedmap.jl

@@ -113,4 +113,12 @@ using Test, LinearMaps, LinearAlgebra
         @test issymmetric(LinearMap(M)) == issymmetric(M)
         @test ishermitian(LinearMap(M)) == ishermitian(M)
     end
+
+    # AbstractQ
+    Q = qr(rand(10, 20)).Q
+    Qmap = @inferred LinearMap(Q)
+    @test size(Qmap) == (10, 10)
+    @test !issymmetric(Qmap)
+    @test !ishermitian(Qmap)
+    @test !isposdef(Qmap)
 end