Julia 1.0 tests, sparse convert, mutating matmul in sparse (#28)

* minor Julia 1.0 edits

* travis test also against Julia v1.0
* removed obsolete precompilation statement
* added sparse matrix convert
* sparse matrix construction uses mutating multiplication

* expanded README and tests

* consistent code style

Author: dkarrasch

.travis.yml

@@ -5,6 +5,7 @@ os:
   - osx
 julia:
   - 0.7
+  - 1.0
   - nightly
 
 matrix:

README.md

@@ -9,6 +9,10 @@
 
 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.2.0
+*   Fully Julia v0.7/v1.0 compatible.
+*   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 `β`.

src/LinearMaps.jl

@@ -1,10 +1,9 @@
-__precompile__(true)
 module LinearMaps
 
 export LinearMap
 
 using LinearAlgebra
-import SparseArrays
+using SparseArrays
 
 abstract type LinearMap{T} end
 
@@ -72,33 +71,34 @@ function Base.Matrix(A::LinearMap)
     return mat
 end
 
-Base.Array(A::LinearMap) = Base.Matrix(A)
-Base.convert(::Type{Matrix}, A:: LinearMap) = Matrix(A)
-Base.convert(::Type{Array}, A:: LinearMap) = Matrix(A)
+Base.Array(A::LinearMap) = Matrix(A)
+Base.convert(::Type{Matrix}, A::LinearMap) = Matrix(A)
+Base.convert(::Type{Array}, A::LinearMap) = Matrix(A)
+Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A)
 
 # sparse: create sparse matrix representation of LinearMap
-function SparseArrays.sparse(A::LinearMap)
+function SparseArrays.sparse(A::LinearMap{T}) where {T}
     M, N = size(A)
-    T = eltype(A)
     rowind = Int[]
     nzval = T[]
     colptr = Vector{Int}(undef, N+1)
     v = fill(zero(T), N)
+    Av = Vector{T}(undef, M)
 
     for i = 1:N
         v[i] = one(T)
-        Lv = A * v
-        js = findall(!iszero, Lv)
+        mul!(Av, A, v)
+        js = findall(!iszero, Av)
         colptr[i] = length(nzval) + 1
         if length(js) > 0
             append!(rowind, js)
-            append!(nzval, Lv[js])
+            append!(nzval, Av[js])
         end
         v[i] = zero(T)
     end
     colptr[N+1] = length(nzval) + 1
 
-    return SparseArrays.SparseMatrixCSC(M, N, colptr, rowind, nzval)
+    return SparseMatrixCSC(M, N, colptr, rowind, nzval)
 end
 
 include("transpose.jl") # transposing linear maps
