Matrix(op) may generate incorrect result if op * x changes x

[andreasvarga]

Matrix(op) may generate incorrect result in the case when op * x changes the values stored in x. Typical example is when op is a linear equation solver which overwrites the result in the provided right hand side. Then the repeated computation of op * ei in the following code may lead to a change of ei after each execution of op * ei. A solution would be to move

ei = zeros(eltype(op), n)

inside the loop. Otherwise, the change of x when computing op * x must be explicitly forbiden.

function Base.Matrix(op :: AbstractLinearOperator)
  (m, n) = size(op)
  A = Array{eltype(op)}(undef, m, n)
  ei = zeros(eltype(op), n)
  for i = 1 : n
    ei[i] = 1
    A[:, i] = op * ei
    ei[i] = 0
  end
  return A
end

[dpo]

In-place operators are not supported at this time, but they would be a welcome addition. Again, I would recommend creating a special type InPlaceOperator <: AbstractLinearOperator and redefining Matrix(op) for them. I'll be happy to review a pull request!

[andreasvarga]

I solved the problem by simply making a copy of the input data. Regarding PR: I never prepared a PR and I am afraid I am not enough skilled to do it now. Perhaps at a later time I will try to propose one, but now I am very busy finishing with my software project for solving matrix equations.

By the way, I have a question regarding linear operators. If L is an arbitrary linear operator, the condition Matrix(L)’ = Matrix(L’) must always be fulfilled? I defined plenty of operators, but for a few of them (e.g., those mapping the upper triangular parts of symmetric matrices), this condition is not fulfilled.

I wonder if this is a deficiency of my definitions.

[dpo]

If L is an arbitrary linear operator, the condition Matrix(L)’ = Matrix(L’) must always be fulfilled?

Yes it should be. If you find a counterexample, please open a new issue. Note also that we're currently revising the design of adjoint and transpose operators in #104.