First version with passing minimal tests

Author: gdalle

src/adjtrans.jl

@@ -141,6 +141,25 @@ function mul!(
   conj!(res)
 end
 
+function mul!(
+  res::AbstractMatrix,
+  op::AdjointLinearOperator{T, S},
+  m::AbstractMatrix,
+  α,
+  β,
+) where {T, S}
+  p = op.parent
+  (size(m, 1) == size(p, 1) && size(res, 1) == size(p, 2) && size(m, 2) == size(res, 2)) ||
+    throw(LinearOperatorException("shape mismatch"))
+  if ishermitian(p)
+    return mul!(res, p, m, α, β)
+  elseif p.ctprod! !== nothing
+    return p.ctprod!(res, m, α, β)
+  else
+    error("Not implemented")
+  end
+end
+
 function mul!(
   res::AbstractVector,
   op::TransposeLinearOperator{T, S},
@@ -188,6 +207,25 @@ function mul!(
   conj!(res)
 end
 
+function mul!(
+  res::AbstractMatrix,
+  op::TransposeLinearOperator{T, S},
+  m::AbstractMatrix,
+  α,
+  β,
+) where {T, S}
+  p = op.parent
+  (size(m, 1) == size(p, 1) && size(res, 1) == size(p, 2) && size(m, 2) == size(res, 2)) ||
+    throw(LinearOperatorException("shape mismatch"))
+  if issymmetric(p)
+    return mul!(res, p, m, α, β)
+  elseif p.tprod! !== nothing
+    return p.tprod!(res, m, α, β)
+  else
+    error("Not implemented")
+  end
+end
+
 function mul!(
   res::AbstractVector,
   op::ConjugateLinearOperator{T, S},
@@ -200,7 +238,17 @@ function mul!(
   conj!(res)
 end
 
-# TODO: overload the above for matrices?
+function mul!(
+  res::AbstractMatrix,
+  op::ConjugateLinearOperator{T, S},
+  v::AbstractMatrix,
+  α,
+  β,
+) where {T, S}
+  p = op.parent
+  mul!(res, p, conj.(v), α, β)
+  conj!(res)
+end
 
 -(op::AdjointLinearOperator) = adjoint(-op.parent)
 -(op::TransposeLinearOperator) = transpose(-op.parent)

src/cat.jl

@@ -32,8 +32,6 @@ function hcat_ctprod!(
   mul!(view(res, (Ancol + 1):nV), B, u, α, β)
 end
 
-# TODO: overload the above for matrices?
-
 function hcat(A::AbstractLinearOperator, B::AbstractLinearOperator)
   size(A, 1) == size(B, 1) || throw(LinearOperatorException("hcat: inconsistent row sizes"))
 
@@ -93,8 +91,6 @@ function vcat_ctprod!(
   mul!(res, B, view(v, (Anrow + 1):nV), α, one(T))
 end
 
-# TODO: overload the above for matrices?
-
 function vcat(A::AbstractLinearOperator, B::AbstractLinearOperator)
   size(A, 2) == size(B, 2) || throw(LinearOperatorException("vcat: inconsistent column sizes"))
 

src/constructors.jl

@@ -102,6 +102,8 @@ op = LinearOperator(Float64, 2, 2, false, false,
                     (res, v) -> mul!(res, A, v),
                     (res, w) -> mul!(res, A', w))
 ```
+
+The 3-args `mul!` also works when applying the operator on a matrix.
 """
 function LinearOperator(
   ::Type{T},

src/kron.jl

@@ -44,8 +44,6 @@ function kron(A::AbstractLinearOperator, B::AbstractLinearOperator)
   return LinearOperator{T}(nrow, ncol, symm, herm, prod!, tprod!, ctprod!)
 end
 
-# TODO: overload the above for matrices?
-
 kron(A::AbstractMatrix, B::AbstractLinearOperator) = kron(LinearOperator(A), B)
 
 kron(A::AbstractLinearOperator, B::AbstractMatrix) = kron(A, LinearOperator(B))

src/operations.jl

@@ -32,19 +32,14 @@ function mul!(res::AbstractVector, op::AbstractLinearOperator{T}, v::AbstractVec
   end
 end
 
-function mul!(res::AbstractMatrix, op::AbstractLinearOperator{T}, m::AbstractMatrix, α, β) where {T}
-  # TODO: how to handle storage?
-  error("5-argument `mul!` is not defined between a linear operator and a matrix.")
+function mul!(res::AbstractMatrix, op::AbstractLinearOperator, m::AbstractMatrix{T}, α, β) where {T}
+  op.prod!(res, m, α, β)
 end
 
-function mul!(res::AbstractVector, op::AbstractLinearOperator, v::AbstractVector{T}) where {T}
+function mul!(res::AbstractVecOrMat, op::AbstractLinearOperator, v::AbstractVecOrMat{T}) where {T}
   mul!(res, op, v, one(T), zero(T))
 end
 
-function mul!(res::AbstractMatrix, op::AbstractLinearOperator, m::AbstractMatrix{T}) where {T}
-  mul!(res, op, m, one(T), zero(T))
-end
-
 # Apply an operator to a vector.
 function *(op::AbstractLinearOperator{T}, v::AbstractVector{S}) where {T, S}
   nrow, ncol = size(op)

src/special-operators.jl

@@ -43,8 +43,6 @@ function mulOpEye!(res, v, α, β::T, n_min) where {T}
   end
 end
 
-# TODO: overload the above for matrices?
-
 """
     opEye(T, n; S = Vector{T})
     opEye(n)
@@ -86,8 +84,6 @@ function mulOpOnes!(res, v, α, β::T) where {T}
   end
 end
 
-# TODO: overload the above for matrices?
-
 """
     opOnes(T, nrow, ncol; S = Vector{T})
     opOnes(nrow, ncol)
@@ -111,8 +107,6 @@ function mulOpZeros!(res, v, α, β::T) where {T}
   end
 end
 
-# TODO: overload the above for matrices?
-
 """
     opZeros(T, nrow, ncol; S = Vector{T})
     opZeros(nrow, ncol)
@@ -136,8 +130,6 @@ function mulSquareOpDiagonal!(res, d, v, α, β::T) where {T}
   end
 end
 
-# TODO: overload the above for matrices?
-
 """
     opDiagonal(d)
 
@@ -157,9 +149,6 @@ function mulOpDiagonal!(res, d, v, α, β::T, n_min) where {T}
   end
   res[(n_min + 1):end] .= 0
 end
-
-# TODO: overload the above for matrices?
-
 """
     opDiagonal(nrow, ncol, d)
 
@@ -184,8 +173,6 @@ function multRestrict!(res, I, u, α, β)
   res[I] = u
 end
 
-# TODO: overload the above for matrices?
-
 """
     Z = opRestriction(I, ncol)
     Z = opRestriction(:, ncol)

test/test_linop.jl

@@ -61,6 +61,18 @@ function test_linop()
         @test(norm(transpose(u) * A - transpose(u) * op) <= rtol * norm(u))
         @test(typeof(u' * op * v) <: Number)
         @test(norm(u' * A * v - u' * op * v) <= rtol * norm(u))
+
+        mv = hcat(v, -2v)
+        mu = hcat(u, -2u)
+        res_mat = similar(mu)
+        res_trans = similar(mv)
+        res_adj = similar(mv)
+        mul!(res_mat, op, mv)
+        mul!(res_trans, transpose(op), mu)
+        mul!(res_adj, op', mu)
+        @test(norm(A * mv - res_mat) <= rtol * norm(mv))
+        @test(norm(transpose(A) * mu - res_trans) <= rtol * norm(mu))
+        @test(norm(A' * mu - res_adj) <= rtol * norm(mu))
       end
 
       A3 = Hermitian(A2' * A2)