Re: [eigen] [RFC] Full pivoting LDLT

```On Tue, 3 Feb 2009, Keir Mierle wrote:

```
```http://codereview.appspot.com/14042/diff/1/4

There are a couple of issues.

1. Pivoting is not optional. Any objections?
2. If the matrix is only semidefinite, then matrixL() gives wrong results
(1's in diagonal where they shouldn't be).

Minus objections, I'll commit this.
```
```
```
I don't understand issue 2. I thought the L in the LDL^T decomposition has per definition 1s on the diagonal?
```
```
It looks like you assume that the matrix is positive semi-definite. However, I believe that the LDL^T decomposition is also useful for symmetric indefinite matrices. Solving Ax = b with A symmetric is twice as fast with LDL^T than with LU (if I remember correctly).
```
```
Why do you store P in two places: m_p and m_transpositions? Is using m_transpositions faster than m_p?
```
```
I was a bit surprised that you use both the triangles in m_matrix during the computation, since one half is just a rescaling of the other half. So I tried using only the lower triangle, but my timings (using g++ 4.1) are inconclusive on whether this is actually faster. I attach the diff for your entertainment, but I don't have any good reason for it.
```
Cheers,
Jitse```
```Index: Eigen/src/Cholesky/LDLT.h
===================================================================
--- Eigen/src/Cholesky/LDLT.h	(revision 921039)
+++ Eigen/src/Cholesky/LDLT.h	(working copy)
@@ -171,13 +171,14 @@
}

if (j == 0) {
-      m_matrix.row(0) = m_matrix.row(0).conjugate();
-      m_matrix.col(0).end(size-1) = m_matrix.row(0).end(size-1) / m_matrix.coeff(0,0);
+      m_matrix.col(0).end(size-1) /= m_matrix.coeff(0,0);
continue;
}

+    _temporary.start(j) = ( m_matrix.row(j).start(j).cwise()
+			    * m_matrix.diagonal().start(j).transpose() ).conjugate();
RealScalar Djj = ei_real(m_matrix.coeff(j,j) - (m_matrix.row(j).start(j)
-                                                  * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
+                                                  * _temporary.start(j)).coeff(0,0));
m_matrix.coeffRef(j,j) = Djj;

// Finish early if the matrix is not full rank or is indefinite.  This
@@ -194,12 +195,10 @@
int endSize = size - j - 1;
if (endSize > 0) {
_temporary.end(endSize) = ( m_matrix.block(j+1,0, endSize, j)
-                                * m_matrix.col(j).start(j).conjugate() ).lazy();
+                                * _temporary.start(j) ).lazy();

-      m_matrix.row(j).end(endSize) = m_matrix.row(j).end(endSize).conjugate()
-                                   - _temporary.end(endSize).transpose();
-
-      m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / Djj;
+      m_matrix.col(j).end(endSize) = ( m_matrix.col(j).end(endSize)
+				       - _temporary.end(endSize)    ) / Djj;
}
}

```

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