scaledmap.jl

"""
Lazy representation of a scaled map `λ * A = A * λ` with real or complex map
`A <: LinearMap{RealOrComplex}` and real or complex scaling factor
`λ <: RealOrComplex`.
"""
struct ScaledMap{T, S<:RealOrComplex, L<:LinearMap} <: LinearMap{T}
    λ::S
    lmap::L
    function ScaledMap{T}(λ::S, A::L) where {T, S <: RealOrComplex, L <: LinearMap{<:RealOrComplex}}
        Tprod = typeof(oneunit(S) * oneunit(eltype(A)))
        promote_type(T, Tprod) == T ||
           error("target type $T vs product of $S and $(eltype(A))")
        new{T,S,L}(λ, A)
    end
end

# constructor
ScaledMap(λ::RealOrComplex, lmap::LinearMap{<:RealOrComplex}) =
    ScaledMap{typeof(oneunit(λ) * oneunit(eltype(lmap)))}(λ, lmap)

# basic methods
Base.size(A::ScaledMap) = size(A.lmap)
Base.axes(A::ScaledMap) = axes(A.lmap)
Base.isreal(A::ScaledMap) = isreal(A.λ) && isreal(A.lmap)
LinearAlgebra.issymmetric(A::ScaledMap) = issymmetric(A.lmap)
LinearAlgebra.ishermitian(A::ScaledMap) = isreal(A.λ) && ishermitian(A.lmap)
LinearAlgebra.isposdef(A::ScaledMap) = isposdef(A.λ) && isposdef(A.lmap)

MulStyle(A::ScaledMap) = MulStyle(A.lmap)

Base.transpose(A::ScaledMap) = A.λ * transpose(A.lmap)
Base.adjoint(A::ScaledMap) = conj(A.λ) * adjoint(A.lmap)

# comparison (sufficient, not necessary)
Base.:(==)(A::ScaledMap, B::ScaledMap) =
    eltype(A) == eltype(B) && A.lmap == B.lmap && A.λ == B.λ

# scalar multiplication and division
Base.:(*)(α::RealOrComplex, A::LinearMap{<:RealOrComplex}) = ScaledMap(α, A)
Base.:(*)(A::LinearMap{<:RealOrComplex}, α::RealOrComplex) = ScaledMap(α, A)

Base.:(*)(α::Number, A::ScaledMap) = (α * A.λ) * A.lmap
Base.:(*)(A::ScaledMap, α::Number) = A.lmap * (A.λ * α)
# needed for disambiguation
Base.:(*)(α::RealOrComplex, A::ScaledMap) = (α * A.λ) * A.lmap
Base.:(*)(A::ScaledMap, α::RealOrComplex) = (A.λ * α) * A.lmap
Base.:(-)(A::LinearMap) = -1 * A

# composition (not essential, but might save multiple scaling operations)
Base.:(*)(A::ScaledMap, B::ScaledMap) = (A.λ * B.λ) * (A.lmap * B.lmap)
Base.:(*)(A::ScaledMap, B::LinearMap) = A.λ * (A.lmap * B)
Base.:(*)(A::LinearMap, B::ScaledMap) = (A * B.lmap) * B.λ

# multiplication with vectors/matrices
for In in (AbstractVector, AbstractMatrix)
    @eval begin
        # commutative case
        _unsafe_mul!(y, A::ScaledMap, x::$In{<:RealOrComplex}) =
            _unsafe_mul!(y, A.lmap, x, A.λ, false)
        _unsafe_mul!(y, A::ScaledMap, x::$In{<:RealOrComplex}, α, β) =
            _unsafe_mul!(y, A.lmap, x, A.λ * α, β)
        # non-commutative case
        _unsafe_mul!(y, A::ScaledMap, x::$In) = lmul!(A.λ, _unsafe_mul!(y, A.lmap, x))
    end
end

_unsafe_mul!(Y, X::ScaledMap, c::Number) = _unsafe_mul!(Y, X.lmap, X.λ*c)
_unsafe_mul!(Y, X::ScaledMap, c::Number, α, β) = _unsafe_mul!(Y, X.lmap, X.λ*c, α, β)