search

mathnet / mathnet-numerics

Math.NET Numerics
http://numerics.mathdotnet.com
MIT License
3.5k stars 896 forks source link

Use gerfs to refine solution of a system of linear equations (MKL, OpenBLAS) #343

Open kjbartel opened 9 years ago

kjbartel commented 9 years ago

When solving a system of linear equations, using the MKL or OpenBLAS native providers, I have often found that the results are not accurate.

For a square matrix the system can be solved calling either, LUFactor followed by LUSolveFactored, or just LUSolve (using Matrix<T>.Solve() will call the former). Both will result in a call to getrf for LU factorisation and then getrs to solve the system. The solution could then be iteratively refined using gerfs, however MathNet doesn't currently do this.

Using gerfs can make solving a system significantly slower and it may not be needed in all cases, particularly for small matrices, but I have found for my application it must be used otherwise the results will be incorrect. The easiest solution would be to add a call to gerfs in lu_solve and lu_solve_factored in lapack.cpp. An alternative would be to add new functions (lu_refine?) so that it's only used when needed and could also return the computed error vectors. Both MKL and OpenBLAS also provide a driver routine, gesvx, which calls getrs, gerfs and gerfs, so this could be used in lu_solve or a new function could be added using gesvx to solve, refine and compute errors.

The changes are quite simple but I wanted some feedback as to the best approach before I make a pull request.

cdrnet commented 9 years ago

What would be typical criterion to control that iteration? absolute error threshold, error improvement threshold and iteration count?

kjbartel commented 9 years ago

According to the source at netlib:

Test stopping criterion. Continue iterating if 1) The residual BERR(J) is larger than machine epsilon, and 2) BERR(J) decreased by at least a factor of 2 during the last iteration, and 3) At most ITMAX iterations tried.

MKL documentation also notes:

For each right-hand side, computation of the backward error involves a minimum of 4n2 floating-point operations (for real flavors) or 16n2 operations (for complex flavors). In addition, each step of iterative refinement involves 6n2 operations (for real flavors) or 24n2 operations (for complex flavors); the number of iterations may range from 1 to 5. Estimating the forward error involves solving a number of systems of linear equations A*x = b with the same coefficient matrix A and different right hand sides b; the number is usually 4 or 5 and never more than 11. Each solution requires approximately 2n2 floating-point operations for real flavors or 8n2 for complex flavors.