Power-user-friendly wrapper for LAPACK?

[jagot]

For my time propagation of PDEs, I need to find the eigenvectors/-values of a Krylov subspace matrix, every step (see e.g. my fork of Krylov.jl). As it is now, e.g. stegr! allocates the work and output arrays, every time. I wonder if there are any plans to provide for those functions who need it, a split version? I mean something in the spirit of

work = stegr_work(...)
for i = 1:steps
    ...
    stegr!(other, args, work...)
end

thereby reducing the amount of (de)allocations.

I could of course implement the functions I need in a library for my own use, but I would prefer to reduce the amount of code duplication, and maybe other people than myself would find this useful?

[andreasnoack]

I think this would be a good idea but a lot of work. It was also brought up in https://github.com/JuliaLang/julia/issues/4262.

The Julia wrapper method that actually calls LAPACK should take all needed arrays (output and workspace) as arguments and only check that they have the right size. There should then also be a method that only takes the actual input arrays as arguments an makes the allocation. In that way the functionality we have now would remain while there would be a way to avoid reallcation.

[jagot]

Ok, but how would I go about it, if I were to start on this? I have never done any development on the Julia source tree before; how do I make changes to base/linalg/lapack.jl without having to recompile the base image, or where do I find info on how to do this? Ideally, without having to restart Julia every time.

Would it be OK to do pull requests on a per-routine basis, or should I create a branch where I convert all routines (or as many as I can muster) and then do a pull request?

[JaredCrean2]

I would definitely find this useful. Anything that makes it possible to do matrix/vector operations without allocating is a step in the right direction (although I am not a core developer, so do take this with a grain of salt).

[andreasnoack]

@jagot I think it's fine to do this routine by routine or in smaller chunks. It shouldn't change the behavior for the present API so it's just incremental added functionality/optimization.

I'm not sure if it is described anywhere in the manual but the easiest way of working on a change like this is to copy the method into a new file and wrap it in a new module. Then you can just include that file over and over as you make changes to it. When you are done, you can copy the methods back.

Feel free to ask questions.

[ViralBShah]

You can also create a new package outside of Base for this to start with.

[jagot]

@andreasnoack

The Julia wrapper method that actually calls LAPACK should take all needed arrays (output and workspace) as arguments and only check that they have the right size.

Actually, I would prefer not to check their size; in my case, I would like to calculate the eigenspectrum of a tridiagonal matrix up to a certain size, say m×m, but the Krylov iteration might terminate before that. I think that if the user is willing to use this interface, he/she should also be sure that the work/output arrays are large enough.

EDIT: I realized that this use case can be solved with a combination of SubArrays and StridedVecOrMat.

[jagot]

I have created a package for this purpose and implemented the splitting for stegr!. All and any help welcome.

[andreasnoack]

I realized that this use case can be solved with a combination of SubArrays and StridedVecOrMat.

Yes exactly.

I have created a package

It's fine if this makes it easier to make progress but I'd really like to see these migrated to base when they are ready.

[ViralBShah]

A package like PowerLAPACK.jl (now 4 yeas old) is the right way forward here.