Potentially modify solving Lyapunov equations

@tholden commented on https://github.com/tholden/dynare/commit/69daaa0460b0ddee97292c39d40355201e316622 that

lyapunov_symn is a utility function provided by Dynare that's likely to be useful in a whole number of places when solving and estimating models. If I recall correctly, I had something in my own code where I needed to solve a Lyapunov equation, so I understandably called the dynare function. Since this wasn't in likelihood initialization, I didn't ever notice the lyapunov=doubling option, and the fact that it led to a different function being called.

Indeed, even if you just compare where lyapunov_symn and disclyap_fast are called within Dynare, you'll see that many bits of other code only call the former, e.g. global identification, so I was not the only person led into the trap of calling the inefficient function.

So this issue looks like it was caused by some bad design. It does not make sense for disclyap_fast to be a separate file that must be called as an alternative to lyapunov_symn, particularly since lyapunov_symn already contains code for multiple methods. Either lyapunov_symn should be a shell function with a big "switch" that dispatches to other methods (each with their own function, not just doubling), or it should directly contain those other method's code. The current set-up is terrible for code reusability.

The issue with sparse matrices is another reusability one. The period doubling algorithm is designed for large matrices, and in most reasonable applications, large matrices are sparse. The cost of a few "full" statements is essentially zero when full matrices are passed in to start with. (Incidentally, in large models, significant performance improvements arise from replacing the matrices in dr with sparse ones, so this is perhaps something dynare ought to do anyway. My dynareOBC does this by default.)

In any case, is there really any situation where the algorithm with doubling is beaten by the one without it? The algorithm without doubling seems to be dominated in all regards in my experiments, so it ought to be the default.

I would thus propose to

• integrate disclyap_fast into lyapunov_symm to have on solver only
• make lyapunov_symm able to work with sparse matrices (at least for the baseline option)
• change the baseline option in Dynare to use the doubling algorithm

@rattoma What do you think? You seem to have the most experience with big models. Would you also prefer the doubling algorithm?