@@ -120,7 +120,7 @@ on length `N` vectors and producing length `M` vectors (with default value `N=M`
 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 (conjugate) transpose of `f`.
+function `fc` can be specified that implements the transpose/adjoint of `f`.
 
 The keyword arguments and their default values for functions `f` are
 *   issymmetric::Bool = false : whether `A` or `f` acts as a symmetric matrix

src/composition.jl

@@ -3,12 +3,12 @@ struct CompositeMap{T, As<:Tuple{Vararg{LinearMap}}} <: LinearMap{T}
     function CompositeMap{T, As}(maps::As) where {T, As}
         N = length(maps)
         for n = 2:N
-            size(maps[n],2) == size(maps[n-1],1) || throw(DimensionMismatch("CompositeMap"))
+            size(maps[n], 2) == size(maps[n-1], 1) || throw(DimensionMismatch("CompositeMap"))
         end
         for n = 1:N
             promote_type(T, eltype(maps[n])) == T || throw(InexactError())
         end
-        new{T,As}(maps)
+        new{T, As}(maps)
     end
 end
 CompositeMap{T}(maps::As) where {T, As<:Tuple{Vararg{LinearMap}}} = CompositeMap{T, As}(maps)
@@ -21,52 +21,52 @@ Base.isreal(A::CompositeMap) = all(isreal, A.maps) # sufficient but not necessar
 function LinearAlgebra.issymmetric(A::CompositeMap)
     N = length(A.maps)
     if isodd(N)
-        issymmetric(A.maps[div(N+1,2)]) || return false
+        issymmetric(A.maps[div(N+1, 2)]) || return false
     end
-    for n = 1:div(N,2)
+    for n = 1:div(N, 2)
         A.maps[n] == transpose(A.maps[N-n+1]) || return false
     end
     return true
 end
 function LinearAlgebra.ishermitian(A::CompositeMap)
     N = length(A.maps)
     if isodd(N)
-        ishermitian(A.maps[div(N+1,2)]) || return false
+        ishermitian(A.maps[div(N+1, 2)]) || return false
     end
-    for n = 1:div(N,2)
+    for n = 1:div(N, 2)
         A.maps[n] == adjoint(A.maps[N-n+1]) || return false
     end
     return true
 end
 function LinearAlgebra.isposdef(A::CompositeMap)
     N = length(A.maps)
     if isodd(N)
-        isposdef(A.maps[div(N+1,2)]) || return false
+        isposdef(A.maps[div(N+1, 2)]) || return false
     end
-    for n = 1:div(N,2)
+    for n = 1:div(N, 2)
         A.maps[n] == adjoint(A.maps[N-n+1]) || return false
     end
     return true
 end
 
 # composition of linear maps
 function Base.:(*)(A1::CompositeMap, A2::CompositeMap)
-    size(A1,2) == size(A2,1) || throw(DimensionMismatch("*"))
-    T = promote_type(eltype(A1),eltype(A2))
+    size(A1, 2) == size(A2, 1) || throw(DimensionMismatch("*"))
+    T = promote_type(eltype(A1), eltype(A2))
     return CompositeMap{T}(tuple(A2.maps..., A1.maps...))
 end
 function Base.:(*)(A1::LinearMap, A2::CompositeMap)
-    size(A1,2) == size(A2,1) || throw(DimensionMismatch("*"))
-    T = promote_type(eltype(A1),eltype(A2))
+    size(A1, 2) == size(A2, 1) || throw(DimensionMismatch("*"))
+    T = promote_type(eltype(A1), eltype(A2))
     return CompositeMap{T}(tuple(A2.maps..., A1))
 end
 function Base.:(*)(A1::CompositeMap, A2::LinearMap)
-    size(A1,2) == size(A2,1) || throw(DimensionMismatch("*"))
-    T = promote_type(eltype(A1),eltype(A2))
+    size(A1, 2) == size(A2, 1) || throw(DimensionMismatch("*"))
+    T = promote_type(eltype(A1), eltype(A2))
     return CompositeMap{T}(tuple(A2, A1.maps...))
 end
 function Base.:(*)(A1::LinearMap, A2::LinearMap)
-    size(A1,2) == size(A2,1) || throw(DimensionMismatch("*"))
+    size(A1, 2) == size(A2, 1) || throw(DimensionMismatch("*"))
     T = promote_type(eltype(A1),eltype(A2))
     return CompositeMap{T}(tuple(A2, A1))
 end
@@ -87,10 +87,10 @@ function A_mul_B!(y::AbstractVector, A::CompositeMap, x::AbstractVector)
         A_mul_B!(dest, A.maps[1], x)
         source = dest
         if N>2
-            dest = Array{T}(undef, size(A.maps[2],1))
+            dest = Array{T}(undef, size(A.maps[2], 1))
         end
         for n=2:N-1
-            resize!(dest, size(A.maps[n],1))
+            resize!(dest, size(A.maps[n], 1))
             A_mul_B!(dest, A.maps[n], source)
             dest, source = source, dest # alternate dest and source
         end

src/functionmap.jl

@@ -1,4 +1,4 @@
-struct FunctionMap{T,F1,F2} <: LinearMap{T}
+struct FunctionMap{T, F1, F2} <: LinearMap{T}
     f::F1
     fc::F2
     M::Int
@@ -12,8 +12,8 @@ function FunctionMap{T}(f::F1, fc::F2, M::Int, N::Int;
     ismutating::Bool  = _ismutating(f),
     issymmetric::Bool = false,
     ishermitian::Bool = (T<:Real && issymmetric),
-    isposdef::Bool    = false) where {T,F1,F2}
-    FunctionMap{T,F1,F2}(f, fc, M, N, ismutating, issymmetric, ishermitian, isposdef)
+    isposdef::Bool    = false) where {T, F1, F2}
+    FunctionMap{T, F1, F2}(f, fc, M, N, ismutating, issymmetric, ishermitian, isposdef)
 end
 
 # additional constructors
