Re: [eigen] still the solve() API debate

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

I have more arguments to debate about the in-place versions (for the
general solver function I still think that opt 3 is a good one).

The ability to deal with sub matrices is much more important for the
in-place versions than for the general ones, and in particular this is
very useful to implement blocked algorithms. That's why I suggested
..re(). However:

- It is very likely that you we'll reuse this sub matrix to perform
some other computations, => -1 for .ref() (e.g. see the implementation
of the blocked LLT or blocked LU)

- but... to make the writing of blocked algorithms easier, I planed to
write a small Partitioner class allowing to easily access the current
and surrounding blocks. In that case the Partitioner class will gives
you temporary Block expressions => + 1 for .ref()

To be more precise, e.g., in the blocked version of LLT, instead of writing:

Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);

llt(A11);
A11.adjoint().template triangularView<UpperTriangular>().template
solveInPlace<OnTheRight>(A21);
A22.template selfadjointView<LowerTriangular>().rankUpdate(A21,-1);

you will declare a partitioner:

Partitioner<MatrixType> P(m, ...);

and then use the shortcuts provided by P:

llt(P.B11());
P.B11().adjoint().template triangularView<UpperTriangular>().template
solveInPlace<OnTheRight>(P.B21().ref());
P.B22().template selfadjointView<LowerTriangular>().rankUpdate(P.B21(),-1);

I know that's another story, but the problem is that now I'm really
puzzled, so please add more arguments into the balance !

gael


