diff --git a/Project.toml b/Project.toml index 8eb6c2f7..627277fa 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "LinearMaps" uuid = "7a12625a-238d-50fd-b39a-03d52299707e" -version = "2.5.2" +version = "2.6" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/README.md b/README.md index e355f49d..64c1948c 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,9 @@ [![Coverage Status](https://coveralls.io/repos/github/Jutho/LinearMaps.jl/badge.svg?branch=master)](https://coveralls.io/github/Jutho/LinearMaps.jl?branch=master) [![codecov.io](http://codecov.io/github/Jutho/LinearMaps.jl/coverage.svg?branch=master)](http://codecov.io/github/Jutho/LinearMaps.jl?branch=master) -A Julia package for defining and working with linear maps, also known as linear transformations or linear operators acting on vectors. The only requirement for a LinearMap is that it can act on a vector (by multiplication) efficiently. +A Julia package for defining and working with linear maps, also known as linear +transformations or linear operators acting on vectors. The only requirement for +a LinearMap is that it can act on a vector (by multiplication) efficiently. ## What's new in v2.6 * New feature: "lazy" Kronecker product, Kronecker sums, and powers thereof @@ -29,91 +31,205 @@ A Julia package for defining and working with linear maps, also known as linear ## What's new in v2.2.0 * Fully Julia v0.7/v1.0 compatible. -* A `convert(SparseMatrixCSC, A::LinearMap)` function, that calls the `sparse` matrix generating function. +* A `convert(SparseMatrixCSC, A::LinearMap)` function, that calls the `sparse` + matrix generating function. ## What's new in v2.1.0 -* Fully Julia v0.7 compatible; dropped compatibility for previous versions of Julia from LinearMaps.jl v2.0.0 on. -* A 5-argument version for `mul!(y, A::LinearMap, x, α=1, β=0)`, which computes `y := α * A * x + β * y` and implements the usual 3-argument `mul!(y, A, x)` for the default `α` and `β`. -* Synonymous `convert(Matrix, A::LinearMap)` and `convert(Array, A::LinearMap)` functions, that call the `Matrix` constructor and return the matrix representation of `A`. +* Fully Julia v0.7 compatible; dropped compatibility for previous versions of + Julia from LinearMaps.jl v2.0.0 on. +* A 5-argument version for `mul!(y, A::LinearMap, x, α=1, β=0)`, which + computes `y := α * A * x + β * y` and implements the usual 3-argument + `mul!(y, A, x)` for the default `α` and `β`. +* Synonymous `convert(Matrix, A::LinearMap)` and `convert(Array, A::LinearMap)` + functions, that call the `Matrix` constructor and return the matrix + representation of `A`. * Multiplication with matrices, interpreted as a block row vector of vectors: - * `mul!(Y::AbstractArray, A::LinearMap, X::AbstractArray, α=1, β=0)`: applies `A` to each column of `X` and stores the result in-place in the corresponding column of `Y`; - * for the out-of-place multiplication, the approach is to compute `convert(Matrix, A * X)`; this is equivalent to applying `A` to each column of `X`. In generic code which handles both `A::AbstractMatrix` and `A::LinearMap`, the additional call to `convert` is a noop when `A` is a matrix. -* Full compatibility with [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl)'s `eigs` and `svds`; previously only `eigs` was working. For more, nicely collaborating packages see the [Example](#example) section. + * `mul!(Y::AbstractArray, A::LinearMap, X::AbstractArray, α=1, β=0)`: + applies `A` to each column of `X` and stores the result in-place in the + corresponding column of `Y`; + * for the out-of-place multiplication, the approach is to compute + `convert(Matrix, A * X)`; this is equivalent to applying `A` to each + column of `X`. In generic code which handles both `A::AbstractMatrix` and + `A::LinearMap`, the additional call to `convert` is a noop when `A` is a + matrix. +* Full compatibility with [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl)'s + `eigs` and `svds`; previously only `eigs` was working. For more, nicely + collaborating packages see the [Example](#example) section. ## Installation -Install with the package manager, i.e., `]add LinearMaps` (or `Pkg.add("LinearMaps")` in Julia versions below 0.7). +Install with the package manager, i.e., `]add LinearMaps` (or `Pkg.add("LinearMaps")` +in Julia versions below 0.7). ## Philosophy -Several iterative linear algebra methods such as linear solvers or eigensolvers only require an efficient evaluation of the matrix vector product, where the concept of a matrix can be formalized / generalized to a linear map (or linear operator in the special case of a square matrix). +Several iterative linear algebra methods such as linear solvers or eigensolvers +only require an efficient evaluation of the matrix vector product, where the +concept of a matrix can be formalized / generalized to a linear map (or linear +operator in the special case of a square matrix). The LinearMaps package provides the following functionality: -1. A `LinearMap` type that shares with the `AbstractMatrix` type that it responds to the functions `size`, `eltype`, `isreal`, `issymmetric`, `ishermitian` and `isposdef`, `transpose` and `adjoint` and multiplication with a vector using both `*` or the in-place version `mul!`. Linear algebra functions that use duck-typing for their arguments can handle `LinearMap` objects similar to `AbstractMatrix` objects, provided that they can be written using the above methods. Unlike `AbstractMatrix` types, `LinearMap` objects cannot be indexed, neither using `getindex` or `setindex!`. - -2. A single method `LinearMap` function that acts as a general purpose constructor (though it is only an abstract type) and allows to construct linear map objects from functions, or to wrap objects of type `AbstractMatrix` or `LinearMap`. The latter functionality is useful to (re)define the properties (`isreal`, `issymmetric`, `ishermitian`, `isposdef`) of the existing matrix or linear map. - -3. A framework for combining objects of type `LinearMap` and of type `AbstractMatrix` using linear combinations, transposition and composition, where the linear map resulting from these operations is never explicitly evaluated but only its matrix vector product is defined (i.e. lazy evaluation). The matrix vector product is written to minimize memory allocation by using a minimal number of temporary vectors. There is full support for the in-place version `mul!`, which should be preferred for higher efficiency in critical algorithms. In addition, it tries to recognize the properties of combinations of linear maps. In particular, compositions such as `A'*A` for arbitrary `A` or even `A'*B*C*B'*A` with arbitrary `A` and `B` and positive definite `C` are recognized as being positive definite and hermitian. In case a certain property of the resulting `LinearMap` object is not correctly inferred, the `LinearMap` method can be called to redefine the properties. +1. A `LinearMap` type that shares with the `AbstractMatrix` type that it + responds to the functions `size`, `eltype`, `isreal`, `issymmetric`, + `ishermitian` and `isposdef`, `transpose` and `adjoint` and multiplication + with a vector using both `*` or the in-place version `mul!`. Linear algebra + functions that use duck-typing for their arguments can handle `LinearMap` + objects similar to `AbstractMatrix` objects, provided that they can be + written using the above methods. Unlike `AbstractMatrix` types, `LinearMap` + objects cannot be indexed, neither using `getindex` or `setindex!`. + +2. A single method `LinearMap` function that acts as a general purpose + constructor (though it is only an abstract type) and allows to construct + linear map objects from functions, or to wrap objects of type + `AbstractMatrix` or `LinearMap`. The latter functionality is useful to + (re)define the properties (`isreal`, `issymmetric`, `ishermitian`, + `isposdef`) of the existing matrix or linear map. + +3. A framework for combining objects of type `LinearMap` and of type + `AbstractMatrix` using linear combinations, transposition and composition, + where the linear map resulting from these operations is never explicitly + evaluated but only its matrix vector product is defined (i.e. lazy + evaluation). The matrix vector product is written to minimize memory + allocation by using a minimal number of temporary vectors. There is full + support for the in-place version `mul!`, which should be preferred for + higher efficiency in critical algorithms. In addition, it tries to recognize + the properties of combinations of linear maps. In particular, compositions + such as `A'*A` for arbitrary `A` or even `A'*B*C*B'*A` with arbitrary `A` + and `B` and positive definite `C` are recognized as being positive definite + and hermitian. In case a certain property of the resulting `LinearMap` + object is not correctly inferred, the `LinearMap` method can be called to + redefine the properties. ## Methods * `LinearMap` - General purpose method to construct `LinearMap` objects of specific types, as described in the [Types](#types) section below + General purpose method to construct `LinearMap` objects of specific types, + as described in the [Types](#types) section below ``` LinearMap{T}(A::AbstractMatrix[; issymmetric::Bool, ishermitian::Bool, isposdef::Bool]) - LinearMap{T}(A::LinearMap[; issymmetric::Bool, ishermitian::Bool, isposdef::Bool]) + LinearMap{T}(A::LinearMap [; issymmetric::Bool, ishermitian::Bool, isposdef::Bool]) ``` - - Create a `WrappedMap` object that will respond to the methods `issymmetric`, `ishermitian`, `isposdef` with the values provided by the keyword arguments, and to `eltype` with the value `T`. The default values correspond to the result of calling these methods on the argument `A`; in particular `{T}` does not need to be specified and is set as `eltype(A)`. This allows to use an `AbstractMatrix` within the `LinearMap` framework and to redefine the properties of an existing `LinearMap`. + Create a `WrappedMap` object that will respond to the methods `issymmetric`, + `ishermitian`, `isposdef` with the values provided by the keyword arguments, + and to `eltype` with the value `T`. The default values correspond to the + result of calling these methods on the argument `A`; in particular `{T}` + does not need to be specified and is set as `eltype(A)`. This allows to use + an `AbstractMatrix` within the `LinearMap` framework and to redefine the + properties of an existing `LinearMap`. ``` LinearMap{T}(f, [fc = nothing], M::Int, [N::Int = M]; ismutating::Bool, issymmetric::Bool, ishermitian::Bool, isposdef::Bool]) ``` + Create a `FunctionMap` instance that wraps an object describing the action + of the linear map on a vector as a function call. Here, `f` can be a + function or any other object for which either the call + `f(src::AbstractVector) -> dest::AbstractVector` (when `ismutating = false`) + or `f(dest::AbstractVector,src::AbstractVector) -> dest` (when + `ismutating = true`) is supported. The value of `ismutating` can be + specified, by default its value is guessed by looking at the number of + arguments of the first method in the method list of `f`. + + A second function or object `fc` can optionally be provided that implements + the action of the adjoint (transposed) linear map. Here, it is always + assumed that this represents the adjoint/conjugate transpose, though this is + of course equivalent to the normal transpose for real linear maps. + Furthermore, the conjugate transpose also enables the use of + `mul!(y, transpose(A), x)` using some extra conjugation calls on the input + and output vector. If no second function is provided, than + `mul!(y, transpose(A), x)` and `mul!(y, adjoint(A), x)` cannot be used with + this linear map, unless it is symmetric or hermitian. + + `M` is the number of rows (length of the output vectors) and `N` the number + of columns (length of the input vectors). When the latter is not specified, + `N = M`. + + Finally, one can specify the `eltype` of the resulting linear map using the + type parameter `T`. If not specified, a default value of `Float64` is + assumed. Use a complex type `T` if the function represents a complex linear + map. + + The keyword arguments and their default values are: + + * `ismutating`: `false` if the function `f` accepts a single vector + argument corresponding to the input, and `true` if it accepts two vector + arguments where the first will be mutated so as to contain the result. + In both cases, the resulting `A::FunctionMap` will support both the + mutating and non-mutating matrix vector multiplication. Default value is + guessed based on the number of arguments for the first method in the + method list of `f`; it is not possible to use `f` and `fc` where only + one of the two is mutating and the other is not. + * `issymmetric [=false]`: whether the function represents the + multiplication with a symmetric matrix. If `true`, this will + automatically enable `A' * x` and `transpose(A) * x`. + * `ishermitian [=T<:Real && issymmetric]`: whether the function represents + the multiplication with a hermitian matrix. If `true`, this will + automatically enable `A' * x` and `transpose(A) * x`. + * `isposdef [=false]`: whether the function represents the multiplication + with a positive definite matrix. - Create a `FunctionMap` instance that wraps an object describing the action of the linear map on a vector as a function call. Here, `f` can be a function or any other object for which either the call `f(src::AbstractVector) -> dest::AbstractVector` (when `ismutating = false`) or `f(dest::AbstractVector,src::AbstractVector) -> dest` (when `ismutating = true`) is supported. The value of `ismutating` can be specified, by default its value is guessed by looking at the number of arguments of the first method in the method list of `f`. - - A second function or object `fc` can optionally be provided that implements the action of the adjoint (transposed) linear map. Here, it is always assumed that this represents the adjoint/conjugate transpose, though this is of course equivalent to the normal transpose for real linear maps. Furthermore, the conjugate transpose also enables the use of `mul!(y, transpose(A), x)` using some extra conjugation calls on the input and output vector. If no second function is provided, than `mul!(y, transpose(A), x)` and `mul!(y, adjoint(A), x)` cannot be used with this linear map, unless it is symmetric or hermitian. - - `M` is the number of rows (length of the output vectors) and `N` the number of columns (length of the input vectors). When the latter is not specified, `N = M`. - - Finally, one can specify the `eltype` of the resulting linear map using the type parameter `T`. If not specified, a default value of `Float64` is assumed. Use a complex type `T` if the function represents a complex linear map. - - In summary, the keyword arguments and their default values are: - - * `ismutating`: `false` if the function `f` accepts a single vector argument corresponding to the input, and `true` if it accepts two vector arguments where the first will be mutated so as to contain the result. In both cases, the resulting `A::FunctionMap` will support both the mutating and non-mutating matrix vector multiplication. Default value is guessed based on the number of arguments for the first method in the method list of `f`; it is not possible to use `f` and `fc` where only one of the two is mutating and the other is not. - * `issymmetric [=false]`: whether the function represents the multiplication with a symmetric matrix. If `true`, this will automatically enable `A' * x` and `transpose(A) * x`. - * `ishermitian [=T<:Real && issymmetric]`: whether the function represents the multiplication with a hermitian matrix. If `true`, this will automatically enable `A' * x` and `transpose(A) * x`. - * `isposdef [=false]`: whether the function represents the multiplication with a positive definite matrix. + ``` + LinearMap(J::UniformScaling, M::Int) + ``` + Create a `UniformScalingMap` instance that corresponds to a uniform scaling + map of size `M`×`M`. * `Array(A::LinearMap)`, `Matrix(A::LinearMap)`, `convert(Matrix, A::LinearMap)` and `convert(Array, A::LinearMap)` - Create a dense matrix representation of the `LinearMap` object, by multiplying it with the successive basis vectors. This is mostly for testing purposes or if you want to have the explicit matrix representation of a linear map for which you only have a function definition (e.g. to be able to use its `transpose` or `adjoint`). This way, one may conveniently make `A` act on the columns of a matrix `X`, instead of interpreting `A * X` as a composed linear map: `Matrix(A * X)`. For generic code, that is supposed to handle both `A::AbstractMatrix` and `A::LinearMap`, it is recommended to use `convert(Matrix, A*X)`. + Create a dense matrix representation of the `LinearMap` object, by + multiplying it with the successive basis vectors. This is mostly for testing + purposes or if you want to have the explicit matrix representation of a + linear map for which you only have a function definition (e.g. to be able to + use its `transpose` or `adjoint`). This way, one may conveniently make `A` + act on the columns of a matrix `X`, instead of interpreting `A * X` as a + composed linear map: `Matrix(A * X)`. For generic code, that is supposed to + handle both `A::AbstractMatrix` and `A::LinearMap`, it is recommended to use + `convert(Matrix, A*X)`. * `SparseArrays.sparse(A::LinearMap)` and `convert(SparseMatrixCSC, A::LinearMap)` - Create a sparse matrix representation of the `LinearMap` object, by multiplying it with the successive basis vectors. This is mostly for testing purposes or if you want to have the explicit sparse matrix representation of a linear map for which you only have a function definition (e.g. to be able to use its `transpose` or `adjoint`). + Create a sparse matrix representation of the `LinearMap` object, by + multiplying it with the successive basis vectors. This is mostly for testing + purposes or if you want to have the explicit sparse matrix representation of + a linear map for which you only have a function definition (e.g. to be able + to use its `transpose` or `adjoint`). * The following matrix multiplication methods. * `A * x`: applies `A` to `x` and returns the result; - * `mul!(y::AbstractVector, A::LinearMap, x::AbstractVector)`: applies `A` to `x` and stores the result in `y`; - * `mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix)`: applies `A` to each column of `X` and stores the results in the corresponding columns of `Y`; - * `mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number=1, β::Number=0)`: computes `α * A * x + β * y` and stores the result in `y`. Analogously for `X,Y::AbstractMatrix`. + * `mul!(y::AbstractVector, A::LinearMap, x::AbstractVector)`: applies `A` to + `x` and stores the result in `y`; + * `mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix)`: applies `A` to + each column of `X` and stores the results in the corresponding columns of + `Y`; + * `mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number=true, β::Number=false)`: + computes `A * x * α + y * β` and stores the result in `y`. Analogously for `X,Y::AbstractMatrix`. - Applying the adjoint or transpose of `A` (if defined) to `x` works exactly as in the usual matrix case: `transpose(A) * x` and `mul!(y, A', x)`, for instance. + Applying the adjoint or transpose of `A` (if defined) to `x` works exactly + as in the usual matrix case: `transpose(A) * x` and `mul!(y, A', x)`, for instance. * Horizontal, vertical, and hv-concatenation as known for regular matrices, where `UniformScaling`s are automatically promoted to `LinearMap`s and their sizes are inferred, if possible. +* Blockdiagonal concatenation works by extension of `Base.cat(As...; dims=(1,2))` + or `SparseArrays.blockdiag(As...)`. `AbstractArray` arguments are + automatically promoted to `LinearMap` provided there is a `LinearMap` + argument among the first 8 arguments. + +* Kronecker product (via extension of `Base.kron`), Kronecker sum (via + `kronsum`), and respective powers (`kron(A, A, A)`, `A::LinearMap` is + equivalently constructed as `A^⊗(3)` and `kronsum(A, A, A)`, + `A::Union{AbstractMatrix,LinearMap}` is equivalently constructed as + `A^⊕(3)`). + ## Types -None of the types below need to be constructed directly; they arise from performing -operations between `LinearMap` objects or by calling the `LinearMap` constructor -described above. +None of the types below need to be constructed directly; they arise from +performing operations between `LinearMap` objects or by calling the `LinearMap` +constructor described above. * `LinearMap` @@ -150,6 +266,14 @@ described above. Used to add and multiply `LinearMap` objects, don't need to be constructed explicitly. +* `KroneckerMap` and `KroneckerSumMap` + + Types for representing Kronecker products and Kronecker sums, resp., lazily. + +* `BlockMap` and `BlockDiagonalMap` + + Types for representing block (diagonal) maps. + ## Example The `LinearMap` object combines well with methods of the following packages: @@ -161,8 +285,9 @@ The `LinearMap` object combines well with methods of the following packages: ``` using LinearMaps -import Arpack, IterativeSolvers, TSVD +import Arpack, IterativeSolvers, KrylovKit, TSVD +# Example 1, 1-dimensional Laplacian with periodic boundary conditions function leftdiff!(y::AbstractVector, x::AbstractVector) # left difference assuming periodic boundary conditions N = length(x) length(y) == N || throw(DimensionMismatch()) @@ -189,4 +314,22 @@ Arpack.svds(D; nsv=3) Σ, L = IterativeSolvers.svdl(D; nsv=3) TSVD.tsvd(D, 3) + +# Example 2, 1-dimensional Laplacian +A = LinearMap(100; issymmetric=true, ismutating=true) do C, B + C[1] = -2B[1] + B[2] + for i in 2:length(B)-1 + C[i] = B[i-1] - 2B[i] + B[i+1] + end + C[end] = B[end-1] - 2B[end] + return C +end + +Arpack.eigs(-A; nev=3, which=:SR) + +# Example 3, 2-dimensional Laplacian +Δ = kronsum(A, A) + +Arpack.eigs(Δ; nev=3, which=:LR) +KrylovKit.eigsolve(x -> Δ*x, size(Δ, 1), 3, :LR) ``` diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl index 99e2b4fa..373cf5b5 100644 --- a/src/LinearMaps.jl +++ b/src/LinearMaps.jl @@ -48,15 +48,17 @@ Base.length(A::LinearMap) = size(A)[1] * size(A)[2] # check dimension consistency for y = A*x and Y = A*X function check_dim_mul(y::AbstractVector, A::LinearMap, x::AbstractVector) - m, n = size(A) - (m == length(y) && n == length(x)) || throw(DimensionMismatch("mul!")) - return nothing - end - function check_dim_mul(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix) - m, n = size(A) - (m == size(Y, 1) && n == size(X, 1) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!")) - return nothing - end + # @info "checked vector dimensions" # uncomment for testing + m, n = size(A) + (m == length(y) && n == length(x)) || throw(DimensionMismatch("mul!")) + return nothing +end +function check_dim_mul(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix) + # @info "checked matrix dimensions" # uncomment for testing + m, n = size(A) + (m == size(Y, 1) && n == size(X, 1) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!")) + return nothing +end # conversion of AbstractMatrix to LinearMap convert_to_lmaps_(A::AbstractMatrix) = LinearMap(A) @@ -66,9 +68,12 @@ convert_to_lmaps(A) = (convert_to_lmaps_(A),) @inline convert_to_lmaps(A, B, Cs...) = (convert_to_lmaps_(A), convert_to_lmaps_(B), convert_to_lmaps(Cs...)...) -Base.:(*)(A::LinearMap, x::AbstractVector) = mul!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x) +function Base.:(*)(A::LinearMap, x::AbstractVector) + size(A, 2) == length(x) || throw(DimensionMismatch("mul!")) + return @inbounds mul!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x) +end function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, α::Number=true, β::Number=false) - length(y) == size(A, 1) || throw(DimensionMismatch("mul!")) + @boundscheck check_dim_mul(y, A, x) if isone(α) iszero(β) && (A_mul_B!(y, A, x); return y) isone(β) && (y .+= A * x; return y) @@ -92,8 +97,8 @@ function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector, end end # the following is of interest in, e.g., subspace-iteration methods -function LinearAlgebra.mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix, α::Number=true, β::Number=false) - (size(Y, 1) == size(A, 1) && size(X, 1) == size(A, 2) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!")) +Base.@propagate_inbounds function LinearAlgebra.mul!(Y::AbstractMatrix, A::LinearMap, X::AbstractMatrix, α::Number=true, β::Number=false) + @boundscheck check_dim_mul(Y, A, X) @inbounds @views for i = 1:size(X, 2) mul!(Y[:, i], A, X[:, i], α, β) end diff --git a/src/blockmap.jl b/src/blockmap.jl index d7acf1d9..f5caad7a 100644 --- a/src/blockmap.jl +++ b/src/blockmap.jl @@ -1,14 +1,14 @@ -struct BlockMap{T,As<:Tuple{Vararg{LinearMap}},Rs<:Tuple{Vararg{Int}}} <: LinearMap{T} +struct BlockMap{T,As<:Tuple{Vararg{LinearMap}},Rs<:Tuple{Vararg{Int}},Rranges<:Tuple{Vararg{UnitRange{Int}}},Cranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T} maps::As rows::Rs - rowranges::Vector{UnitRange{Int}} - colranges::Vector{UnitRange{Int}} + rowranges::Rranges + colranges::Cranges function BlockMap{T,R,S}(maps::R, rows::S) where {T, R<:Tuple{Vararg{LinearMap}}, S<:Tuple{Vararg{Int}}} for A in maps promote_type(T, eltype(A)) == T || throw(InexactError()) end rowranges, colranges = rowcolranges(maps, rows) - return new{T,R,S}(maps, rows, rowranges, colranges) + return new{T,R,S,typeof(rowranges),typeof(colranges)}(maps, rows, rowranges, colranges) end end @@ -28,28 +28,28 @@ Determines the range of rows for each block row and the range of columns for eac map in `maps`, according to its position in a virtual matrix representation of the block linear map obtained from `hvcat(rows, maps...)`. """ -function rowcolranges(maps, rows)::Tuple{Vector{UnitRange{Int}},Vector{UnitRange{Int}}} - rowranges = Vector{UnitRange{Int}}(undef, length(rows)) - colranges = Vector{UnitRange{Int}}(undef, length(maps)) +function rowcolranges(maps, rows) + rowranges = () + colranges = () mapind = 0 rowstart = 1 - for rowind in 1:length(rows) - xinds = vcat(1, map(a -> size(a, 2), maps[mapind+1:mapind+rows[rowind]])...) + for row in rows + xinds = vcat(1, map(a -> size(a, 2), maps[mapind+1:mapind+row])...) cumsum!(xinds, xinds) mapind += 1 rowend = rowstart + size(maps[mapind], 1) - 1 - rowranges[rowind] = rowstart:rowend - colranges[mapind] = xinds[1]:xinds[2]-1 - for colind in 2:rows[rowind] + rowranges = (rowranges..., rowstart:rowend) + colranges = (colranges..., xinds[1]:xinds[2]-1) + for colind in 2:row mapind +=1 - colranges[mapind] = xinds[colind]:xinds[colind+1]-1 + colranges = (colranges..., xinds[colind]:xinds[colind+1]-1) end rowstart = rowend + 1 end - return rowranges, colranges + return rowranges::NTuple{length(rows), UnitRange{Int}}, colranges::NTuple{length(maps), UnitRange{Int}} end -Base.size(A::BlockMap) = (last(A.rowranges[end]), last(A.colranges[end])) +Base.size(A::BlockMap) = (last(last(A.rowranges)), last(last(A.colranges))) ############ # concatenation @@ -299,75 +299,82 @@ LinearAlgebra.transpose(A::BlockMap) = TransposeMap(A) LinearAlgebra.adjoint(A::BlockMap) = AdjointMap(A) ############ -# multiplication with vectors +# multiplication helper functions ############ -function A_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) - require_one_based_indexing(y, x) - m, n = size(A) - @boundscheck (m == length(y) && n == length(x)) || throw(DimensionMismatch("A_mul_B!")) +@inline function _blockmul!(y, A::BlockMap, x, α, β) maps, rows, yinds, xinds = A.maps, A.rows, A.rowranges, A.colranges mapind = 0 - @views @inbounds for rowind in 1:length(rows) - yrow = y[yinds[rowind]] + @views @inbounds for (row, yi) in zip(rows, yinds) + yrow = selectdim(y, 1, yi) mapind += 1 - A_mul_B!(yrow, maps[mapind], x[xinds[mapind]]) - for colind in 2:rows[rowind] + mul!(yrow, maps[mapind], selectdim(x, 1, xinds[mapind]), α, β) + for _ in 2:row mapind +=1 - mul!(yrow, maps[mapind], x[xinds[mapind]], true, true) + mul!(yrow, maps[mapind], selectdim(x, 1, xinds[mapind]), α, true) end end return y end -function At_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) - require_one_based_indexing(y, x) - m, n = size(A) - @boundscheck (n == length(y) && m == length(x)) || throw(DimensionMismatch("At_mul_B!")) +@inline function _transblockmul!(y, A::BlockMap, x, α, β, transform) maps, rows, xinds, yinds = A.maps, A.rows, A.rowranges, A.colranges - mapind = 0 - # first block row (rowind = 1) of A, meaning first block column of A', fill all of y @views @inbounds begin - xcol = x[xinds[1]] - for colind in 1:rows[1] - mapind +=1 - A_mul_B!(y[yinds[mapind]], transpose(maps[mapind]), xcol) + # first block row (rowind = 1) of A, meaning first block column of A', fill all of y + xcol = selectdim(x, 1, first(xinds)) + for rowind in 1:first(rows) + mul!(selectdim(y, 1, yinds[rowind]), transform(maps[rowind]), xcol, α, β) end - # subsequent block rows of A, add results to corresponding parts of y - for rowind in 2:length(rows) - xcol = x[xinds[rowind]] - for colind in 1:rows[rowind] + mapind = first(rows) + # subsequent block rows of A (block columns of A'), + # add results to corresponding parts of y + # TODO: think about multithreading + for (row, xi) in zip(Base.tail(rows), Base.tail(xinds)) + xcol = selectdim(x, 1, xi) + for _ in 1:row mapind +=1 - mul!(y[yinds[mapind]], transpose(maps[mapind]), xcol, true, true) + mul!(selectdim(y, 1, yinds[mapind]), transform(maps[mapind]), xcol, α, true) end end end return y end -function Ac_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) - require_one_based_indexing(y, x) - m, n = size(A) - @boundscheck (n == length(y) && m == length(x)) || throw(DimensionMismatch("At_mul_B!")) - maps, rows, xinds, yinds = A.maps, A.rows, A.rowranges, A.colranges - mapind = 0 - # first block row (rowind = 1) of A, fill all of y - @views @inbounds begin - xcol = x[xinds[1]] - for colind in 1:rows[1] - mapind +=1 - A_mul_B!(y[yinds[mapind]], adjoint(maps[mapind]), xcol) - end - # subsequent block rows of A, add results to corresponding parts of y - for rowind in 2:length(rows) - xcol = x[xinds[rowind]] - for colind in 1:rows[rowind] - mapind +=1 - mul!(y[yinds[mapind]], adjoint(maps[mapind]), xcol, true, true) - end +############ +# multiplication with vectors & matrices +############ + +Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) = + mul!(y, A, x) + +Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::TransposeMap{<:Any,<:BlockMap}, x::AbstractVector) = + mul!(y, A, x) + +Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) = + mul!(y, transpose(A), x) + +Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::AdjointMap{<:Any,<:BlockMap}, x::AbstractVector) = + mul!(y, A, x) + +Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockMap, x::AbstractVector) = + mul!(y, adjoint(A), x) + +for Atype in (AbstractVector, AbstractMatrix) + @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, A::BlockMap, x::$Atype, + α::Number=true, β::Number=false) + require_one_based_indexing(y, x) + @boundscheck check_dim_mul(y, A, x) + return _blockmul!(y, A, x, α, β) + end + + for (maptype, transform) in ((:(TransposeMap{<:Any,<:BlockMap}), :transpose), (:(AdjointMap{<:Any,<:BlockMap}), :adjoint)) + @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, wrapA::$maptype, x::$Atype, + α::Number=true, β::Number=false) + require_one_based_indexing(y, x) + @boundscheck check_dim_mul(y, wrapA, x) + return _transblockmul!(y, wrapA.lmap, x, α, β, $transform) end end - return y end ############ @@ -388,3 +395,91 @@ end # show(io, T) # print(io, '}') # end + +############ +# BlockDiagonalMap +############ + +struct BlockDiagonalMap{T,As<:Tuple{Vararg{LinearMap}},Ranges<:Tuple{Vararg{UnitRange{Int}}}} <: LinearMap{T} + maps::As + rowranges::Ranges + colranges::Ranges + function BlockDiagonalMap{T,As}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} + for A in maps + promote_type(T, eltype(A)) == T || throw(InexactError()) + end + # row ranges + inds = vcat(1, size.(maps, 1)...) + cumsum!(inds, inds) + rowranges = ntuple(i -> inds[i]:inds[i+1]-1, Val(length(maps))) + # column ranges + inds[2:end] .= size.(maps, 2) + cumsum!(inds, inds) + colranges = ntuple(i -> inds[i]:inds[i+1]-1, Val(length(maps))) + return new{T,As,typeof(rowranges)}(maps, rowranges, colranges) + end +end + +BlockDiagonalMap{T}(maps::As) where {T,As<:Tuple{Vararg{LinearMap}}} = + BlockDiagonalMap{T,As}(maps) +BlockDiagonalMap(maps::LinearMap...) = + BlockDiagonalMap{promote_type(map(eltype, maps)...)}(maps) + +for k in 1:8 # is 8 sufficient? + Is = ntuple(n->:($(Symbol(:A,n))::AbstractMatrix), Val(k-1)) + # yields (:A1, :A2, :A3, ..., :A(k-1)) + L = :($(Symbol(:A,k))::LinearMap) + # yields :Ak + mapargs = ntuple(n -> :(LinearMap($(Symbol(:A,n)))), Val(k-1)) + # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1))) + + @eval begin + SparseArrays.blockdiag($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...) = + BlockDiagonalMap($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...) + function Base.cat($(Is...), $L, As::Union{LinearMap,AbstractMatrix}...; dims::Dims{2}) + if dims == (1,2) + return BlockDiagonalMap($(mapargs...), $(Symbol(:A,k)), convert_to_lmaps(As...)...) + else + throw(ArgumentError("dims keyword in cat of LinearMaps must be (1,2)")) + end + end + end +end + +Base.size(A::BlockDiagonalMap) = (last(A.rowranges[end]), last(A.colranges[end])) + +LinearAlgebra.issymmetric(A::BlockDiagonalMap) = all(issymmetric, A.maps) +LinearAlgebra.ishermitian(A::BlockDiagonalMap{<:Real}) = all(issymmetric, A.maps) +LinearAlgebra.ishermitian(A::BlockDiagonalMap) = all(ishermitian, A.maps) + +LinearAlgebra.adjoint(A::BlockDiagonalMap{T}) where {T} = BlockDiagonalMap{T}(map(adjoint, A.maps)) +LinearAlgebra.transpose(A::BlockDiagonalMap{T}) where {T} = BlockDiagonalMap{T}(map(transpose, A.maps)) + +Base.:(==)(A::BlockDiagonalMap, B::BlockDiagonalMap) = (eltype(A) == eltype(B) && A.maps == B.maps) + +Base.@propagate_inbounds A_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) = + mul!(y, A, x, true, false) + +Base.@propagate_inbounds At_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) = + mul!(y, transpose(A), x, true, false) + +Base.@propagate_inbounds Ac_mul_B!(y::AbstractVector, A::BlockDiagonalMap, x::AbstractVector) = + mul!(y, adjoint(A), x, true, false) + +for Atype in (AbstractVector, AbstractMatrix) + @eval Base.@propagate_inbounds function LinearAlgebra.mul!(y::$Atype, A::BlockDiagonalMap, x::$Atype, + α::Number=true, β::Number=false) + require_one_based_indexing(y, x) + @boundscheck check_dim_mul(y, A, x) + return _blockscaling!(y, A, x, α, β) + end +end + +@inline function _blockscaling!(y, A::BlockDiagonalMap, x, α, β) + maps, yinds, xinds = A.maps, A.rowranges, A.colranges + # TODO: think about multi-threading here + @views @inbounds for i in eachindex(yinds, maps, xinds) + mul!(selectdim(y, 1, yinds[i]), maps[i], selectdim(x, 1, xinds[i]), α, β) + end + return y +end diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl index 212c0de7..49ee410d 100644 --- a/src/wrappedmap.jl +++ b/src/wrappedmap.jl @@ -30,28 +30,13 @@ Base.:(==)(A::MatrixMap, B::MatrixMap) = (eltype(A)==eltype(B) && A.lmap==B.lmap && A._issymmetric==B._issymmetric && A._ishermitian==B._ishermitian && A._isposdef==B._isposdef) -if VERSION ≥ v"1.3.0-alpha.115" - -LinearAlgebra.mul!(y::AbstractVector, A::WrappedMap, x::AbstractVector, α::Number=true, β::Number=false) = - mul!(y, A.lmap, x, α, β) - -LinearAlgebra.mul!(Y::AbstractMatrix, A::MatrixMap, X::AbstractMatrix, α::Number=true, β::Number=false) = - mul!(Y, A.lmap, X, α, β) - -else - -LinearAlgebra.mul!(Y::AbstractMatrix, A::MatrixMap, X::AbstractMatrix) = - mul!(Y, A.lmap, X) - -end # VERSION - # properties Base.size(A::WrappedMap) = size(A.lmap) LinearAlgebra.issymmetric(A::WrappedMap) = A._issymmetric LinearAlgebra.ishermitian(A::WrappedMap) = A._ishermitian LinearAlgebra.isposdef(A::WrappedMap) = A._isposdef -# multiplication with vector +# multiplication with vectors & matrices A_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = A_mul_B!(y, A.lmap, x) Base.:(*)(A::WrappedMap, x::AbstractVector) = *(A.lmap, x) @@ -61,6 +46,23 @@ At_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = Ac_mul_B!(y::AbstractVector, A::WrappedMap, x::AbstractVector) = ishermitian(A) ? A_mul_B!(y, A.lmap, x) : Ac_mul_B!(y, A.lmap, x) +if VERSION ≥ v"1.3.0-alpha.115" + for Atype in (AbstractVector, AbstractMatrix) + @eval Base.@propagate_inbounds LinearAlgebra.mul!(y::$Atype, A::WrappedMap, x::$Atype, + α::Number=true, β::Number=false) = + mul!(y, A.lmap, x, α, β) + end +else +# This is somewhat suboptimal, because the absence of 5-arg mul! for MatrixMaps +# doesn't allow to define a 5-arg mul! for WrappedMaps which do have a 5-arg mul! +# I'd assume, however, that 5-arg mul! becomes standard in Julia v≥1.3 anyway +# the idea is to let the fallback handle 5-arg calls + for Atype in (AbstractVector, AbstractMatrix) + @eval Base.@propagate_inbounds LinearAlgebra.mul!(Y::$Atype, A::WrappedMap, X::$Atype) = + mul!(Y, A.lmap, X) + end +end # VERSION + # combine LinearMap and Matrix objects: linear combinations and map composition Base.:(+)(A₁::LinearMap, A₂::AbstractMatrix) = +(A₁, WrappedMap(A₂)) Base.:(+)(A₁::AbstractMatrix, A₂::LinearMap) = +(WrappedMap(A₁), A₂) diff --git a/test/blockmap.jl b/test/blockmap.jl index 734a9edc..e44119fb 100644 --- a/test/blockmap.jl +++ b/test/blockmap.jl @@ -1,4 +1,4 @@ -using Test, LinearMaps, LinearAlgebra +using Test, LinearMaps, LinearAlgebra, SparseArrays @testset "block maps" begin @testset "hcat" begin @@ -79,6 +79,15 @@ using Test, LinearMaps, LinearAlgebra @test size(L) == (30, 30) @test Matrix(L) ≈ A @test L * x ≈ A * x + y = randn(elty, size(L, 1)) + for α in (0, 1, rand(elty)), β in (0, 1, rand(elty)) + @test mul!(copy(y), L, x, α, β) ≈ y*β .+ A*x*α + end + X = rand(elty, 30, 10) + Y = randn(elty, size(L, 1), 10) + for α in (0, 1, rand(elty)), β in (0, 1, rand(elty)) + @test mul!(copy(Y), L, X, α, β) ≈ Y*β .+ A*X*α + end A = rand(elty, 10,10); LA = LinearMap(A) B = rand(elty, 20,30); LB = LinearMap(B) @test [LA LA LA; LB] isa LinearMaps.BlockMap{elty} @@ -140,7 +149,50 @@ using Test, LinearMaps, LinearAlgebra Lt = @inferred transform(L) @test Lt isa LinearMaps.LinearMap{elty} @test Lt * x ≈ transform(A) * x - @test Matrix(Lt) ≈ Matrix(transform(A)) + @test Matrix(Lt) ≈ Matrix(transform(LinearMap(L))) ≈ Matrix(transform(A)) + @test Matrix(transform(LinearMap(L))+transform(LinearMap(L))) ≈ 2Matrix(transform(A)) + X = rand(elty, size(L, 1), 10) + Y = randn(elty, size(L, 2), 10) + for α in (0, 1, rand(elty)), β in (0, 1, rand(elty)) + @test mul!(copy(Y), Lt, X, α, β) ≈ Y*β .+ transform(A)*X*α + end + end + end + + @testset "block diagonal maps" begin + for elty in (Float64, ComplexF64) + m = 5; n = 6 + M1 = 10*(1:m) .+ (1:(n+1))'; L1 = LinearMap(M1) + M2 = randn(elty, m, n+2); L2 = LinearMap(M2) + M3 = randn(elty, m, n+3); L3 = LinearMap(M3) + + # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse: + Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...)) + x = randn(elty, size(Md, 2)) + Bd = @inferred blockdiag(L1, L2, L3, L2, L1) + @test Matrix(@inferred blockdiag(L1)) == M1 + @test Matrix(@inferred blockdiag(L1, L2)) == blockdiag(sparse.((M1, M2))...) + Bd2 = @inferred cat(L1, L2, L3, L2, L1; dims=(1,2)) + @test_throws ArgumentError cat(L1, L2, L3, L2, L1; dims=(2,2)) + @test Bd == Bd2 + @test Bd == blockdiag(L1, M2, M3, M2, M1) + @test size(Bd) == (25, 39) + @test !issymmetric(Bd) + @test !ishermitian(Bd) + @test @inferred Bd * x ≈ Md * x + for transform in (identity, adjoint, transpose) + @test Matrix(@inferred transform(Bd)) == transform(Md) + @test Matrix(@inferred transform(LinearMap(Bd))) == transform(Md) + end + y = randn(elty, size(Md, 1)) + for α in (0, 1, rand(elty)), β in (0, 1, rand(elty)) + @test mul!(copy(y), Bd, x, α, β) ≈ y*β .+ Md*x*α + end + X = randn(elty, size(Md, 2), 10) + Y = randn(elty, size(Md, 1), 10) + for α in (0, 1, rand(elty)), β in (0, 1, rand(elty)) + @test mul!(copy(Y), Bd, X, α, β) ≈ Y*β .+ Md*X*α + end end end end