Brainstorm about speeding up TJLF

[orso82]

For typical use-cases each eigensolve tends to be fast (grows with number of species and number of basis). The main issue is that TJLF does a lot of eigensolves! (order ~200)

Many of the eigensolves are used to find the width of the modes. Can such width be determined only for few Ky's and not for each one of them? For example, when width is entered by user, the same width is used for all Ky's. Perhaps we can calculate only few widths (eg at lowest and highest Kw) and interpolate between them?

[bclyons12]

Ah, okay, this is the stuff I don't know anything about. Each call to eigen allocates some workspace (I believe) and then calls a non-allocating version:

function eigen(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
    S = promote_type(eigtype(TA), TB)
    eigen!(eigencopy_oftype(A, S), eigencopy_oftype(B, S); kws...)
end

If we can preallocate the workspace and call eigen! directly, maybe that would have some improvement. I'm not sure though. It does look like most of the time is spent in ggev!

[tomneiser]

The width depends on the type of mode. Above a certain ky we can expect to be only dealing with ETG's, so there could be some improvement for say k_y > 10. Otherwise we would need to run a form of mode ID routine that guesses the best width based on inputs or initial eigenvalues.

[orso82]

@bclyons12 there could certainly be some allocation savings if we preallocated the workspace. As you know, we just need to be careful if we eventually desire to multi-thread the calculation per k_y that each thread uses its own work matrix.

[orso82]

@tomneiser should we expect those widths to change significantly between iterations of the FluxMatcher? If not, we could happily pay the price of calculating the optimal width once, and then use those same widths for all subsequent FluxMatcher iterations.

[tomneiser]

I think this is a good idea. The widths don't change every iteration, maybe every 30 iterations we can re-evaluate them