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

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

2009/7/28 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>
> In case some peoples are interested by high performance on large matrices,
> I've almost finished to implement fast routines covering the complete BLAS
> API. Those include:
>
> * fast triangular * general matrix/vector (TRVM / TRMM), e.g.:
>   m3 += s1 * m1.adjoint().triangularView<UnitUpperTriangular>() * m2;
>   m3 -= s1 * m2 * m1.triangularView<LowerTriangular>();
>
> * fast triangular solver for multiple right hand sides (TRSV / TRSM):
>   m1.adjoint().triangularView<UnitUpperTriangular>().solveInPlace(m2);
>
> * fast selfadjoint * general matrix/vector (SYMV / SYMM / HEMM), e.g.:
>   m3 += s1 * m1.selfadjointView<UpperTriangular>() * m2;
>   m3 -= s1 * m2 * m1.selfadjointView<LowerTriangular>();
>
> * fast rank 1 / rank K update of a selfadjoint matrix (SYR / SYRK), e.g.:
>   m1.selfadjointView<UpperTriangular>().rankUpdate(u,s); // m1 += s * u *
> u^T
>   m1.selfadjointView<UpperTriangular>().rankUpdate(m2.adjoint(),s); // m1 +=
> s * m2.adjoint() * m2
>
> * fast rank 2 update of a selfadjoint matrix (SYR2), e.g.:
>   m1.selfadjointView<UpperTriangular>().rankUpdate(u,v,s); // m1 += s * u *
> v^T + s * v * u^T
>
>
> Some remarks/questions:
>
> * all these new routines are just high level blocking algorithms built on
> top of a single highly optimized product kernel.
>
> * they all work at full speed in all configurations (storage orders).
>
> * are you ok with the new api to deal with triangular/selfadjoint matrices ?

I'm OK, at least I cant think of a better one

>
> * do you prefer the shorter triangular<> / selfadjoint<> instead of
> triangularView<> / selfadjointView<> ?

I think the View helps potentially disambiguate

>
> * are you ok with the rankUpdate() API ? do you have any better idea to
> expose this ?

No experience with this so I cant say if that sounds natural to a
specialist. At least i cant think of  a better name.

>
> * 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.

>
> * these new special products are lazy by default.

Why? No risk of aliasing effects, somehow?

Benoit