Re: [eigen] [RFC] Full pivoting LDLT |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] [RFC] Full pivoting LDLT*From*: Keir Mierle <mierle@xxxxxxxxx>*Date*: Thu, 5 Feb 2009 10:45:46 -0800*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 :content-transfer-encoding; bh=PZOdflX4GH66eOyQiMrVYbVPQSKbCEBixhBY4XeYBFY=; b=I/2h8uJr+STeteXjQg+mUy6KvM08nNd5UmHHyVbDy4Bw86Z6dwgvNT8VfCEH6nnh3y O+I4/JtEJxY4PQS8I0r6dq6Cb4W7F2xsub3kDy11uGtyZDAditvpReyqcAO7UYr+7/u1 gE3Cc4mQAYuYjQm3odO1JYL/IrPcpPoZTUiEM=*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:content-transfer-encoding; b=GxhENEOL/9RhxrNMKd9DX8TsP3bvgmO83/r1EgMaQpWqz6RymzjP2oC5ml/zROUVsP HEuY8av0cGnl52koEpvXnJqlNpxVDS+vLmUm+CWa3QeeITNy9Gc9U77PopxegsAlpufo R+mSJkFuH/0GiDifjhj8rvcTji5Tv021Tw+GM=

On Thu, Feb 5, 2009 at 10:31 AM, Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx> wrote: > On Wed, 4 Feb 2009, Keir Mierle wrote: > >> If A is symmetric indefinite, then a very different algorithm is >> required because D is block-diagonal as 1x1 and 2x2 blocks. > > I don't understand why D would have 2x2 blocks. Could you please explain > this, or tell me where you read this? > >> I left that in for the non-solve-in-place case. If you are solving in >> place, you need to apply the transpositions via m_transpositions. If >> you are solving into a new matrix, then m_p should be used instead >> (see LU.h's use of m_p in solve()). I haven't implemented solve() >> using m_p yet though (feel free!). m_p is a direct mapping from one >> permutation to another; transpositions is the series of swaps to apply >> to get to the other permutation. > > Fair enough, thank you for your explanation. > > I'm not sure transpositions are the fastest way to do this. Neither am I :) I did this because I didn't want to make a sizeof(b.rows()) copy when solving in place. It is probably faster to do that in some cases though. I'm open to measurements. It's an easy change to make. Actually, this brings up something else I was thinking of: We should make efficient permutation matrix class, which has a specialization for matrix multiply which permutes the matrix. This way we could return a meaningful matrixP() from the pivoting decompositions, and it would remove code duplication (both LU and LDLT do pivoting now, so will QR when I get a chance). > For instance, > suppose you have a 3x3 matrix, and your pivots are (1,1), (2,2) and (0,0). > Thus m_transpositions is (1,2,2) so you permute the vector b on the > right-hand side by doing two swaps, thus six reads and writes. However, it's > possible to do it in four: > tmp = b(0); b(0) = b(1); b(1) = b(2); b(2) = tmp; > > On the other hand, the triangular solve is O(n^2) while the permutation is > O(n) so it probably doesn't make sense to optimize the permutation. Again, for the in-place solve I didn't want to allocate a O(n) temporary. However, it is probably worthwhile. We should measure this. Keir > > Cheers, > Jitse > > >

**References**:**[eigen] [RFC] Full pivoting LDLT***From:*Keir Mierle

**Re: [eigen] [RFC] Full pivoting LDLT***From:*Jitse Niesen

**Re: [eigen] [RFC] Full pivoting LDLT***From:*Keir Mierle

**Re: [eigen] [RFC] Full pivoting LDLT***From:*Jitse Niesen

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] [RFC] Full pivoting LDLT** - Next by Date:
**[eigen] Re: [PATCH] Reverse expression (trying again)** - Previous by thread:
**Re: [eigen] [RFC] Full pivoting LDLT** - Next by thread:
**[eigen] Machine precision<> too coarse**

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