I switched libmv to Eigen2 from FLENS. I got to delete quite a bit of code in the process. Also, I ported the FLENS tests over to Eigen2 to ensure the functionality we are using works; which it all did with no local modifications. I am excited about not needing local modifications! The only non-superficial change I had to make was for libmv's Nullspace() call, which finds the nullspace of a possibly wide matrix for solving homogenous equations, to pad the matrix A in Ax=0 with zero rows to make A square, when A is wide. Benoit suggested the LU's kernel() function, which I will investigate next though I am unsure if it is as stable as the SVD method.
Thoughts: The compile and errors from Eigen2 were much better than with FLENS. The fixed size and dynamically sized matrices worked together flawlessly, requiring no fixes to Eigen2. There were no surprising aliasing effects; Eigen2 Does The Right Thing with respect to possibly aliased matricies. The array type functionality is quite convenient. Eigen2 expressions are concise and readable; much more so than with FLENS where it is necessary to make explicit temporaries (rather than letting the compiler decide if temporaries are necessary). For example, it's really nice that you can do (A*x + b).norm2() inline, without the boilerplate. I expected compile times to be really bad, they weren't. One of the best things about using Eigen2 is that we get better than MKL or ATLAS performance, while eliminating dependencies on crufty old FORTRAN libraries that aren't installed everywhere. That's really fantastic; getting the right LAPACK version is hard.
Overall I am quite impressed. Great job guys!
I spent a bit of time hacking the SVD code to extract the first part (bidiagonalization) but didn't end up finishing it because I wasn't sure what the right way to avoid copying the resulting U and V matrices. If you make a function Bidiagonalize(Matrix *A, Matrix *U, Vector *diagonal, Vector *super_diagonal, Matrix *V), where U and V can be passed as NULL to avoid computing them, you can't do 'local' typedefs to remove some of the template complexity like you can with a class. So I am not sure what the right approach is here. Also, the current SVD algorithm supports avoiding computing U and V, but does not expose this; I am not sure what the best way to extend the API is.
One issue with Eigen2 is that the default formatting for output is not very nice. I may make local modifications so that the array display is matlab- and python- copy & pastable (and has aligned columns). I am not sure what benefit the current format has over a more conventional display format. Also, some of the examples in the docs are hard to understand because the example output matrices are not aligned (does that have 3 or 4 columns?).