Re: [eigen] TriangularView::solve() interface

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

As I'm going to explain ASAP in a mail to this list, I believe that we
can drop altogether the "bool success" in all solve() methods.

There remains the choice of an API between
API 1: a.solve(b,&result)
and
API 2: result = a.solve(b)

The main argument in favor of API 1 is the one Keir mentions.

But what about the wonderful little ReturnByValue class that Gael
wrote? I suggest we investigate that, at the moment I don't see why it
wouldn't work here! It would allow to have API 2 with the guarantee of
it unfolding to the same as API 1.

Maybe Gael can comment.
Meanwhile, is there consensus that, if it works, then we want it?
Again, this is assuming that the bool return value are removed, i'm
going to write to the list to advocate that.

Benoit


2009/8/25 Keir Mierle <mierle@xxxxxxxxx>:
> I had some thoughts about this, so before I forget:
>
> 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
>>
>>
>
>