@@ -22,11 +22,11 @@ FunctionMap{T}(f, M::Int, N::Int; kwargs...) where {T} = FunctionMap{T}(f, nothi
 FunctionMap{T}(f, fc, M::Int; kwargs...) where {T}     = FunctionMap{T}(f, fc, M, M; kwargs...)
 
 # show
-function Base.show(io::IO, A::FunctionMap{T,F,Nothing}) where {T,F}
-    print(io,"LinearMaps.FunctionMap{$T}($(A.f), $(A.M), $(A.N); ismutating=$(A._ismutating), issymmetric=$(A._issymmetric), ishermitian=$(A._ishermitian), isposdef=$(A._isposdef))")
+function Base.show(io::IO, A::FunctionMap{T, F, Nothing}) where {T, F}
+    print(io, "LinearMaps.FunctionMap{$T}($(A.f), $(A.M), $(A.N); ismutating=$(A._ismutating), issymmetric=$(A._issymmetric), ishermitian=$(A._ishermitian), isposdef=$(A._isposdef))")
 end
 function Base.show(io::IO, A::FunctionMap{T}) where {T}
-    print(io,"LinearMaps.FunctionMap{$T}($(A.f), $(A.fc), $(A.M), $(A.N); ismutating=$(A._ismutating), issymmetric=$(A._issymmetric), ishermitian=$(A._ishermitian), isposdef=$(A._isposdef))")
+    print(io, "LinearMaps.FunctionMap{$T}($(A.f), $(A.fc), $(A.M), $(A.N); ismutating=$(A._ismutating), issymmetric=$(A._issymmetric), ishermitian=$(A._ishermitian), isposdef=$(A._isposdef))")
 end
 
 # properties
@@ -61,7 +61,7 @@ function At_mul_B!(y::AbstractVector, A::FunctionMap, x::AbstractVector)
         if !isreal(A)
             x = conj(x)
         end
-        (ismutating(A) ? A.fc(y,x) : copyto!(y, A.fc(x)))
+        (ismutating(A) ? A.fc(y, x) : copyto!(y, A.fc(x)))
         if !isreal(A)
             conj!(y)
         end

src/linearcombination.jl

@@ -1,15 +1,15 @@
 struct LinearCombination{T, As<:Tuple{Vararg{LinearMap}}, Ts<:Tuple} <: LinearMap{T}
     maps::As
     coeffs::Ts
-    function LinearCombination{T,As,Ts}(maps::As, coeffs::Ts) where {T,As,Ts}
+    function LinearCombination{T, As, Ts}(maps::As, coeffs::Ts) where {T, As, Ts}
         N = length(maps)
         N == length(coeffs) || error("Number of coefficients doesn't match number of terms")
         sz = size(maps[1])
         for n = 1:N
             size(maps[n]) == sz || throw(DimensionMismatch("LinearCombination"))
             promote_type(T, eltype(maps[n]), typeof(coeffs[n])) == T || throw(InexactError())
         end
-        new{T,As,Ts}(maps, coeffs)
+        new{T, As, Ts}(maps, coeffs)
     end
 end
 
@@ -32,7 +32,7 @@ function Base.:(+)(A1::LinearMap, A2::LinearCombination)
     T = promote_type(eltype(A1), eltype(A2))
     return LinearCombination{T}(tuple(A1, A2.maps...), tuple(one(T), A2.coeffs...))
 end
-Base.:(+)(A1::LinearCombination, A2::LinearMap) = +(A2,A1)
+Base.:(+)(A1::LinearCombination, A2::LinearMap) = +(A2, A1)
 function Base.:(+)(A1::LinearMap, A2::LinearMap)
     size(A1)==size(A2) || throw(DimensionMismatch("+"))
     T = promote_type(eltype(A1), eltype(A2))
@@ -53,7 +53,7 @@ function Base.:(-)(A1::LinearCombination, A2::LinearMap)
     T = promote_type(eltype(A1), eltype(A2))
     return LinearCombination{T}(tuple(A1.maps..., A2), tuple(A1.coeffs..., -one(T)))
 end
-function Base.:(-)(A1::LinearMap,A2::LinearMap)
+function Base.:(-)(A1::LinearMap, A2::LinearMap)
     size(A1) == size(A2) || throw(DimensionMismatch("-"))
     T = promote_type(eltype(A1), eltype(A2))
     return LinearCombination{T}(tuple(A1, A2), tuple(one(T), -one(T)))
@@ -83,7 +83,7 @@ function Base.:(\)(α::Number, A::LinearCombination)
     T = promote_type(eltype(α), eltype(A))
     return LinearCombination{T}(A.maps, map(x->α\x, A.coeffs))
 end
-Base.:(/)(A::LinearCombination, α::Number) = \(α,A)
+Base.:(/)(A::LinearCombination, α::Number) = \(α, A)
 
 # comparison of LinearCombination objects
 Base.:(==)(A::LinearCombination, B::LinearCombination) = (eltype(A)==eltype(B) && A.maps==B.maps && A.coeffs==B.coeffs)

src/transpose.jl

@@ -14,7 +14,7 @@ LinearAlgebra.adjoint(A::LinearMap{<:Real}) = transpose(A)
 LinearAlgebra.adjoint(A::LinearMap) = AdjointMap(A)
 
 # properties
-Base.size(A::Union{TransposeMap, AdjointMap}) = (size(A.lmap,2), size(A.lmap,1))
+Base.size(A::Union{TransposeMap, AdjointMap}) = (size(A.lmap, 2), size(A.lmap, 1))
 LinearAlgebra.issymmetric(A::Union{TransposeMap, AdjointMap}) = issymmetric(A.lmap)
 LinearAlgebra.ishermitian(A::Union{TransposeMap, AdjointMap}) = ishermitian(A.lmap)
 LinearAlgebra.isposdef(A::Union{TransposeMap, AdjointMap}) = isposdef(A.lmap)

src/wrappedmap.jl

@@ -8,13 +8,13 @@ function WrappedMap(lmap::Union{AbstractMatrix{T}, LinearMap{T}};
     issymmetric::Bool = issymmetric(lmap),
     ishermitian::Bool = ishermitian(lmap),
     isposdef::Bool = isposdef(lmap)) where {T}
-    WrappedMap{T,typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef)
+    WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef)
 end
 function WrappedMap{T}(lmap::Union{AbstractMatrix, LinearMap};
     issymmetric::Bool = issymmetric(lmap),
     ishermitian::Bool = ishermitian(lmap),
     isposdef::Bool = isposdef(lmap)) where {T}
-    WrappedMap{T,typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef)
+    WrappedMap{T, typeof(lmap)}(lmap, issymmetric, ishermitian, isposdef)
 end
 
 # properties

test/runtests.jl

@@ -95,7 +95,8 @@ B = LinearMap(Hermitian(rand(ComplexF64, 10, 10)))
 @test B == B'
 
 # test sparse conversions
-@test sparse(M) == SparseArrays.sparse(Array(M))
+@test sparse(M) == sparse(Array(M))
+@test convert(SparseMatrixCSC, M) == sparse(Array(M))
 
 B = copy(A)
 B[rand(1:length(A), 30)] .= 0.