don't use mul! inplace (#482)

* don't use mul! inplace

* accept temp args in inplace mul_coefficients

Author: jishnub

src/Operators/SubOperator.jl

@@ -379,14 +379,16 @@ function mul_coefficients(A::SubOperator{<:Any,<:Any,NTuple{2,UnitRange{Int}}},
     if size(A,2) == length(b)
         AbstractMatrix(A)*b
     else
-        view(A,:,1:length(b))*b
+        AbstractMatrix(view(A,:,1:length(b)))*b
     end
 end
-function mul_coefficients!(A::SubOperator{<:Any,<:Any,NTuple{2,UnitRange{Int}}}, b)
+function mul_coefficients!(A::SubOperator{<:Any,<:Any,NTuple{2,UnitRange{Int}}}, b,
+        temp = similar(b, promote_type(eltype(A), eltype(b)), size(A,1)))
     if size(A,2) == length(b)
-        mul!(b, AbstractMatrix(A), b)
+        mul!(temp, AbstractMatrix(A), b)
     else
-        mul!(b, view(A,:,1:length(b)), b)
+        mul!(temp, AbstractMatrix(view(A,:,1:length(b))), b)
     end
+    b .= temp
     return b
 end

src/Operators/banded/ConstantOperator.jl

@@ -78,7 +78,7 @@ function mul_coefficients(A::ConstantOperator, b)
         A.λ * b
     end
 end
-function mul_coefficients!(A::ConstantOperator, b)
+function mul_coefficients!(A::ConstantOperator, b, args...)
     λ = A.λ
     if !isone(λ)
         b .*= λ

src/Operators/general/algebra.jl

@@ -694,27 +694,27 @@ end
 ## Operations
 for mulcoeff in [:mul_coefficients, :mul_coefficients!]
     @eval begin
-        function $mulcoeff(A::Operator, b)
-            mulcoeff_maybeconvert($mulcoeff, A, b)
+        function $mulcoeff(A::Operator, b, args...)
+            mulcoeff_maybeconvert($mulcoeff, A, b, args...)
         end
 
-        function $mulcoeff(A::TimesOperator, b)
+        function $mulcoeff(A::TimesOperator, b, args...)
             ret = b
             for k = length(A.ops):-1:1
-                ret = $mulcoeff(A.ops[k], ret)
+                ret = $mulcoeff(A.ops[k], ret, args...)
             end
 
             ret
         end
     end
 end
 
-function _mulcoeff_maybeconvert(f::F, A, b) where {F}
+function _mulcoeff_maybeconvert(f::F, A, b, args...) where {F}
     n = size(b, 1)
-    n > 0 ? f(view(A, FiniteRange, 1:n), b) : b
+    n > 0 ? f(view(A, FiniteRange, 1:n), b, args...) : b
 end
 
-mulcoeff_maybeconvert(f, A, b) = _mulcoeff_maybeconvert(f, A, b)
+mulcoeff_maybeconvert(args...) = _mulcoeff_maybeconvert(args...)
 
 function mulcoeff_maybeconvert(f::typeof(mul_coefficients), A, b::Vector)
     v = _mulcoeff_maybeconvert(f, A, b)

src/Operators/ldiv.jl

@@ -47,7 +47,7 @@ end
 \(A::Operator,B::MatrixFun;kwds...) = \(A,Array(B);kwds...)
 
 ldiv_coefficients(A::Operator,b;kwds...) = ldiv_coefficients(qr(A),b;kwds...)
-ldiv_coefficients!(A::Operator,b;kwds...) = ldiv_coefficients!(qr(A),b;kwds...)
+ldiv_coefficients!(A::Operator,b,args...;kwds...) = ldiv_coefficients!(qr(A),b,args...;kwds...)
 
 \(A::Operator,B::Operator) = TimesOperator(inv(A),B)
 

src/Operators/qr.jl

@@ -143,14 +143,18 @@ det(A::Operator) = det(qr(A))
 
 mul_coefficients(At::Transpose{T,<:QROperatorQ{T}},B::AbstractVector{T}) where {T<:Real} = parent(At)'*B
 mul_coefficients(At::Transpose{T,<:QROperatorQ{T}},B::AbstractMatrix{T}) where {T<:Real} = parent(At)'*B
-mul_coefficients!(At::Transpose{T,<:QROperatorQ{T}},B::AbstractVector{T}) where {T<:Real} = mul!(B, parent(At)', B)
-mul_coefficients!(At::Transpose{T,<:QROperatorQ{T}},B::AbstractMatrix{T}) where {T<:Real} = mul!(B, parent(At)', B)
+function mul_coefficients!(At::Transpose{T,<:QROperatorQ{T}},B::AbstractVector{T},args...) where {T<:Real}
+    mul_coefficients!(parent(At)', B, args...)
+end
+function mul_coefficients!(At::Transpose{T,<:QROperatorQ{T}},B::AbstractMatrix{T},args...) where {T<:Real}
+    mul_coefficients!(parent(At)', B, args...)
+end
 
 function mul_coefficients(Ac::Adjoint{T,<:QROperatorQ{<:Any,T}},B::AbstractVector{T};
             tolerance=eps(eltype(Ac))/10,maxlength=1000000) where {T}
     mulpars(Ac,B,tolerance,maxlength)
 end
-function mul_coefficients!(Ac::Adjoint{T,<:QROperatorQ{<:Any,T}},B::AbstractVector{T};
+function mul_coefficients!(Ac::Adjoint{T,<:QROperatorQ{<:Any,T}},B::AbstractVector{T}, args...;
             tolerance=eps(eltype(Ac))/10,maxlength=1000000) where {T}
     mulpars(Ac,B,tolerance,maxlength,Val(true))
 end
@@ -171,7 +175,7 @@ end
 
 
 ldiv_coefficients(A::QROperatorQ, B; opts...) = mul_coefficients(A', B; opts...)
-ldiv_coefficients!(A::QROperatorQ, B; opts...) = mul_coefficients!(A', B; opts...)
+ldiv_coefficients!(A::QROperatorQ, B, args...; opts...) = mul_coefficients!(A', B, args...; opts...)
 \(A::QROperatorQ, B::Fun; opts...) = *(A', B; opts...)
 
 
@@ -187,7 +191,7 @@ function ldiv_coefficients(R::QROperatorR, b::AbstractVector)
     U = uppertriangularview!(R, length(b))
     U \ b
 end
-function ldiv_coefficients!(R::QROperatorR, b::AbstractVector)
+function ldiv_coefficients!(R::QROperatorR, b::AbstractVector, args...)
     U = uppertriangularview!(R, length(b))
     ldiv!(U, b)
 end
@@ -226,8 +230,8 @@ ldiv_coefficients(QR::QROperator{<:Any,<:Any,<:Complex}, b::AbstractVector{<:Rea
     ldiv_coefficients(QR, Vector{eltype(QR)}(b);kwds...)
 
 # fallback methods that may not be very efficient
-function ldiv_coefficients!(QR::QROperator{<:Any,<:Any,<:Real}, b::AbstractVector{<:Number}; kwds...)
-    ldiv_coefficients!(convert(Operator{eltype(b)}, QR), b)
+function ldiv_coefficients!(QR::QROperator{<:Any,<:Any,<:Real}, b::AbstractVector{<:Number}, args...; kwds...)
+    ldiv_coefficients!(convert(Operator{eltype(b)}, QR), b, args...; kwds...)
 end
 
 \(A::QROperator,b::Fun;kwds...) =