{"action":"created","issue":{"url":"https://api.github.com/repos/mathnet/mathnet-numerics/issues/343","labels_url":"https://api.github.com/repos/mathnet/mathnet-numerics/issues/343/labels{/name}","comments_url":"https://api.github.com/repos/mathnet/mathnet-numerics/issues/343/comments","events_url":"https://api.github.com/repos/mathnet/mathnet-numerics/issues/343/events","html_url":"https://github.com/mathnet/mathnet-numerics/issues/343","id":109216420,"number":343,"title":"Use gerfs to refine solution of a system of linear equations (MKL, OpenBLAS)","user":{"login":"kjbartel","id":4158630,"avatar_url":"https://avatars.githubusercontent.com/u/4158630?v=3","gravatar_id":"","url":"https://api.github.com/users/kjbartel","html_url":"https://github.com/kjbartel","followers_url":"https://api.github.com/users/kjbartel/followers","following_url":"https://api.github.com/users/kjbartel/following{/other_user}","gists_url":"https://api.github.com/users/kjbartel/gists{/gist_id}","starred_url":"https://api.github.com/users/kjbartel/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/kjbartel/subscriptions","organizations_url":"https://api.github.com/users/kjbartel/orgs","repos_url":"https://api.github.com/users/kjbartel/repos","events_url":"https://api.github.com/users/kjbartel/events{/privacy}","received_events_url":"https://api.github.com/users/kjbartel/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":2,"created_at":"2015-10-01T02:19:10Z","updated_at":"2015-10-01T05:35:22Z","closed_at":null,"body":"When solving a system of linear equations, using the MKL or OpenBLAS native providers, I have often found that the results are not accurate.\r\n\r\nFor 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`](https://software.intel.com/en-us/node/520877#E4779E02-346C-4670-92AB-C67BD8559051) for LU factorisation and then [`getrs`](https://software.intel.com/en-us/node/520892#642A8C07-088C-408D-BC89-D0F2A6E75416) to solve the system. The solution could then be iteratively refined using [`gerfs`](https://software.intel.com/en-us/node/520926#B82F227A-98F0-4453-AAD4-1EF922954E2A), however MathNet doesn't currently do this.\r\n\r\nUsing `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`](https://software.intel.com/en-us/node/520974) to solve, refine and compute errors.\r\n\r\nThe changes are quite simple but I wanted some feedback as to the best approach before I make a pull request."},"comment":{"url":"https://api.github.com/repos/mathnet/mathnet-numerics/issues/comments/144625951","html_url":"https://github.com/mathnet/mathnet-numerics/issues/343#issuecomment-144625951","issue_url":"https://api.github.com/repos/mathnet/mathnet-numerics/issues/343","id":144625951,"user":{"login":"kjbartel","id":4158630,"avatar_url":"https://avatars.githubusercontent.com/u/4158630?v=3","gravatar_id":"","url":"https://api.github.com/users/kjbartel","html_url":"https://github.com/kjbartel","followers_url":"https://api.github.com/users/kjbartel/followers","following_url":"https://api.github.com/users/kjbartel/following{/other_user}","gists_url":"https://api.github.com/users/kjbartel/gists{/gist_id}","starred_url":"https://api.github.com/users/kjbartel/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/kjbartel/subscriptions","organizations_url":"https://api.github.com/users/kjbartel/orgs","repos_url":"https://api.github.com/users/kjbartel/repos","events_url":"https://api.github.com/users/kjbartel/events{/privacy}","received_events_url":"https://api.github.com/users/kjbartel/received_events","type":"User","site_admin":false},"created_at":"2015-10-01T05:35:21Z","updated_at":"2015-10-01T05:35:21Z","body":"According to the source at [netlib](http://www.netlib.org/lapack/explore-html/df/d32/dgerfs_8f_source.html):\r\n\r\n> Test stopping criterion. Continue iterating if\r\n> 1) The residual BERR(J) is larger than machine epsilon, and\r\n> 2) BERR(J) decreased by at least a factor of 2 during the last iteration, and\r\n> 3) At most ITMAX iterations tried.\r\n\r\nMKL [documentation also notes](https://software.intel.com/en-us/node/520926#B82F227A-98F0-4453-AAD4-1EF922954E2A):\r\n\r\n> 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."}}