WIP: use mutating matrix multiplication in eigs

[Jutho]

Currently, the arpack wrapper of eigs allocates way too much memory. Firstly, no mutating methods are used, so every matrix multiplication creates (allocates) a new output, which is then copied to the position requested by Arpack. Secondly, for every new matrix multiplication, a new inputvector x=workd[load_idx] is created (allocated). So that is two vector allocations per matrix multiplication, both of which could be avoided. The second one will automatically go away when indexing creates a view. The first one requires to switch to mutating methods whenever possible (e.g. whenever defined for the type of A).

This is a work in progress, since so far only the standard mode is covered (i.e. no shift-invert, no generalised eigenvalue problems, so this will fail the tests). Also, I currently use pointer and pointer_to_array to create a view, but a simple slice would probably have been ok.

[ViralBShah]

Thanks for tackling this. This is much needed!

[ViralBShah]

Does using SubArray not work instead of pointer_to_array? It seems that the various mul methods do work with SubArray, so it should work out.

[Jutho]

Yes it would. I was just doing a lot of pointer work in some code of mine, where sub or slice would not work, and while doing so I also changed this and did not even think of using SubArray.

I will fix this if I complete this PR. I guess there is no need to hurry with this, i.e. it will be a 0.4 thing, together with eigs accepting function handles? I don't have much time for it right now. I just patched the standard case now since that's what I needed.

[ViralBShah]

Yes, definitely 0.4. Main things are this performance enhancement, operator functions instead of matrix inputs, and then svds. That would complete this whole feature set. Once we have multi-threading, the matvecs could be easily parallel by default.

[ViralBShah]

@Jutho Do you think this PR should be updated to use SubArrays or should we merge it as is? I am assuming the performance difference is notable.

[Jutho]

This is far from finished. If I have some time, I will try to finish it, but it will probably not be until next week.

[ViralBShah]

I believe that with the recent subarray work and indexing operators defaulting to views, this should no longer be necessary.

[Jutho]

Correct.

[Jutho]

Although, with the current implementation, every matrix vector multiplication will still allocate a new vector for the result (but no longer for the input, which is obtained by indexing into a larger array), and only then copy it to the correct position. So there would still be the benefit of using A_mul_B! for matrices and these types which support it.

[ViralBShah]

I see - matvecA will allocate. We really must use A_mul_B! and A_div_B!.

[ViralBShah]

Some svds discussion for reference: https://groups.google.com/forum/#!topic/julia-users/YEg_k2JjtmY

[Jutho]

This is on my todo list. I'll make some time for this after I have finished the improved cat PR. It will also be interesting to see how much can still be gained from this now that we have the new GC.

[ViralBShah]

I was wondering exactly the same thing myself about the new GC. Perhaps it is worth setting up some performance benchmarks, which will otherwise be good to have too - before doing a whole bunch of work here.

[ViralBShah]

For this problem, I see almost the same time and memory allocation:

Julia 0.3: 28 seconds (5.4GB allocated)
Julia 0.4-dev: 25 seconds (5.1GB allocated)

julia> s = sprand(10^5,10^5,5/10^5);
julia> @time eigs(s);
elapsed time: 25.280660173 seconds (5181 MB allocated, 1.70% gc time in 234 pauses with 0 full sweep)

Octave took 24 seconds, but I am not sure if Julia and octave are running the same number of iterations, and are comparable.

[tkelman]

closing in favor of #16294