Dear George,
the next step in the C implementation of the estimation routine is to have a compiled version of the Lyapunov equation solver. The corresponding Matlab code is in lyapunov_symm.m or around lines 185-242 of DsgeLikelihood.m. I probably got the idea of the algorithm in "Numerical Recepees in C++" or "Matrix Computations", this last one by Golub and van Loan.
lyapunov_symm.m returns the covariance matrix (and the list of stationary variables) while in DsgeLikelihood, we keep the Schur decomposition of the state space transition matrix in order to have a quasi-triangular transition matrix in the Kalman filter (matrix T) and the Schur vectors to build matrix Z. Currently we only use it for the diffuse filter, but there is no reason not to use for the regular filter if it is more efficient.
I see two ways to implement it. Write it from scratch in C++ using direct calls to the necessary Lapack functions and/or use RECSY http://www.cs.umu.se/~isak/recsy/ (we need of course to evaluate it first).
The problem/advantage of RECSY is that it uses parallel computing via OpenMP, but currently (as of 2008 versions) MKL distributed with Matlab is incompatible with OpenMP (see the discussion lead by Stephane Adjemian a few months ago).
My feeling is that building the estimation DLL will be easier with a Lyapunov solver built from scratch and that we can keep dealing with RESY for further improvments of efficiency.
Note that there must be a Lyapunov solver in dynare++, but I'm afraid that it will be plagued with the same problems as the Kalman filter and my feeling is that it takes longer to modify it than to write clean code from the start.
Please tell me your thoughts,
Best
Michel