Skip to content

Commit d00a916

Browse files
dkarraschJutho
authored andcommitted
Julia 0.7 upgrade (#23)
* Julia 0.7 upgrade * adjusted README * simplify Transpose/AdjointMap constructors * Update travis.yml for Julia 0.7 compliance * expanded qualifiers, added a few tests * rm Ax_mul_B, typos, code style, more tests * removed all `Ax_mul_B` (w/o !) since they are neither exported nor called internally * added a few more tests to increase coverage * improved code consistency: brackets, spaces, etc. * constructors, code style * removed (::Type{LinearMap{T}} in constructors * polished code
1 parent 1dd072f commit d00a916

File tree

12 files changed

+267
-336
lines changed

12 files changed

+267
-336
lines changed

.DS_Store

6 KB
Binary file not shown.

.travis.yml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ os:
44
- linux
55
- osx
66
julia:
7-
- 0.6
7+
- 0.7
88
- nightly
99

1010
matrix:

README.md

Lines changed: 12 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
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.
99

1010
## What's new.
11-
* Fully julia v0.6 compatible; dropped compatibility for previous versions of Julia.
11+
* Fully julia v0.7 compatible; dropped compatibility for previous versions of Julia.
1212

1313
* `LinearMap` is now the name of the abstract type on top (used to be `AbstractLinearMap` because you could not have a function/constructor with the same name as an abstract type in older julia versions)
1414

@@ -24,11 +24,11 @@ Several iterative linear algebra methods such as linear solvers or eigensolvers
2424

2525
The LinearMaps package provides the following functionality:
2626

27-
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 `ctranspose` and multiplication with a vector using both `*` or the in-place version `A_mul_B!`. Depending on the subtype, also `At_mul_B`, `At_mul_B!`, `Ac_mul_B` and `Ac_mul_B!` are supported. Linear algebra functions that uses duck-typing for its 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!`.
27+
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 uses duck-typing for its 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!`.
2828

2929
2. A single method `LinearMap` function that acts as a general purpose constructor (though it 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.
3030

31-
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 `A_mul_B!`, 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.
31+
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.
3232

3333
## Methods
3434

@@ -49,7 +49,7 @@ The LinearMaps package provides the following functionality:
4949
5050
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 spefified, by default its value is guessed by looking at the number of arguments of the first method in the method list of `f`.
5151
52-
A second function or object can optionally be provided that implements the action of the adjoint (transposed) linear map. Here, it is always assumed that this represents the 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 `At_mul_B(!)` using some extra conjugation calls on the input and output vector. If no second function is provided, than `At_mul_B(!)` and `Ac_mul_B(!)` cannot be used with this linear map, unless it is symmetric or hermitian.
52+
A second function or object 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.
5353
5454
`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`.
5555
@@ -62,9 +62,13 @@ The LinearMaps package provides the following functionality:
6262
* `ishermitian [=T<:Real && issymmetric]`: whether the function represents the multiplication with a hermitian matrix. If `true`, this will automatically enable `A'*x` and `A.'*x`.
6363
* `isposdef [=false]`: whether the function represents the multiplication with a positive definite matrix.
6464
65-
* `Base.full(linearmap)`
65+
* `Base.Array(linearmap)`
6666
67-
Creates a full 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 `(c)transpose`).
67+
Creates 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`).
68+
69+
* `SparseArrays.sparse(linearmap)`
70+
71+
Creates 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`).
6872
6973
* All matrix multiplication methods and the corresponding mutating versions.
7074
@@ -88,13 +92,13 @@ None of the types below need to be constructed directly; they arise from perform
8892
8993
Type for representing the identity map of a certain size `M=N`, obtained simply as `IdentityMap{T}(M)`, `IdentityMap(T,M)=IdentityMap(T,M,N)=IdentityMap(T,(M,N))` or even `IdentityMap(M)=IdentityMap(M,N)=IdentityMap((M,N))`. If `T` is not specified, `Bool` is assumed, since operations between `Bool` and any other `Number` will always be converted to the type of the other `Number`. If `M!=N`, an error is returned. An `IdentityMap` of the correct size and element type will automatically be created if `LinearMap` objects are combined with `I`, Julia's built in identity (`UniformScaling`).
9094
91-
* `LinearCombination`, `CompositeMap`, `TransposeMap` and `CTransposeMap`
95+
* `LinearCombination`, `CompositeMap`, `TransposeMap` and `AdjointMap`
9296
9397
Used to add and multiply `LinearMap` objects, don't need to be constructed explicitly.
9498
9599
## Examples
96100
97-
The `LinearMap` object combines well with the iterative eigensolver `eigs`, which is the Julia wrapper for Arpack.
101+
The `LinearMap` object combines well with the iterative eigensolver `eigs` from [Arpack.jl](https://github.com/JuliaLinearAlgebra/Arpack.jl) and with iterative solvers from [IterativeSolvers.jl](https://github.com/JuliaMath/IterativeSolvers.jl).
98102
99103
```
100104
using LinearMaps

REQUIRE

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
julia 0.6
1+
julia 0.7-

src/LinearMaps.jl

Lines changed: 37 additions & 85 deletions
Original file line numberDiff line numberDiff line change
@@ -1,116 +1,76 @@
11
__precompile__(true)
22
module LinearMaps
33

4-
export LinearMap, AbstractLinearMap
4+
export LinearMap
55

6-
import Base: +, -, *, \, /, ==, transpose
7-
8-
if VERSION >= v"0.7.0-DEV.1415"
9-
const adjoint = Base.adjoint
10-
else
11-
const adjoint = Base.ctranspose
12-
end
6+
using LinearAlgebra
7+
import SparseArrays
138

9+
import Base: +, -, *, \, /, ==
1410

1511
abstract type LinearMap{T} end
1612

17-
const AbstractLinearMap = LinearMap # will be deprecated
18-
1913
Base.eltype(::LinearMap{T}) where {T} = T
20-
Base.eltype(::Type{L}) where {T,L<:LinearMap{T}} = T
2114

2215
Base.isreal(A::LinearMap) = eltype(A) <: Real
23-
Base.issymmetric(::LinearMap) = false # default assumptions
24-
Base.ishermitian(A::LinearMap{<:Real}) = issymmetric(A)
25-
Base.ishermitian(::LinearMap) = false # default assumptions
26-
Base.isposdef(::LinearMap) = false # default assumptions
16+
LinearAlgebra.issymmetric(::LinearMap) = false # default assumptions
17+
LinearAlgebra.ishermitian(A::LinearMap{<:Real}) = issymmetric(A)
18+
LinearAlgebra.ishermitian(::LinearMap) = false # default assumptions
19+
LinearAlgebra.isposdef(::LinearMap) = false # default assumptions
2720

2821
Base.ndims(::LinearMap) = 2
2922
Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions"))
3023
Base.length(A::LinearMap) = size(A)[1] * size(A)[2]
3124

32-
# any LinearMap subtype will have to overwrite at least one of the two following methods to avoid running in circles
33-
*(A::LinearMap, x::AbstractVector) = Base.A_mul_B!(similar(x, promote_type(eltype(A),eltype(x)), size(A,1)), A, x)
34-
Base.A_mul_B!(y::AbstractVector, A::LinearMap, x::AbstractVector) = begin
35-
length(y) == size(A,1) || throw(DimensionMismatch("A_mul_B!"))
36-
copy!(y, A*x)
25+
Base.:(*)(A::LinearMap, x::AbstractVector) = mul!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x)
26+
function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap, x::AbstractVector)
27+
length(y) == size(A, 1) || throw(DimensionMismatch("mul!"))
28+
A_mul_B!(y, A, x)
3729
end
3830

39-
# the following for multiplying with transpose and adjoint map are optional:
40-
# subtypes can overwrite nonmutating methods, implement mutating methods or do nothing
41-
function Base.At_mul_B(A::LinearMap, x::AbstractVector)
42-
l = methods(Base.At_mul_B!,Tuple{AbstractVector, typeof(A), AbstractVector})
43-
if length(l) > 0 && first(l.ms).sig.parameters[3] != LinearMap
44-
Base.At_mul_B!(similar(x, promote_type(eltype(A), eltype(x)), size(A,2)), A, x)
45-
else
46-
throw(MethodError(Base.At_mul_B, (A, x)))
47-
end
48-
end
49-
function Base.At_mul_B!(y::AbstractVector, A::LinearMap, x::AbstractVector)
50-
length(y) == size(A, 2) || throw(DimensionMismatch("At_mul_B!"))
51-
l = methods(Base.At_mul_B,Tuple{typeof(A), AbstractVector})
52-
if length(l) > 0 && first(l.ms).sig.parameters[2] != LinearMap
53-
copy!(y, Base.At_mul_B(A, x))
54-
else
55-
throw(MethodError(Base.At_mul_B!, (y, A, x)))
56-
end
57-
end
58-
function Base.Ac_mul_B(A::LinearMap,x::AbstractVector)
59-
l = methods(Base.Ac_mul_B!,Tuple{AbstractVector, typeof(A), AbstractVector})
60-
if length(l) > 0 && first(l.ms).sig.parameters[3] != LinearMap
61-
Base.Ac_mul_B!(similar(x, promote_type(eltype(A), eltype(x)), size(A,2)), A, x)
62-
else
63-
throw(MethodError(Base.Ac_mul_B, (A, x)))
64-
end
65-
end
66-
function
67-
Base.Ac_mul_B!(y::AbstractVector, A::LinearMap, x::AbstractVector)
68-
length(y)==size(A,2) || throw(DimensionMismatch("At_mul_B!"))
69-
l = methods(Base.Ac_mul_B,Tuple{typeof(A), AbstractVector})
70-
if length(l) > 0 && first(l.ms).sig.parameters[2] != LinearMap
71-
copy!(y, Base.Ac_mul_B(A, x))
72-
else
73-
throw(MethodError(Base.Ac_mul_B!, (y, A, x)))
74-
end
75-
end
31+
A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
32+
At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x)
33+
Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x)
7634

77-
# full: create matrix representation of LinearMap
78-
function Base.full(A::LinearMap)
35+
# Matrix: create matrix representation of LinearMap
36+
function Base.Matrix(A::LinearMap)
7937
M, N = size(A)
8038
T = eltype(A)
81-
mat = zeros(T, (M, N))
82-
v = zeros(T, N)
39+
mat = Matrix{T}(undef, (M, N))
40+
v = fill(zero(T), N)
8341
for i = 1:N
8442
v[i] = one(T)
85-
A_mul_B!(view(mat,:,i), A, v)
43+
mul!(view(mat, :, i), A, v)
8644
v[i] = zero(T)
8745
end
8846
return mat
8947
end
9048

91-
# sparse: create sparse matrix representtion of LinearMap
92-
function Base.sparse(A::LinearMap)
49+
Base.Array(A::LinearMap) = Base.Matrix(A)
50+
51+
# sparse: create sparse matrix representation of LinearMap
52+
function SparseArrays.sparse(A::LinearMap)
9353
M, N = size(A)
9454
T = eltype(A)
9555
rowind = Int[]
9656
nzval = T[]
97-
colptr = Vector{Int}(N+1)
98-
v = zeros(T, N)
57+
colptr = Vector{Int}(undef, N+1)
58+
v = fill(zero(T), N)
9959

10060
for i = 1:N
10161
v[i] = one(T)
102-
Lv = A*v
103-
js = find(Lv)
104-
colptr[i] = length(nzval)+1
62+
Lv = A * v
63+
js = findall(!iszero, Lv)
64+
colptr[i] = length(nzval) + 1
10565
if length(js) > 0
10666
append!(rowind, js)
10767
append!(nzval, Lv[js])
10868
end
10969
v[i] = zero(T)
11070
end
111-
colptr[N+1] = length(nzval)+1
112-
113-
return SparseMatrixCSC(M, N, colptr, rowind, nzval)
71+
colptr[N+1] = length(nzval) + 1
72+
73+
return SparseArrays.SparseMatrixCSC(M, N, colptr, rowind, nzval)
11474
end
11575

11676
include("transpose.jl") # transposing linear maps
@@ -127,7 +87,7 @@ include("functionmap.jl") # using a function as linear map
12787
Construct a linear map object, either from an existing `LinearMap` or `AbstractMatrix` `A`,
12888
with the purpose of redefining its properties via the keyword arguments `kwargs`, or
12989
from a function or callable object `f`. In the latter case, one also needs to specialize
130-
the size of the equivalent matrix representation `(M,N)`, i.e. for functions `f` acting
90+
the size of the equivalent matrix representation `(M, N)`, i.e. for functions `f` acting
13191
on length `N` vectors and producing length `M` vectors (with default value `N=M`). Preferably,
13292
also the `eltype` `T` of the corresponding matrix representation needs to be specified, i.e.
13393
whether the action of `f` on a vector will be similar to e.g. multiplying by numbers of type `T`.
@@ -139,7 +99,7 @@ The keyword arguments and their default values for functions `f` are
13999
* ishermitian::Bool = issymmetric & T<:Real : whether `A` or `f` acts as a Hermitian matrix
140100
* isposdef::Bool = false : whether `A` or `f` acts as a positive definite matrix.
141101
For existing linear maps or matrices `A`, the default values will be taken by calling
142-
`issymmetric`, `ishermitian` and `isposdef` on the exising object `A`.
102+
`issymmetric`, `ishermitian` and `isposdef` on the existing object `A`.
143103
144104
For functions `f`, there is one more keyword arguments
145105
* ismutating::Bool : flags whether the function acts as a mutating matrix multiplication
@@ -148,21 +108,13 @@ For functions `f`, there is one more keyword arguments
148108
The default value is guessed by looking at the number of arguments of the first occurence
149109
of `f` in the method table.
150110
"""
151-
LinearMap(A::Union{AbstractMatrix,LinearMap}; kwargs...) = WrappedMap(A; kwargs...)
111+
LinearMap(A::Union{AbstractMatrix, LinearMap}; kwargs...) = WrappedMap(A; kwargs...)
152112
LinearMap(f, M::Int; kwargs...) = LinearMap{Float64}(f, M; kwargs...)
153113
LinearMap(f, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, M, N; kwargs...)
154114
LinearMap(f, fc, M::Int; kwargs...) = LinearMap{Float64}(f, fc, M; kwargs...)
155115
LinearMap(f, fc, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, fc, M, N; kwargs...)
156116

157-
(::Type{LinearMap{T}})(A::Union{AbstractMatrix,LinearMap}; kwargs...) where {T} = WrappedMap{T}(A; kwargs...)
158-
(::Type{LinearMap{T}})(f, args...; kwargs...) where {T} = FunctionMap{T}(f, args...; kwargs...)
159-
160-
@deprecate LinearMap(f, T::Type, args...; kwargs...) LinearMap{T}(f, args...; kwargs...)
161-
@deprecate LinearMap(f, fc, T::Type, args...; kwargs...) LinearMap{T}(f, fc, args...; kwargs...)
162-
163-
@deprecate LinearMap(f, M::Int, T::Type; kwargs...) LinearMap{T}(f, M; kwargs...)
164-
@deprecate LinearMap(f, M::Int, N::Int, T::Type; kwargs...) LinearMap{T}(f, M, N; kwargs...)
165-
@deprecate LinearMap(f, fc, M::Int, T::Type; kwargs...) LinearMap{T}(f, fc, M; kwargs...)
166-
@deprecate LinearMap(f, fc, M::Int, N::Int, T::Type; kwargs...) LinearMap{T}(f, fc, M, N; kwargs...)
117+
LinearMap{T}(A::Union{AbstractMatrix, LinearMap}; kwargs...) where {T} = WrappedMap{T}(A; kwargs...)
118+
LinearMap{T}(f, args...; kwargs...) where {T} = FunctionMap{T}(f, args...; kwargs...)
167119

168120
end # module

0 commit comments

Comments
 (0)