Re: [eigen] Re: recent improvements in the products

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

ok committed with version 1.

here is an overview of my next commit: a very efficient (8/9 GFLops on
2GHz CPU) block LLT implementation making use of it:

template<typename MatrixType>
void block_llt_lo(MatrixType& m)
{
  int BlockSize = 128;

  for (int k=0; k<m.rows(); k+=BlockSize)
  {
    int bs = std::min(BlockSize, m.rows()-k);
    int rs = size - k - bs;

    Block<MatrixType,Dynamic,Dynamic> A00(m,0,0,k,k);
    Block<MatrixType,Dynamic,Dynamic> A10(m,k,0,bs,k);
    Block<MatrixType,Dynamic,Dynamic> A11(m,k,k,bs,bs);

    A00.template triangularView<LowerTriangular>().template
solveInPlace<OnTheRight>(A10); // bottleneck
    A11.template selfadjointView<LowerTriangular>().rankUpdate(A10,-1);
    ei_inplace_llt_lo(A11);
  }
}

cheers,
Gael

On Thu, Jul 30, 2009 at 3:53 PM, Benoit Jacob<jacob.benoit.1@xxxxxxxxx> wrote:
> 2009/7/30 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>> On Tue, Jul 28, 2009 at 8:24 AM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
>>> > * currently we can only solve for  X = M^1 * Y, but we also need X = Y *
>>> > M^-1. Implementation wise we just have to transpose everything, so actually,
>>> > the user can probably already do:
>>> >    m1.transpose().triangularView<UpperTriangular>(m2.transpose());
>>> >  but that's not very nice API wise. The only solution I see is to add an
>>> > option to TriangularView::solveInPlace(), e.g., solveInPlace(m2, Right);
>>>
>>> Yes and then i'd call that OnTheRight.
>>
>> ok, finally it turns out that this "side" parameter should be a
>> compile time parameter, otherwise we would have:
>>
>> if(side==OnTheRight)
>>  call version_1
>> else
>>  call version_2
>>
>> that means for a template library twice the compilation time. So:
>> 1 - it could be a template parameter of the solve function with
>> overloads for mimic a default parameter:
>>
>>    m.triangularView<UnitUpperTriangular>().solveInPlace<OnTheRight>(x);
>>
>> 2 - or we can add new solve functions:
>>
>>    m.triangularView<UnitUpperTriangular>().solveInPlaceOnTheRight(x);
>>
>> 3 - or we can define OnTheRight/OnTheLeft as instances of different
>> types allowing:
>>
>>    m.triangularView<UnitUpperTriangular>().solveInPlace(x, OnTheRight);
>>
>>
>> any preference ?
>
> I think I prefer 1 or 2, but not 3. I would reserve 3 for cases where
> the special value replaces what can otherwise be a normal parameter:
> see resize(int,int) and NoChange.
>
> Benoit
>
>
>