| Re: [eigen] [RFC] Full pivoting LDLT |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
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. 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.
Cheers,
Jitse