Matrix-vector product with preallocation

[dpo]

Is it possible currently to specify an array for the result of A * x? In Krylov methods where many matrix-vector products will be computed, it's often possible to save time and memory by preallocating the result and overwriting it. Something like A.matvec(x, y), which should be equivalent to y = A * x. It's a little bit un-python-esque, but combined with LinearOperators, it's possible to define a custom operator that uses a preallocated array transparently, e.g.,

def OperatorWithPreallocation(A):
  y = np.empty(A.shape[0])
  def matvec(x):
    A.matvec(x, y)
    return y

  ...

and then

B = OperatorWithPreallocation(A)
B * x  # python-esque!

Currently, no.

Will implement it "soon".

Can be done without "too much work".

Do you want this for LL, CSC and CSR?

[dpo]

Ideally ;-)

Doable.

I guess you only use a contiguous y vector? @dpo

[dpo]

Is that a strong requirement? Why wouldn't it work with general indexing?

Not a strong requirement, just more work to do this efficiently. ;-)

So I guess we'll go for the more general case...

[dpo]

Sorry you asked whether y should be contiguous. I guess it's ok. I can't think of a situation where a non-contiguous y would be useful, but I suppose it will happen one day.