On Fri, Sep 11, 2009 at 10:25 AM, Gael Guennebaud
<gael.guennebaud@xxxxxxxxx> wrote:
> First of all I don't like the use of pointers here, because it is not
> a standard c++ practice, and therefore users we'll start wondering: "
> why do they use a pointer here ? maybe because this parameter is
> optional ? hm does not make sense... unless if x==0, then the
> computation is performed in place overwritten b ? or the solve
> function creates the matrix object and then I have to delete it myself
> ? etc. ". OK, maybe I'm exaggerating a bit, but that's the idea.
>
> ok, now let's summarize, ideally we want an API which satisfies the
> following constraints:
>
> a - concise (not too verbose)
> b - readable
> c - simple behavior
> d - it should return the result so that we can directly use it in an expression
> e - does not create unnecessary temporary
> f - directly works with temporary proxy object
> g - compliant with standard c++ practices
>
> since I changed my mind about ReturnByValue inheriting MatrixBase, I
> think opt 3 is the one which best satisfies all these constraints.
>
>
> However, I (we) forgot to include in the discussion the API of the
> in-place solve functions which should be as close as possible to the
> generic solve().  Currently this only apply to the triangular solver,
> LLT and LDLT, but I think this could also apply to PartialLU. Here the
> problem is a bit different because the result is also and input and so
> only option 1 and 4 could work:
>
> opt 1: A.solveInPlace(&b_and_x);
>
> opt 4: A.solveInPlace(b_and_x.ref());
>
> Since opt 1 does not satisfy the constraint f and that this constraint
> is particularly important for the in-place triangular solver, it
> currently remains the opt 4 only.
>
> In order to satisfy the constraint d, the solveInPlace functions could
> return a reference to the inout argument. So at that point one could
> do either:
>
> Z = 2 * M * A.solveInPlace(b.ref());    // no temporary, and b is overwritten
>
> Z = 2 * M * A.solve(b);    // one temporary but b is preserved
>
> So yes you got me right, I'm suggesting to rely on both the
> ReturnByValue<> and Reference<> mechanisms ! and since they have not
> been extensively tested yet I cannot promise that we won't discover
> any issue with them (especially the Reference<> one).
>
>
>
> Also, note that a variant of opt 4, is to replace .ref() by the (more
> explicit ?) pair .out() and .inout(), in which case the "InPlace"
> suffix may be removed, e.g.:
>
> Z = 2 * M * A.solve(b.inout())
>
> but note that how subtle the difference is with:
>
> Z = 2 * M * A.solve(b);
>
> so maybe it is better to stuck with .ref() and the InPlace suffix.
>
>
> gael.
>
> On Fri, Sep 11, 2009 at 1:10 AM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
>> 2009/9/10 Keir Mierle <mierle@xxxxxxxxx>:
>>> On Thu, Sep 10, 2009 at 2:47 PM, Gael Guennebaud <gael.guennebaud@xxxxxxxxxx>
>>> wrote:
>>>>
>>>> oh, actually we can even write a single
>>>> SolveReturnType<DecomositionType, DerivedB> class inheriting
>>>> ReturnByValue which would be used for all x = solve(b) functions.
>>>
>>> I am not convinced the added complexity justifies the benefits. What about
>>> using the result in an expression?
>>> y = M*A.solve(b); // Parts of my code could use this today.
>>
>> My understanding is that Gael is now proposing a ReturnByValue class
>> that inherits MatrixBase, so that one could use that in any
>> expression. That's a departure from my assumption that he was
>> proposing ReturnByValue not inheriting MatrixBase.
>>
>> Gael : is that right?
>>
>>> x.start(3) = A.solve(b) // Never used this.
>>> Not saying that the 2nd use case isn't justified, just that I personally
>>> haven't needed it. Can others speak to this?
>>
>> Same here, I don't think use cases are common. I can think of one:
>> when implementing a blocked algorithm (e.g. block LU), sometimes some
>> of the blocks can be obtained as a matrix solve.
>>
>> However, for such complex use cases, I think that it wasn't too bad
>> (although not wonderful either) to force using named xpr objects
>> instead of temporaries.
>>
>> As Keir notes, here we're giving a lot of weight to a quite unusual
>> use case. I don't think it's worth going for option 4 just for that
>> use case, as that would make the API heavier for everybody just to
>> make it lighter in that unusual use case.
>>
>>>> On Thu, Sep 10, 2009 at 11:37 PM, Gael Guennebaud
>>>> <gael.guennebaud@xxxxxxxxx> wrote:
>>>> > you forgot the main advantage of sol 3 and 4 that is to be able to use
>>>> > temporary proxy object without having to declare them, e.g.:
>>>> >
>>>> > opt 3:  x.start(n) = A.solve(b);        // works without temporary
>>
>> You're right, I was forgetting about that unusual use case, but again,
>> i've never met it outside of the implementation of blocked algorithms.
>>
>> Now it's true that it has intrinsic value, and that it's a good goal
>> for our API design to fully support the seamless implementation of
>> blocked algorithms.
>>
>> My main concern, like Keir, was that this option 3 wouldn't allow the
>> result of solve() to be used in expressions, because you earlier said
>> that ReturnByValue shouldn't inherit MatrixBase. Now I didn't think as
>> much as you did about this, but now you say,
>>
>>>> >
>>>> > So why not giving a try to opt. 3 ? Indeed, actually the multiple
>>>> > matrix products are now implemented using something very very similar
>>>> > to opt 3 with a ReturnByValue inheriting MatrixBase, and they are not
>>>> > that heavy.
>>
>> Ah OK, so now we're again discussing a ReturnByValue that inherits
>> MatrixBase. So option 3 indeed looks pretty good, can be used in any
>> expression.
>>
>>>> > What I mean is that all the limitations of opt.3 with a
>>>> > ReturnByValue inheriting MatrixBase already exists with the matrix
>>>> > products...  (the limitation is tha the NestByValue flag *must* be
>>>> > honored, otherwise you have a weird (currently) compilation error).
>>
>> OK, that's not too bad, especially as, as you said, this issue already
>> exists with all the products (and there it really is the best solution
>> anyway).
>>
>> At least we'll have some uniformity. What's needed then is a good dox
>> page on writing functions that take MatrixBase as argument.
>>
>>>> > Now, honestly I don't really like to use pointers for return types,
>>>> > unless they are optional.
>>
>> Do you mean: as opposed to natural syntax, result = function() ? As in option 3.
>> Then I can understand. Me too, i prefer natural syntax!
>>
>> Or do you mean: as opposed to function(result.ref()) or passing by c++
>> reference? Then I need to ask, why? What is bad about pointer, as
>> opposed to references or custom structs? OK custom structs allow
>> passing temporaries, but are you also saying that you prefer c++
>> references over pointers? Then I really would like to know why.
>>
>> Benoit
>>
>>
>>
>