RFC: enabled selctg to reorder eigen values from within gees/gges

[sglyon]

The Schur and generalized Schur routines in LAPACK have the ability to sort the result based on a callback function (parameter named selctg) that analyzes the eigenvalues.

This RFC enables this feature for each of the gees and gges routines as well the corresponding shurfact, schurfact!, and schur methods.

I have implemented basic tests that make sure the factorization still holds, but we should write better tests to see if the reordering is happening as expected. The problem is I don't know exactly what the answer should be, so it is hard to know if this is doing the right thing. Do we have any LAPACK experts out there who could come up with a good test case?

I was considering implementing the test case for the scipy qz function, but I am not getting the same answers as scipy. Could the difference possibly be attributed to differences in LAPACK implementation or platform dependence?

CC: @cc7768

[tkelman]

See Travis failure

    From worker 2:       * linalg1
    From worker 3:       * linalg2
exception on 2: ERROR: type: selctg2c: in typeassert, expected Int32, got Int64
 in selctg2c at linalg/lapack.jl:3367
 in gees! at linalg/lapack.jl:3415
 in schurfact at linalg/factorization.jl:694
 in anonymous at no file:198
 in runtests at /private/tmp/julia/share/julia/test/testdefs.jl:5
 in anonymous at multi.jl:852
 in run_work_thunk at multi.jl:603
 in anonymous at task.jl:852
while loading linalg1.jl, in expression starting on line 25

[sglyon]

Tests are passing on my machine.

Can anyone reproduce this locally? Not sure what the problem could be.

[andreasnoack]

I think that selctg returns 0 for all pairs of a and b. Couldn't that be the problem? Do you have a reason to define selctg(a, b) = abs(a) == 0.0 ? 1 : 0 like that. E.g. I think the opposite works, i.e. selctg(a, b) = abs(a) != 0.0 ? 1 : 0.

[cc7768]

I think the original intention was to match the scipy's results, but it doesn't look like they had it for the complex version. I'm going to try the less than 1 example and see what happens.

[cc7768]

Wrote something on this gist.

I want it to grab eigenvalues with modulus less than 1 and put those in the "upper left" of the matrix. They both give me the same matrix, but there are definitely some generalized eigenvalues with modulus less than 1 (can verify with sfd_ns[:values] .* conj(sfd_ns[:values]))

[ViralBShah]

Bump. What's the next step here?

[cc7768]

@spencerlyon2

We brought in another LAPACK function that let us order the eigenvalues in the way we wanted and implemented this in PR #9701 which was enough to finish what we wanted to do at that point in time. Obviously, implementing the selctg functionality would be better in the long run.

A few weeks ago, we came across some work done by Dynare.jl that looks like they have already implemented this functionality in their own package -- See this file. Neither of us have had a chance to look at this more closely due to the semester being in full gear, but if what they have works then it is already done and it is just a matter of writing it up.

[sglyon]

As @cc7768 noted, I think we can do something similar to what was done here if we want to get this up and running.

I may need to enable this functionality soon, so I may revisit this early next week.

[andreasnoack]

I think we can do better now that we have Ref and it is also possible to avoid the global variable. With Ref we can use normal Julia functions instead of special Julia functions where the arguments are pointers. I've taken a stab in anj/gees. On that branch, you can do

julia> A = randn(5,5);

julia> abs(eigvals(A))
5-element Array{Float64,1}:
 1.72023 
 1.72023 
 1.0195  
 0.704244
 0.101243

julia> testf(x,y) = Int(hypot(x,y) < 1);

julia> abs(schurfact(A, testf)[:values])
5-element Array{Float64,1}:
 0.704244
 0.101243
 1.72023 
 1.72023 
 1.0195

[vtjnash]

Code has been reorganized quite a lot since this PR, so I am going to assume someone will re-do this if they need it