Hi Gael & Eigeners,
I would like to contribute a feature that I find lacking in Eigen's linear algebra code: Efficient condition number estimation. I have implemented Higham's improved version of Hager's L1 norm estimation algorithm described in:
N. J. Higham, "FORTRAN Codes for Estimating the One-Norm of a Real or
Complex Matrix with Applications to Condition Estimation", ACM
Transactions on Mathematical Software, Vol. 14, No. 4, December 1988,
Pages 381-396.
The algorithm performs at most 10 matrix-vector multiplies (often much less). It's main application is for condition number estimation and it is used by the xyyCON condition estimation routines in LAPACK. The motivation is that the best trade-off between speed and reliability when solving general linear systems is often to use LU with partial pivoting (which has a fast blocked implementation in Eigen) and then do this fast O(n^2) condition estimation to verify that it is "safe" to use. The L1-estimator can also be used to compute accurate forward and backward error estimates for the solutions.
My implementation uses Eigen but resides in a separate code base. I am starting to move it into Eigen itself, but rather than just submitting a patch, I would like to discuss what would be the cleanest way.
The first step that I would like you opinion on, is to how to add functionality to the non-symmetric decompositions (specifically FullPivLU and PartialPivLU) for solving the system adjoint system A^* x = rhs. Here I see two approaches:
a) Add new methods "solveAdjoint" in the LU classes that return a new type solve_adjoint_retval internal type in Solve.h. (I have this implemented already. Benoit Jacob reviewed it.)
b) Add a template parameter "bool Adjoint" (with default value false) to existing solve methods in the LU classes (and pass it to solve_retval etc.).
The second step would be to add a method "Scalar rcond()" to all the decomposition classes returning the estimate of the reciprocal condition number.
Let me know what you think/prefer.
Thanks,
Rasmus