On Mon, Aug 24, 2009 at 7:36 AM, Jitse Niesen
<jitse@xxxxxxxxxxxxxxxxx> wrote:
Hello,
While rewriting part of the tutorial, I noticed that the new TriangularView::solve() function returns the solution vector, while the other solve functions pass the solution via one of the parameters. Compare
x = A.triangularView<UpperTriangular>().solve(b);
With this interface, are we guaranteed that NVRO will work? In other words, are we sure this won't force memory copies and reallocations? With the pointer interface, I am confident that eigen will do the right thing if my matrix is the right size. With the return value, I have less control. In particular, I believe you are only guaranteed that the following:
Matrix<...> x = A.lu().solve(b)
will directly construct x instead of via a temporary. Consider:
Matrix<double, HUGE, HUGE> x = A.lu().solve(b); // Ok, NVRO works here
// ...
// ...
// form new b
x = A.lu()..solve(b) // Oops! x is already constructed, so here we must create a temporary and copy.
Many numerical problems involve solving systems with the same dimensions repeatedly; it would be a shame if there was no ability to avoid the price of large temporaries.
with
A.partialLu().solve(b, &x);
An alternative, which I am not necessarily advocating because of the reasons above, is the following:
bool success;
x = A.partialLU().solve(b, &success);
This has the benefit that success can be an optional parameter, and the API becomes uniform.
Yet another thought:
If we make the solver interfaces like the following:
bool maybe_success = A.lu().solve(b, &x);
and changed the meaning of the return value to the following:
If the return value is true, then the solve *may* have succeeded. If the return value is false, then the system *definitely* did not solve. This way we can use the same API everywhere, and ignore the return value in cases where it's known that the algorithm always returns true. Each solve() implementation will declare what behaviour it has regarding success, with the default behaviour to always return true.
Thoughts?
Keir
Was this done on purpose?
The docs on the eigen website [1] were magically updated, but some of the latex was not parsed correctly. Am I doing anything wrong? Everything works fine on my own computer.
[1] http://eigen.tuxfamily.org/dox-devel/TutorialAdvancedLinearAlgebra.html
Something to perhaps add to the to-do list: a routine to estimate the condition number of a matrix based on the LU decomposition.
Cheers,
Jitse