Re: [eigen] Problem with LU and Cholesky inversion

[ Thread Index | Date Index | More Archives ]

Quoting Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:

Hi all,

as I said the problem is in solveTriangular(), there is a FIXME. In fact, it
assumes the right hand side is column major. If it is not we should copy it
to a column major one (allows much faster code), but I don't have time for
that right now....

Aaah, I didn't understand that it assumed column-major. That explains everything, thanks.

I sure can fix it if you don't have time, during this weekend (I also can't do it now).

Meanwhile, Timothy, please evaluate your transpose to a column-major matrix.

a.transpose().eval() doesn't work because it evals to a row-major matrix.
You need to do something like:
MatrixType a_transpose = a.transpose();
Indeed, unless you explicitly require otherwise, any MatrixType has column-major storage.



On Fri, Oct 10, 2008 at 5:22 PM, Timothy Hunter <tjhunter@xxxxxxxxxxxx>wrote:

At one point, I was trying to figure out what was wrong in the
.solve() method of the cholesky decomposition and I decomposed the
solveTriangular(diag inverse(solveTriangular())) to check each step.
It had an assert failure about the row major bit.
Let me check if I can get this code back.

On Fri, Oct 10, 2008 at 7:35 AM, Benoît Jacob <jacob@xxxxxxxxxxxxxxx>
> The spamming continues...
> The LU::solve() actually does evaluate its argument b into a matrix c and
> only c from that point on... in particular it calls solveTriangular on c
> only... so I really don't understand why it has trouble.
> Moreover I added this test in test/lu.cpp:
>  // unit-test for a bad bug that escaped us for a long time:
>  // make sure that LU::solve() works well when passed an expression
>  m4 = (m1*m2).transpose();
>  lu.solve(m4.transpose(), &m2);
>  VERIFY_IS_APPROX(m4.transpose(), m1*m2);
> This test succeeds ! This is strange as it seems the same as Timothy's
> which makes LU::solve() fail.
> I'm puzzled!
> Benoit
> On Friday 10 October 2008 15:20:21 Benoît Jacob wrote:
>> Hm I read too fast: indeed we also have the same problem in LU.... which
>> my baby i.e. my fault.
>> Still not decided whether the right approach is to evaluate or to
>> nest-by-value...
>> Cheers,
>> Benoit
>> On Friday 10 October 2008 15:04:01 Benoît Jacob wrote:
>> > On Friday 10 October 2008 14:58:18 Benoît Jacob wrote:
>> > > Did you forget a nestByValue() ?
>> >
>> > Replying to myself: hm, rather a .eval() of course as a solver
>> > uses each coeff many times.
>> >
>> > By the way: the solve() API in Cholesky is not homogeneous with the
>> > solvers, it returns by value while the other solvers have a C-style
>> > Can you harmonize that by the way?
>> >
>> > Cheers,
>> > Benoit
>> >
>> > ---
>> ---
> ---

This message was sent using IMP, the Internet Messaging Program.


Mail converted by MHonArc 2.6.19+