Merge pull request #175 from vpuri3/resize

overload Base.resize<bang>

Author: xtalax

src/basic.jl

@@ -21,6 +21,7 @@ Base.size(ii::IdentityOperator) = (ii.len, ii.len)
 Base.adjoint(A::IdentityOperator) = A
 Base.transpose(A::IdentityOperator) = A
 Base.conj(A::IdentityOperator) = A
+
 LinearAlgebra.opnorm(::IdentityOperator, p::Real=2) = true
 for pred in (
              :issymmetric, :ishermitian, :isposdef,
@@ -230,6 +231,7 @@ for op in (
     @eval Base.$op(L::ScaledOperator) = ScaledOperator($op(L.λ), $op(L.L))
 end
 Base.conj(L::ScaledOperator) = conj(L.λ) * conj(L.L)
+Base.resize!(L::ScaledOperator, n::Integer) = (resize!(L.L, n); L)
 LinearAlgebra.opnorm(L::ScaledOperator, p::Real=2) = abs(L.λ) * opnorm(L.L, p)
 
 getops(L::ScaledOperator) = (L.λ, L.L,)
@@ -377,6 +379,12 @@ for op in (
     @eval Base.$op(L::AddedOperator) = AddedOperator($op.(L.ops)...)
 end
 Base.conj(L::AddedOperator) = AddedOperator(conj.(L.ops))
+function Base.resize!(L::AddedOperator, n::Integer)
+    for op in L.ops
+        resize!(op, n)
+    end
+    L
+end
 
 getops(L::AddedOperator) = L.ops
 islinear(L::AddedOperator) = all(islinear, getops(L))
@@ -509,6 +517,19 @@ for op in (
                                                           )
 end
 Base.conj(L::ComposedOperator) = ComposedOperator(conj.(L.ops); cache=L.cache)
+function Base.resize!(L::ComposedOperator, n::Integer)
+
+    for op in L.ops
+        resize!(op, n)
+    end
+
+    for v in L.cache
+        resize!(v, n)
+    end
+
+    L
+end
+
 LinearAlgebra.opnorm(L::ComposedOperator) = prod(opnorm, L.ops)
 
 getops(L::ComposedOperator) = L.ops
@@ -667,6 +688,13 @@ Base.size(L::InvertedOperator) = size(L.L) |> reverse
 Base.transpose(L::InvertedOperator) = InvertedOperator(transpose(L.L); cache = iscached(L) ? L.cache' : nothing)
 Base.adjoint(L::InvertedOperator) = InvertedOperator(adjoint(L.L); cache = iscached(L) ? L.cache' : nothing)
 Base.conj(L::InvertedOperator) = InvertedOperator(conj(L.L); cache=L.cache)
+function Base.resize!(L::InvertedOperator, n::Integer)
+
+    resize!(L.L, n)
+    resize!(L.cache, n)
+
+    L
+end
 
 getops(L::InvertedOperator) = (L.L,)
 islinear(L::InvertedOperator) = islinear(L.L)

src/func.jl

@@ -291,6 +291,23 @@ function Base.inv(L::FunctionOperator)
                     )
 end
 
+function Base.resize!(L::FunctionOperator, n::Integer)
+
+    for op in getops(L)
+        if static_hasmethod(resize!, typeof((op, n)))
+            resize!(op, n)
+        end
+    end
+
+    for v in L.cache
+        resize!(v, n)
+    end
+
+    L.traits = (; L.traits..., size = (n, n),)
+
+    L
+end
+
 function LinearAlgebra.opnorm(L::FunctionOperator, p)
     L.traits.opnorm === nothing && error("""
       M.opnorm is nothing, please define opnorm as a function that takes one

src/interface.jl

@@ -243,6 +243,10 @@ function Base.getindex(L::AbstractSciMLOperator, I::Vararg{Int, N}) where {N}
     convert(AbstractMatrix,L)[I...]
 end
 
+function Base.resize!(L::AbstractSciMLOperator, n::Integer)
+    throw(MethodError(resize!, typeof.((L, n))))
+end
+
 LinearAlgebra.exp(L::AbstractSciMLOperator) = exp(Matrix(L))
 LinearAlgebra.opnorm(L::AbstractSciMLOperator, p::Real=2) = opnorm(convert(AbstractMatrix,L), p)
 for pred in (

src/left.jl

@@ -32,8 +32,8 @@ end
 
 function LinearAlgebra.ldiv!(u::AbstractVecOrMat, L::AbstractSciMLOperator)
     op = (u isa Transpose) ? transpose : adjoint
-    ldiv!(op(v), op(L), op(u))
-    v
+    ldiv!(op(L), op(u))
+    u
 end
 
 ###
@@ -60,7 +60,7 @@ AbstractAdjointVecOrMat    = Adjoint{  T,<:AbstractVecOrMat} where{T}
 AbstractTransposedVecOrMat = Transpose{T,<:AbstractVecOrMat} where{T}
 
 has_adjoint(::AdjointOperator) = true
-#has_adjoint(::TransposedOperator) = ??
+has_adjoint(L::TransposedOperator) = isreal(L) & has_adjoint(L.L)
 
 islinear(L::AdjointOperator) = islinear(L.L)
 islinear(L::TransposedOperator) = islinear(L.L)
@@ -79,6 +79,7 @@ for (op, LType, VType) in (
 
     # traits
     @eval Base.size(L::$LType) = size(L.L) |> reverse
+    @eval Base.resize!(L::$LType, n::Integer) = (resize!(L.L, n); L)
     @eval Base.$op(L::$LType) = L.L
 
     @eval getops(L::$LType) = (L.L,)

src/matrix.jl

@@ -172,6 +172,7 @@ Base.size(L::InvertibleOperator) = size(L.F)
 Base.transpose(L::InvertibleOperator) = InvertibleOperator(transpose(L.F))
 Base.adjoint(L::InvertibleOperator) = InvertibleOperator(L.F')
 Base.conj(L::InvertibleOperator) = InvertibleOperator(conj(L.F))
+Base.resize!(L::InvertibleOperator, n::Integer) = (resize!(L.F, n); L)
 LinearAlgebra.opnorm(L::InvertibleOperator{T}, p=2) where{T} = one(T) / opnorm(L.F)
 LinearAlgebra.issuccess(L::InvertibleOperator) = issuccess(L.F)
 
@@ -278,6 +279,15 @@ islinear(::AffineOperator) = false
 
 Base.size(L::AffineOperator) = size(L.A)
 Base.iszero(L::AffineOperator) = all(iszero, getops(L))
+function Base.resize!(L::AffineOperator, n::Integer)
+
+    resize!(L.A, n)
+    resize!(L.B, n)
+    resize!(L.b, n)
+
+    L
+end
+
 has_adjoint(L::AffineOperator) = all(has_adjoint, L.ops)
 has_mul(L::AffineOperator) = has_mul(L.A)
 has_mul!(L::AffineOperator) = has_mul!(L.A)

test/basic.jl

@@ -37,6 +37,7 @@ K = 12
     @test size(Id) == (N, N)
     @test Id' isa IdentityOperator
     @test isconstant(Id)
+    @test_throws MethodError resize!(Id, N)
 
     for op in (
                *, \,
@@ -69,6 +70,8 @@ end
     @test islinear(Z)
     @test NullOperator(u) isa NullOperator
     @test isconstant(Z)
+    @test_throws MethodError resize!(Z, N)
+
     @test zero(A) isa NullOperator
     @test convert(AbstractMatrix, Z) == zeros(size(Z))
 

test/matrix.jl

@@ -39,6 +39,8 @@ K = 19
     @test isconstant(FF)
     @test isconstant(FFt)
 
+    @test_throws MethodError resize!(AA, N)
+
     @test eachindex(A)  === eachindex(AA)
     @test eachindex(A') === eachindex(AAt) === eachindex(MatrixOperator(At))
 
@@ -101,6 +103,7 @@ end
 
     L = DiagonalOperator(d)
     @test isconstant(L)
+    @test_throws MethodError resize!(L, N)
 
     @test issquare(L)
     @test islinear(L)

test/total.jl

@@ -97,4 +97,68 @@ end
     v=rand(N2,K); @test mul!(v, op, u) ≈ op * u
     v=rand(N2,K); w=copy(v); @test mul!(v, op, u, α, β) ≈ α*(op * u) + β * w
 end
+
+@testset "Resize! test" begin
+    M1 = 4
+    M2 = 12
+
+    u = rand(N)
+    u1 = rand(M1)
+    u2 = rand(M2)
+
+    f(u, p, t) = 2 * u
+    f(v, u, p, t) = (copy!(v, u); lmul!(2, v))
+
+    fi(u, p, t) = 0.5 * u
+    fi(v, u, p, t) = (copy!(v, u); lmul!(0.5, v))
+
+    F = FunctionOperator(f, u, u; islinear = true, op_inverse = fi, issymmetric = true)
+
+    multest(L, u) = @test mul!(zero(u), L, u) ≈ L * u
+
+    function multest(L::SciMLOperators.AdjointOperator, u)
+        @test mul!(adjoint(zero(u)), adjoint(u), L) ≈ adjoint(u) * L
+    end
+
+    function multest(L::SciMLOperators.TransposedOperator, u)
+        @test mul!(transpose(zero(u)), transpose(u), L) ≈ transpose(u) * L
+    end
+
+    function multest(L::SciMLOperators.InvertedOperator, u)
+        @test ldiv!(zero(u), L, u) ≈ L \ u
+    end
+
+    for (L, LT) in (
+                    (F, FunctionOperator),
+                    (F + F, SciMLOperators.AddedOperator),
+                    (F * 2, SciMLOperators.ScaledOperator),
+                    (F ∘ F, SciMLOperators.ComposedOperator),
+                    (AffineOperator(F, F, u), AffineOperator),
+                    (SciMLOperators.AdjointOperator(F), SciMLOperators.AdjointOperator),
+                    (SciMLOperators.TransposedOperator(F), SciMLOperators.TransposedOperator),
+                    (SciMLOperators.InvertedOperator(F), SciMLOperators.InvertedOperator),
+                    (SciMLOperators.InvertibleOperator(F), SciMLOperators.InvertibleOperator),
+                   )
+
+        @info "$LT"
+
+        L = deepcopy(L)
+        L = cache_operator(L, u)
+
+        @test L isa LT
+        @test size(L) == (N, N)
+        multest(L, u)
+
+        resize!(L, M1); @test size(L) == (M1, M1)
+        multest(L, u1)
+
+        resize!(L, M2); @test size(L) == (M2, M2)
+        multest(L, u2)
+
+    end
+
+    # InvertedOperator
+    # AffineOperator
+    # FunctionOperator
+end
 #