Re: [eigen] LDLt and LLt fixes |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] LDLt and LLt fixes*From*: Bill Greene <w.h.greene@xxxxxxxxx>*Date*: Mon, 30 Mar 2009 21:47:11 -0400*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:in-reply-to:references :date:message-id:subject:from:to:content-type; bh=RO2uYqZHminP3pxbLWIhlAy62/YTArmrcQvkd6qFGcs=; b=o5cneZ3QSd2znK984vGjXN32+5vPZtZt4PdLrTYUCe/DTS9Eqn2Ww0xkvg9z7NmvED wfmS0mmuPlQAYIPcte0tMYsoWD6+mSjsY7nYM2Q4h1t4BmYRXRFSB9ThWSteOhG4NPC8 hvbemmzSP57FYMBnU9c4NdFoG12Mn8dSlmUDE=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type; b=uv6u8V/oqstNyDAevHmch2CY+Ol//R1hgQemW27dQs5vcYFyQvddh6zo6aMnikrLm5 98+PeeX2JPK4N/Lcig43+/jhIciYQBlmbwyeYJssvT8aym6RiqzrnIiY7n3yh1f3i8dh aKDiM2T+aWOdAQHHAVaT+p0NluyOPRWw6bgaU=

Benoit,

I think I accidentally sent you a test case that actually *requires* pivoting because the (0,0) term in the

matrix is zero. So please ignore it if you like.

Bill

On Mon, Mar 30, 2009 at 9:12 PM, Bill Greene <w.h.greene@xxxxxxxxx> wrote:

Benoit,

Wow, thanks for the quick work!

I just downloaded the latest code from svn and tried a quick test. Surprisingly, it didn't work. After some

debugging and investigation, I think this may be an unusual case. My other test cases work correctly--

exactly as I expect. Here is the code for my failed case

void bathe2DOF()

{

const int n = 2;

Eigen::Matrix2d K, M;

K << 6., -2., -2., 4.;

M << 2., 0., 0., 1.;

const double omega2 = 3.0;

Eigen::Matrix2d A = K - omega2*M;

Eigen::LDLT<Eigen::Matrix2d> factA(A);

std::cout << "Is invertible: " << factA.isInvertible() << std::endl;

Eigen::Vector2d b, x;

b << 0.0, 10.0;

factA.solve(b, &x);

prtMat(x, "x");

}

The correct value for x, by the way, is [-5.0, 0].

The lapack routines dsytrf and dsytrs compute this solution. I promised to provide more details of the lapack

solution-- and I intend to-- but I'm somewhat puzzled by the factored matrix returned by dsytrf. Even though

the solution from dsytrs is correct, the factored matrix is not what I expected-- or perhaps its just that the

documentation is so cryptic I don't know what I'm seeing. At any rate, I'll provide more lapack details as

soon as I understand them.

BillOn Mon, Mar 30, 2009 at 5:55 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

Hi,

the thread with Bill inspired me... so now (trunk revision 947097),

LDLt supports the negative semidefinite case.

There's a couple of things to discuss.

First, the floating-point comparisons in both LLt and LDLt had issues. LLt did

if(x < eps)

and LDLt did

if(Djj <= 0)

and replacing these by proper fuzzy comparisons (experimenting with a

few different ones) resulted in a big improvement in the accuracy of

some methods, like LDLt::rank() which used to be untested and had

typically 3% imprecision, and now is exact in 95% of cases (However

for an exact rank computation, it's still not nearly as reliable as

the full-pivoting LU which can be explained by the different between

bi-directional pivoting and the diagonal pivoting that's all one can

do with Cholesky).

The other issue is that LLt and LDLt have claims in their dox that

they use only a triangular half of the matrix and don't need the other

half to be initialized, well that's not currently the case, one can

check that by uncommenting a line I added to the cholesky.cpp

unit-test. Help fixing this welcome...

Cheers,

Benoit

**Follow-Ups**:**Re: [eigen] LDLt and LLt fixes***From:*Benoit Jacob

**References**:**[eigen] LDLt and LLt fixes***From:*Benoit Jacob

**Re: [eigen] LDLt and LLt fixes***From:*Bill Greene

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] LDLt and LLt fixes** - Next by Date:
**Re: [eigen] LDLt and LLt fixes** - Previous by thread:
**Re: [eigen] LDLt and LLt fixes** - Next by thread:
**Re: [eigen] LDLt and LLt fixes**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |