add a minimum degree ordering routine based on CSparse (LGPL) and a new built-in sparse cholesky decomposition

7 files changed, 1174 insertions, 85 deletions

diff --git a/unsupported/Eigen/SparseExtra b/unsupported/Eigen/SparseExtra
index 8b4afa3a8..45b154339 100644
--- a/unsupported/Eigen/SparseExtra
+++ b/unsupported/Eigen/SparseExtra
@@ -54,8 +54,12 @@ enum {
 #include "src/misc/Solve.h"
 
 #include "src/SparseExtra/RandomSetter.h"
+#include "src/SparseExtra/Solve.h"
+#include "src/SparseExtra/Amd.h"
+#include "src/SparseExtra/SimplicialCholesky.h"
+
 #include "src/SparseExtra/SparseLLT.h"
-#include "src/SparseExtra/SparseLDLT.h"
+#include "src/SparseExtra/SparseLDLTLegacy.h"
 #include "src/SparseExtra/SparseLU.h"
 
 } // namespace Eigen
diff --git a/unsupported/Eigen/src/SparseExtra/Amd.h b/unsupported/Eigen/src/SparseExtra/Amd.h
new file mode 100644
index 000000000..a716305fd
--- /dev/null
+++ b/unsupported/Eigen/src/SparseExtra/Amd.h
@@ -0,0 +1,487 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+/*
+
+NOTE: this routine has been adapted from the CSparse library:
+
+Copyright (c) 2006, Timothy A. Davis.
+http://www.cise.ufl.edu/research/sparse/CSparse
+
+CSparse is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+CSparse is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with this Module; if not, write to the Free Software
+Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+
+*/
+
+#ifndef EIGEN_SPARSE_AMD_H
+#define EIGEN_SPARSE_AMD_H
+
+namespace internal {
+  
+
+#define CS_FLIP(i) (-(i)-2)
+#define CS_UNFLIP(i) (((i) < 0) ? CS_FLIP(i) : (i))
+#define CS_MARKED(w,j) (w[j] < 0)
+#define CS_MARK(w,j) { w[j] = CS_FLIP (w[j]); }
+
+/* clear w */
+template<typename Index>
+static int cs_wclear (Index mark, Index lemax, Index *w, Index n)
+{
+  Index k;
+  if(mark < 2 || (mark + lemax < 0))
+  {
+    for(k = 0; k < n; k++)
+      if(w[k] != 0)
+        w[k] = 1;
+    mark = 2;
+  }
+  return (mark);     /* at this point, w[0..n-1] < mark holds */
+}
+
+/* depth-first search and postorder of a tree rooted at node j */
+template<typename Index>
+Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Index *stack)
+{
+  int i, p, top = 0;
+  if(!head || !next || !post || !stack) return (-1);    /* check inputs */
+  stack[0] = j;                 /* place j on the stack */
+  while (top >= 0)                /* while (stack is not empty) */
+  {
+    p = stack[top];           /* p = top of stack */
+    i = head[p];              /* i = youngest child of p */
+    if(i == -1)
+    {
+      top--;                 /* p has no unordered children left */
+      post[k++] = p;        /* node p is the kth postordered node */
+    }
+    else
+    {
+      head[p] = next[i];   /* remove i from children of p */
+      stack[++top] = i;     /* start dfs on child node i */
+    }
+  }
+  return k;
+}
+
+/** keeps off-diagonal entries; drops diagonal entries */
+template<typename Index, typename Scalar>
+struct keep_diag {
+  inline bool operator() (const Index& row, const Index& col, const Scalar&) const
+  {
+    return row!=col;
+  }
+};
+
+/** p = amd(A+A') if symmetric is true, or amd(A'A) otherwise */
+template<typename Scalar, int Options, typename Index>
+int *minimum_degree_ordering(int order, const SparseMatrix<Scalar,Options,Index>& A)  /* order 0:natural, 1:Chol, 2:LU, 3:QR */
+{
+  typedef SparseMatrix<Scalar,ColMajor,Index> CCS;
+  
+  int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
+      k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
+      ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t;
+  unsigned int h;
+  /* --- Construct matrix C ----------------------------------------------- */
+  if(order <= 0 || order > 3) return (NULL); /* check */
+  
+  Index m = A.rows();
+  Index n = A.cols();
+  dense = std::max<Index> (16, 10 * sqrt ((double) n));   /* find dense threshold */
+  dense = std::min<Index> (n-2, dense);
+  CCS C;
+  if(order == 1 && n == m) // Cholesky
+  {
+    C = A + SparseMatrix<Scalar,Options,Index>(A.adjoint());
+  }
+  else if(order == 2)  // LU
+  {
+    CCS AT = A.adjoint();
+    // drop dense columns from AT
+    Index* ATp = AT._outerIndexPtr();
+    Index* ATi = AT._innerIndexPtr();
+    Index p2;
+    Index j;
+    for(p2 = 0, j = 0; j < m; j++)
+    {
+      Index p = ATp[j];                         // column j of AT starts here
+      ATp[j] = p2;                              // new column j starts here
+      if(ATp[j+1] - p > dense) continue;        // skip dense col j
+      for(; p < ATp[j+1]; p++)
+        ATi[p2++] = ATi[p];
+    }
+    ATp[m] = p2;                                // finalize AT
+    // TODO this could be implemented using a sparse filter expression
+    // TODO do a cheap selfadjoint rank update
+    C = AT * AT.adjoint();                    // C=A'*A with no dense rows
+  }
+  else  // QR
+  {
+    C = A.adjoint() * A;
+  }
+    
+  C.prune(keep_diag<Index,Scalar>());
+  
+  
+  Index cnz = A.nonZeros();
+  Index* P = new Index[n+1];     /* allocate result */
+  t = cnz + cnz/5 + 2*n;                 /* add elbow room to C */
+  C.resizeNonZeros(t);
+  
+  Index* W       = new Index[8*(n+1)]; /* get workspace */
+  Index* len     = W;
+  Index* nv      = W +   (n+1);
+  Index* next    = W + 2*(n+1);
+  Index* head    = W + 3*(n+1);
+  Index* elen    = W + 4*(n+1);
+  Index* degree  = W + 5*(n+1);
+  Index* w       = W + 6*(n+1);
+  Index* hhead   = W + 7*(n+1);
+  Index* last    = P;                              /* use P as workspace for last */
+  
+  /* --- Initialize quotient graph ---------------------------------------- */
+  Index* Cp = C._outerIndexPtr();
+  Index* Ci = C._innerIndexPtr();
+  for(k = 0; k < n; k++)
+    len[k] = Cp[k+1] - Cp[k];
+  len[n] = 0;
+  nzmax = t;
+  
+  for(i = 0; i <= n; i++)
+  {
+    head[i]   = -1;                     // degree list i is empty
+    last[i]   = -1;
+    next[i]   = -1;
+    hhead[i]  = -1;                     // hash list i is empty 
+    nv[i]     = 1;                      // node i is just one node
+    w[i]      = 1;                      // node i is alive
+    elen[i]   = 0;                      // Ek of node i is empty
+    degree[i] = len[i];                 // degree of node i
+  }
+  mark = cs_wclear (0, 0, w, n);         /* clear w */
+  elen[n] = -2;                         /* n is a dead element */
+  Cp[n] = -1;                           /* n is a root of assembly tree */
+  w[n] = 0;                             /* n is a dead element */
+  
+  /* --- Initialize degree lists ------------------------------------------ */
+  for(i = 0; i < n; i++)
+  {
+    d = degree[i];
+    if(d == 0)                         /* node i is empty */
+    {
+      elen[i] = -2;                 /* element i is dead */
+      nel++;
+      Cp[i] = -1;                   /* i is a root of assembly tree */
+      w[i] = 0;
+    }
+    else if(d > dense)                 /* node i is dense */
+    {
+      nv[i] = 0;                    /* absorb i into element n */
+      elen[i] = -1;                 /* node i is dead */
+      nel++;
+      Cp[i] = CS_FLIP (n);
+      nv[n]++;
+    }
+    else
+    {
+      if(head[d] != -1) last[head[d]] = i;
+      next[i] = head[d];           /* put node i in degree list d */
+      head[d] = i;
+    }
+  }
+  
+  while (nel < n)                         /* while (selecting pivots) do */
+  {
+    /* --- Select node of minimum approximate degree -------------------- */
+    for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++);
+    if(next[k] != -1) last[next[k]] = -1;
+    head[mindeg] = next[k];          /* remove k from degree list */
+    elenk = elen[k];                  /* elenk = |Ek| */
+    nvk = nv[k];                      /* # of nodes k represents */
+    nel += nvk;                        /* nv[k] nodes of A eliminated */
+    
+    /* --- Garbage collection ------------------------------------------- */
+    if(elenk > 0 && cnz + mindeg >= nzmax)
+    {
+      for(j = 0; j < n; j++)
+      {
+        if((p = Cp[j]) >= 0)      /* j is a live node or element */
+        {
+          Cp[j] = Ci[p];          /* save first entry of object */
+          Ci[p] = CS_FLIP (j);    /* first entry is now CS_FLIP(j) */
+        }
+      }
+      for(q = 0, p = 0; p < cnz; ) /* scan all of memory */
+      {
+        if((j = CS_FLIP (Ci[p++])) >= 0)  /* found object j */
+        {
+          Ci[q] = Cp[j];       /* restore first entry of object */
+          Cp[j] = q++;          /* new pointer to object j */
+          for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++];
+        }
+      }
+      cnz = q;                       /* Ci[cnz...nzmax-1] now free */
+    }
+    
+    /* --- Construct new element ---------------------------------------- */
+    dk = 0;
+    nv[k] = -nvk;                     /* flag k as in Lk */
+    p = Cp[k];
+    pk1 = (elenk == 0) ? p : cnz;      /* do in place if elen[k] == 0 */
+    pk2 = pk1;
+    for(k1 = 1; k1 <= elenk + 1; k1++)
+    {
+      if(k1 > elenk)
+      {
+        e = k;                     /* search the nodes in k */
+        pj = p;                    /* list of nodes starts at Ci[pj]*/
+        ln = len[k] - elenk;      /* length of list of nodes in k */
+      }
+      else
+      {
+        e = Ci[p++];              /* search the nodes in e */
+        pj = Cp[e];
+        ln = len[e];              /* length of list of nodes in e */
+      }
+      for(k2 = 1; k2 <= ln; k2++)
+      {
+        i = Ci[pj++];
+        if((nvi = nv[i]) <= 0) continue; /* node i dead, or seen */
+        dk += nvi;                 /* degree[Lk] += size of node i */
+        nv[i] = -nvi;             /* negate nv[i] to denote i in Lk*/
+        Ci[pk2++] = i;            /* place i in Lk */
+        if(next[i] != -1) last[next[i]] = last[i];
+        if(last[i] != -1)         /* remove i from degree list */
+        {
+          next[last[i]] = next[i];
+        }
+        else
+        {
+          head[degree[i]] = next[i];
+        }
+      }
+      if(e != k)
+      {
+        Cp[e] = CS_FLIP (k);      /* absorb e into k */
+        w[e] = 0;                 /* e is now a dead element */
+      }
+    }
+    if(elenk != 0) cnz = pk2;         /* Ci[cnz...nzmax] is free */
+    degree[k] = dk;                   /* external degree of k - |Lk\i| */
+    Cp[k] = pk1;                      /* element k is in Ci[pk1..pk2-1] */
+    len[k] = pk2 - pk1;
+    elen[k] = -2;                     /* k is now an element */
+    
+    /* --- Find set differences ----------------------------------------- */
+    mark = cs_wclear (mark, lemax, w, n);  /* clear w if necessary */
+    for(pk = pk1; pk < pk2; pk++)    /* scan 1: find |Le\Lk| */
+    {
+      i = Ci[pk];
+      if((eln = elen[i]) <= 0) continue;/* skip if elen[i] empty */
+      nvi = -nv[i];                      /* nv[i] was negated */
+      wnvi = mark - nvi;
+      for(p = Cp[i]; p <= Cp[i] + eln - 1; p++)  /* scan Ei */
+      {
+        e = Ci[p];
+        if(w[e] >= mark)
+        {
+          w[e] -= nvi;          /* decrement |Le\Lk| */
+        }
+        else if(w[e] != 0)        /* ensure e is a live element */
+        {
+          w[e] = degree[e] + wnvi; /* 1st time e seen in scan 1 */
+        }
+      }
+    }
+    
+    /* --- Degree update ------------------------------------------------ */
+    for(pk = pk1; pk < pk2; pk++)    /* scan2: degree update */
+    {
+      i = Ci[pk];                   /* consider node i in Lk */
+      p1 = Cp[i];
+      p2 = p1 + elen[i] - 1;
+      pn = p1;
+      for(h = 0, d = 0, p = p1; p <= p2; p++)    /* scan Ei */
+      {
+        e = Ci[p];
+        if(w[e] != 0)             /* e is an unabsorbed element */
+        {
+          dext = w[e] - mark;   /* dext = |Le\Lk| */
+          if(dext > 0)
+          {
+            d += dext;         /* sum up the set differences */
+            Ci[pn++] = e;     /* keep e in Ei */
+            h += e;            /* compute the hash of node i */
+          }
+          else
+          {
+            Cp[e] = CS_FLIP (k);  /* aggressive absorb. e->k */
+            w[e] = 0;             /* e is a dead element */
+          }
+        }
+      }
+      elen[i] = pn - p1 + 1;        /* elen[i] = |Ei| */
+      p3 = pn;
+      p4 = p1 + len[i];
+      for(p = p2 + 1; p < p4; p++) /* prune edges in Ai */
+      {
+        j = Ci[p];
+        if((nvj = nv[j]) <= 0) continue; /* node j dead or in Lk */
+        d += nvj;                  /* degree(i) += |j| */
+        Ci[pn++] = j;             /* place j in node list of i */
+        h += j;                    /* compute hash for node i */
+      }
+      if(d == 0)                     /* check for mass elimination */
+      {
+        Cp[i] = CS_FLIP (k);      /* absorb i into k */
+        nvi = -nv[i];
+        dk -= nvi;                 /* |Lk| -= |i| */
+        nvk += nvi;                /* |k| += nv[i] */
+        nel += nvi;
+        nv[i] = 0;
+        elen[i] = -1;             /* node i is dead */
+      }
+      else
+      {
+        degree[i] = std::min<Index> (degree[i], d);   /* update degree(i) */
+        Ci[pn] = Ci[p3];         /* move first node to end */
+        Ci[p3] = Ci[p1];         /* move 1st el. to end of Ei */
+        Ci[p1] = k;               /* add k as 1st element in of Ei */
+        len[i] = pn - p1 + 1;     /* new len of adj. list of node i */
+        h %= n;                    /* finalize hash of i */
+        next[i] = hhead[h];      /* place i in hash bucket */
+        hhead[h] = i;
+        last[i] = h;              /* save hash of i in last[i] */
+      }
+    }                                   /* scan2 is done */
+    degree[k] = dk;                   /* finalize |Lk| */
+    lemax = std::max<Index>(lemax, dk);
+    mark = cs_wclear (mark+lemax, lemax, w, n);    /* clear w */
+    
+    /* --- Supernode detection ------------------------------------------ */
+    for(pk = pk1; pk < pk2; pk++)
+    {
+      i = Ci[pk];
+      if(nv[i] >= 0) continue;         /* skip if i is dead */
+      h = last[i];                      /* scan hash bucket of node i */
+      i = hhead[h];
+      hhead[h] = -1;                    /* hash bucket will be empty */
+      for(; i != -1 && next[i] != -1; i = next[i], mark++)
+      {
+        ln = len[i];
+        eln = elen[i];
+        for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark;
+        jlast = i;
+        for(j = next[i]; j != -1; ) /* compare i with all j */
+        {
+          ok = (len[j] == ln) && (elen[j] == eln);
+          for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++)
+          {
+            if(w[Ci[p]] != mark) ok = 0;    /* compare i and j*/
+          }
+          if(ok)                     /* i and j are identical */
+          {
+            Cp[j] = CS_FLIP (i);  /* absorb j into i */
+            nv[i] += nv[j];
+            nv[j] = 0;
+            elen[j] = -1;         /* node j is dead */
+            j = next[j];          /* delete j from hash bucket */
+            next[jlast] = j;
+          }
+          else
+          {
+            jlast = j;             /* j and i are different */
+            j = next[j];
+          }
+        }
+      }
+    }
+    
+    /* --- Finalize new element------------------------------------------ */
+    for(p = pk1, pk = pk1; pk < pk2; pk++)   /* finalize Lk */
+    {
+      i = Ci[pk];
+      if((nvi = -nv[i]) <= 0) continue;/* skip if i is dead */
+      nv[i] = nvi;                      /* restore nv[i] */
+      d = degree[i] + dk - nvi;         /* compute external degree(i) */
+      d = std::min<Index> (d, n - nel - nvi);
+      if(head[d] != -1) last[head[d]] = i;
+      next[i] = head[d];               /* put i back in degree list */
+      last[i] = -1;
+      head[d] = i;
+      mindeg = std::min<Index> (mindeg, d);       /* find new minimum degree */
+      degree[i] = d;
+      Ci[p++] = i;                      /* place i in Lk */
+    }
+    nv[k] = nvk;                      /* # nodes absorbed into k */
+    if((len[k] = p-pk1) == 0)         /* length of adj list of element k*/
+    {
+      Cp[k] = -1;                   /* k is a root of the tree */
+      w[k] = 0;                     /* k is now a dead element */
+    }
+    if(elenk != 0) cnz = p;           /* free unused space in Lk */
+  }
+  
+  /* --- Postordering ----------------------------------------------------- */
+  for(i = 0; i < n; i++) Cp[i] = CS_FLIP (Cp[i]);/* fix assembly tree */
+  for(j = 0; j <= n; j++) head[j] = -1;
+  for(j = n; j >= 0; j--)              /* place unordered nodes in lists */
+  {
+    if(nv[j] > 0) continue;          /* skip if j is an element */
+    next[j] = head[Cp[j]];          /* place j in list of its parent */
+    head[Cp[j]] = j;
+  }
+  for(e = n; e >= 0; e--)              /* place elements in lists */
+  {
+    if(nv[e] <= 0) continue;         /* skip unless e is an element */
+    if(Cp[e] != -1)
+    {
+      next[e] = head[Cp[e]];      /* place e in list of its parent */
+      head[Cp[e]] = e;
+    }
+  }
+  for(k = 0, i = 0; i <= n; i++)       /* postorder the assembly tree */
+  {
+    if(Cp[i] == -1) k = cs_tdfs (i, k, head, next, P, w);
+  }
+
+  delete[] W;
+  return P;
+}
+
+} // namespace internal
+
+#endif // EIGEN_SPARSE_AMD_H
diff --git a/unsupported/Eigen/src/SparseExtra/CholmodSupport.h b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
index b0f1c8e42..14f18d30f 100644
--- a/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
+++ b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
@@ -27,56 +27,6 @@
 
 namespace internal {
 
-template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval_base;
-template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval;
-  
-template<typename DecompositionType, typename Rhs>
-struct traits<sparse_solve_retval_base<DecompositionType, Rhs> >
-{
-  typedef typename DecompositionType::MatrixType MatrixType;
-  typedef SparseMatrix<typename Rhs::Scalar, Rhs::Options, typename Rhs::Index> ReturnType;
-};
-
-template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval_base
- : public ReturnByValue<sparse_solve_retval_base<_DecompositionType, Rhs> >
-{
-  typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned;
-  typedef _DecompositionType DecompositionType;
-  typedef ReturnByValue<sparse_solve_retval_base> Base;
-  typedef typename Base::Index Index;
-
-  sparse_solve_retval_base(const DecompositionType& dec, const Rhs& rhs)
-    : m_dec(dec), m_rhs(rhs)
-  {}
-
-  inline Index rows() const { return m_dec.cols(); }
-  inline Index cols() const { return m_rhs.cols(); }
-  inline const DecompositionType& dec() const { return m_dec; }
-  inline const RhsNestedCleaned& rhs() const { return m_rhs; }
-
-  template<typename Dest> inline void evalTo(Dest& dst) const
-  {
-    static_cast<const sparse_solve_retval<DecompositionType,Rhs>*>(this)->evalTo(dst);
-  }
-
-  protected:
-    const DecompositionType& m_dec;
-    const typename Rhs::Nested m_rhs;
-};
-
-#define EIGEN_MAKE_SPARSE_SOLVE_HELPERS(DecompositionType,Rhs) \
-  typedef typename DecompositionType::MatrixType MatrixType; \
-  typedef typename MatrixType::Scalar Scalar; \
-  typedef typename MatrixType::RealScalar RealScalar; \
-  typedef typename MatrixType::Index Index; \
-  typedef Eigen::internal::sparse_solve_retval_base<DecompositionType,Rhs> Base; \
-  using Base::dec; \
-  using Base::rhs; \
-  using Base::rows; \
-  using Base::cols; \
-  sparse_solve_retval(const DecompositionType& dec, const Rhs& rhs) \
-    : Base(dec, rhs) {}
-
 template<typename Scalar, typename CholmodType>
 void cholmod_configure_matrix(CholmodType& mat)
 {
@@ -243,8 +193,8 @@ class CholmodDecomposition
       cholmod_finish(&m_cholmod);
     }
     
-    int cols() const { return m_cholmodFactor->n; }
-    int rows() const { return m_cholmodFactor->n; }
+    inline Index cols() const { return m_cholmodFactor->n; }
+    inline Index rows() const { return m_cholmodFactor->n; }
     
     void setMode(CholmodMode mode)
     {
diff --git a/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h b/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h
new file mode 100644
index 000000000..6b22d1502
--- /dev/null
+++ b/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h
@@ -0,0 +1,457 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+/*
+
+NOTE: the _symbolic, and _numeric functions has been adapted from
+      the LDL library:
+
+LDL Copyright (c) 2005 by Timothy A. Davis.  All Rights Reserved.
+
+LDL License:
+
+    Your use or distribution of LDL or any modified version of
+    LDL implies that you agree to this License.
+
+    This library is free software; you can redistribute it and/or
+    modify it under the terms of the GNU Lesser General Public
+    License as published by the Free Software Foundation; either
+    version 2.1 of the License, or (at your option) any later version.
+
+    This library is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    Lesser General Public License for more details.
+
+    You should have received a copy of the GNU Lesser General Public
+    License along with this library; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301
+    USA
+
+    Permission is hereby granted to use or copy this program under the
+    terms of the GNU LGPL, provided that the Copyright, this License,
+    and the Availability of the original version is retained on all copies.
+    User documentation of any code that uses this code or any modified
+    version of this code must cite the Copyright, this License, the
+    Availability note, and "Used by permission." Permission to modify
+    the code and to distribute modified code is granted, provided the
+    Copyright, this License, and the Availability note are retained,
+    and a notice that the code was modified is included.
+ */
+
+#ifndef EIGEN_SIMPLICIAL_CHOLESKY_H
+#define EIGEN_SIMPLICIAL_CHOLESKY_H
+
+enum SimplicialCholeskyMode {
+  SimplicialCholeskyLLt,
+  SimplicialCholeskyLDLt
+};
+
+/** \brief A direct sparse Cholesky factorization
+  *
+  * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization.
+  * The sparse matrix A must be selfajoint and positive definite. The vectors or matrices
+  * X and B can be either dense or sparse.
+  *
+  * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
+  * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
+  *               or Upper. Default is Lower.
+  *
+  */
+template<typename _MatrixType, int _UpLo = Lower>
+class SimplicialCholesky
+{
+  public:
+    typedef _MatrixType MatrixType;
+    enum { UpLo = _UpLo };
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    typedef typename MatrixType::Index Index;
+    typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
+    typedef Matrix<Scalar,MatrixType::ColsAtCompileTime,1> VectorType;
+
+  public:
+
+    SimplicialCholesky()
+      : m_info(Success), m_isInitialized(false), m_LDLt(true)
+    {}
+
+    SimplicialCholesky(const MatrixType& matrix)
+      : m_info(Success), m_isInitialized(false), m_LDLt(true)
+    {
+      compute(matrix);
+    }
+
+    ~SimplicialCholesky()
+    {
+    }
+    
+    inline Index cols() const { return m_matrix.cols(); }
+    inline Index rows() const { return m_matrix.rows(); }
+  
+    SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
+    {
+      switch(mode)
+      {
+        case SimplicialCholeskyLLt:
+          m_LDLt = false;
+          break;
+        case SimplicialCholeskyLDLt:
+          m_LDLt = true;
+          break;
+        default:
+          break;
+      }
+      
+      return *this;
+    }
+    
+    /** \brief Reports whether previous computation was successful.
+      *
+      * \returns \c Success if computation was succesful,
+      *          \c NumericalIssue if the matrix.appears to be negative.
+      */
+    ComputationInfo info() const
+    {
+      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
+      return m_info;
+    }
+
+    /** Computes the sparse Cholesky decomposition of \a matrix */
+    SimplicialCholesky& compute(const MatrixType& matrix)
+    {
+      analyzePattern(matrix);
+      factorize(matrix);
+      return *this;
+    }
+    
+    /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
+      *
+      * \sa compute()
+      */
+    template<typename Rhs>
+    inline const internal::solve_retval<SimplicialCholesky, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      eigen_assert(m_isInitialized && "SimplicialCholesky is not initialized.");
+      eigen_assert(rows()==b.rows()
+                && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
+      return internal::solve_retval<SimplicialCholesky, Rhs>(*this, b.derived());
+    }
+    
+    /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
+      *
+      * \sa compute()
+      */
+//     template<typename Rhs>
+//     inline const internal::sparse_solve_retval<CholmodDecomposition, Rhs>
+//     solve(const SparseMatrixBase<Rhs>& b) const
+//     {
+//       eigen_assert(m_isInitialized && "SimplicialCholesky is not initialized.");
+//       eigen_assert(rows()==b.rows()
+//                 && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
+//       return internal::sparse_solve_retval<SimplicialCholesky, Rhs>(*this, b.derived());
+//     }
+    
+    /** Performs a symbolic decomposition on the sparcity of \a matrix.
+      *
+      * This function is particularly useful when solving for several problems having the same structure.
+      * 
+      * \sa factorize()
+      */
+    void analyzePattern(const MatrixType& a)
+    {
+      eigen_assert(a.rows()==a.cols());
+      const Index size = a.rows();
+      m_matrix.resize(size, size);
+      m_parent.resize(size);
+      m_nonZerosPerCol.resize(size);
+      
+      Index* tags = ei_aligned_stack_new(Index, size);
+      
+      // TODO allows to configure the permutation
+      const Index* P = internal::minimum_degree_ordering(1, a);
+      const Index* Pinv = 0;
+      if(P)
+      {
+        m_P.indices() = VectorXi::Map(P,size);
+        m_Pinv        = m_P.inverse();
+        Pinv          = m_Pinv.indices().data();
+      }
+      else
+      {
+        m_P.resize(0);
+        m_Pinv.resize(0);
+      }
+
+      for(Index k = 0; k < size; ++k)
+      {
+        /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
+        m_parent[k] = -1;             /* parent of k is not yet known */
+        tags[k] = k;                  /* mark node k as visited */
+        m_nonZerosPerCol[k] = 0;      /* count of nonzeros in column k of L */
+        Index kk = P ? P[k] : k;      /* kth original, or permuted, column */
+        for(typename MatrixType::InnerIterator it(a,kk); it; ++it)
+        {
+          /* A (i,k) is nonzero (original or permuted A) */
+          Index i = Pinv ? Pinv[it.index()] : it.index();
+          if(i < k)
+          {
+            /* follow path from i to root of etree, stop at flagged node */
+            for(; tags[i] != k; i = m_parent[i])
+            {
+              /* find parent of i if not yet determined */
+              if (m_parent[i] == -1)
+                m_parent[i] = k;
+              ++m_nonZerosPerCol[i];        /* L (k,i) is nonzero */
+              tags[i] = k;                  /* mark i as visited */
+            }
+          }
+        }
+      }
+      
+      // release worspace
+      ei_aligned_stack_delete(Index, tags, size);
+      
+      /* construct Lp index array from m_nonZerosPerCol column counts */
+      Index* Lp = m_matrix._outerIndexPtr();
+      Lp[0] = 0;
+      for(Index k = 0; k < size; ++k)
+        Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (m_LDLt ? 0 : 1);
+
+      m_matrix.resizeNonZeros(Lp[size]);
+      
+      m_isInitialized     = true;
+      m_info              = Success;
+      m_analysisIsOk      = true;
+      m_factorizationIsOk = false;
+    }
+    
+    /** Performs a numeric decomposition of \a matrix
+      *
+      * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
+      *
+      * \sa analyzePattern()
+      */
+    void factorize(const MatrixType& a)
+    {
+      eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
+      eigen_assert(a.rows()==a.cols());
+      const Index size = a.rows();
+      eigen_assert(m_parent.size()==size);
+      eigen_assert(m_nonZerosPerCol.size()==size);
+
+      const Index* Lp = m_matrix._outerIndexPtr();
+      Index* Li = m_matrix._innerIndexPtr();
+      Scalar* Lx = m_matrix._valuePtr();
+
+      Scalar* y       = ei_aligned_stack_new(Scalar, size);
+      Index* pattern  = ei_aligned_stack_new(Index, size);
+      Index* tags     = ei_aligned_stack_new(Index, size);
+      
+      Index* P    = 0;
+      Index* Pinv = 0;
+      
+      if(m_P.size()==size)
+      {
+        P     = m_P.indices().data();
+        Pinv  = m_Pinv.indices().data();
+      }
+      
+      bool ok = true;
+      
+      m_diag.resize(m_LDLt ? size : 0);
+      
+      for(Index k = 0; k < size; ++k)
+      {
+        /* compute nonzero pattern of kth row of L, in topological order */
+        y[k] = 0.0;                     /* Y(0:k) is now all zero */
+        Index top = size;               /* stack for pattern is empty */
+        tags[k] = k;                    /* mark node k as visited */
+        m_nonZerosPerCol[k] = 0;        /* count of nonzeros in column k of L */
+        Index kk = (P) ? (P[k]) : (k);  /* kth original, or permuted, column */
+        for(typename MatrixType::InnerIterator it(a,kk); it; ++it)
+        {
+          Index i = Pinv ? Pinv[it.index()] : it.index();
+          if(i <= k)
+          {
+            y[i] += internal::conj(it.value());            /* scatter A(i,k) into Y (sum duplicates) */
+            Index len;
+            for(len = 0; tags[i] != k; i = m_parent[i])
+            {
+              pattern[len++] = i;     /* L(k,i) is nonzero */
+              tags[i] = k;            /* mark i as visited */
+            }
+            while(len > 0)
+              pattern[--top] = pattern[--len];
+          }
+        }
+
+        /* compute numerical values kth row of L (a sparse triangular solve) */
+        Scalar d = y[k];                  // get D(k,k) and clear Y(k)
+        y[k] = 0.0;
+        for(; top < size; ++top)
+        {
+          if(1)
+          {
+            Index i = pattern[top];       /* pattern[top:n-1] is pattern of L(:,k) */
+            Scalar yi = y[i];             /* get and clear Y(i) */
+            y[i] = 0.0;
+            
+            /* the nonzero entry L(k,i) */
+            Scalar l_ki;
+            if(m_LDLt)
+              l_ki = yi / m_diag[i];       
+            else
+              yi = l_ki = yi / Lx[Lp[i]];
+            
+            Index p2 = Lp[i] + m_nonZerosPerCol[i];
+            Index p;
+            for(p = Lp[i] + (m_LDLt ? 0 : 1); p < p2; ++p)
+              y[Li[p]] -= internal::conj(Lx[p]) * yi;
+            d -= l_ki * internal::conj(yi);
+            Li[p] = k;                          /* store L(k,i) in column form of L */
+            Lx[p] = l_ki;
+            ++m_nonZerosPerCol[i];              /* increment count of nonzeros in col i */
+          }
+        }
+        if(m_LDLt)
+          m_diag[k] = d;
+        else
+        {
+          Index p = Lp[k]+m_nonZerosPerCol[k]++;
+          Li[p] = k ;                /* store L(k,k) = sqrt (d) in column k */
+          Lx[p] = internal::sqrt(d) ;
+        }
+        if(d == Scalar(0))
+        {
+          ok = false;                         /* failure, D(k,k) is zero */
+          break;
+        }
+      }
+
+      // release workspace
+      ei_aligned_stack_delete(Scalar, y, size);
+      ei_aligned_stack_delete(Index, pattern, size);
+      ei_aligned_stack_delete(Index, tags, size);
+      
+      m_info = ok ? Success : NumericalIssue;
+      m_factorizationIsOk = true;
+    }
+    
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
+    /** \internal */
+    template<typename Rhs,typename Dest>
+    void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
+    {
+      eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
+      eigen_assert(m_matrix.rows()==b.rows());
+      
+      if(m_info!=Success)
+        return;
+    
+      if(m_P.size()>0)
+        dest = m_Pinv * b;
+      else
+        dest = b;
+      
+      if(m_LDLt)
+      {
+        if(m_matrix.nonZeros()>0) // otherwise L==I
+          m_matrix.template triangularView<UnitLower>().solveInPlace(dest);
+      
+        dest = m_diag.asDiagonal().inverse() * dest;
+      
+        if (m_matrix.nonZeros()>0) // otherwise L==I
+          m_matrix.adjoint().template triangularView<UnitUpper>().solveInPlace(dest);
+      }
+      else
+      {
+        if(m_matrix.nonZeros()>0) // otherwise L==I
+          m_matrix.template triangularView<Lower>().solveInPlace(dest);
+      
+        if (m_matrix.nonZeros()>0) // otherwise L==I
+          m_matrix.adjoint().template triangularView<Upper>().solveInPlace(dest);
+      }
+      
+      if(m_P.size()>0)
+        dest = m_P * dest;
+    }
+    
+    /** \internal */
+    /*
+    template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
+    void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
+    {
+      // TODO
+    }
+    */
+    #endif // EIGEN_PARSED_BY_DOXYGEN
+
+  protected:
+    mutable ComputationInfo m_info;
+    bool m_isInitialized;
+    bool m_factorizationIsOk;
+    bool m_analysisIsOk;
+    bool m_LDLt;
+    
+    CholMatrixType m_matrix;
+    VectorType m_diag;                  // the diagonal coefficients in case of a LDLt decomposition
+    VectorXi m_parent;                  // elimination tree
+    VectorXi m_nonZerosPerCol;
+    PermutationMatrix<Dynamic> m_P;     // the permutation
+    PermutationMatrix<Dynamic> m_Pinv;  // the inverse permutation
+};
+
+namespace internal {
+  
+template<typename _MatrixType, int _UpLo, typename Rhs>
+struct solve_retval<SimplicialCholesky<_MatrixType,_UpLo>, Rhs>
+  : solve_retval_base<SimplicialCholesky<_MatrixType,_UpLo>, Rhs>
+{
+  typedef SimplicialCholesky<_MatrixType,_UpLo> Dec;
+  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    dec()._solve(rhs(),dst);
+  }
+};
+
+template<typename _MatrixType, int _UpLo, typename Rhs>
+struct sparse_solve_retval<SimplicialCholesky<_MatrixType,_UpLo>, Rhs>
+  : sparse_solve_retval_base<SimplicialCholesky<_MatrixType,_UpLo>, Rhs>
+{
+  typedef SimplicialCholesky<_MatrixType,_UpLo> Dec;
+  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    dec()._solve(rhs(),dst);
+  }
+};
+
+}
+
+#endif // EIGEN_SIMPLICIAL_CHOLESKY_H
diff --git a/unsupported/Eigen/src/SparseExtra/Solve.h b/unsupported/Eigen/src/SparseExtra/Solve.h
new file mode 100644
index 000000000..19449e9de
--- /dev/null
+++ b/unsupported/Eigen/src/SparseExtra/Solve.h
@@ -0,0 +1,82 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_SPARSE_SOLVE_H
+#define EIGEN_SPARSE_SOLVE_H
+
+namespace internal {
+
+template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval_base;
+template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval;
+  
+template<typename DecompositionType, typename Rhs>
+struct traits<sparse_solve_retval_base<DecompositionType, Rhs> >
+{
+  typedef typename DecompositionType::MatrixType MatrixType;
+  typedef SparseMatrix<typename Rhs::Scalar, Rhs::Options, typename Rhs::Index> ReturnType;
+};
+
+template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval_base
+ : public ReturnByValue<sparse_solve_retval_base<_DecompositionType, Rhs> >
+{
+  typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned;
+  typedef _DecompositionType DecompositionType;
+  typedef ReturnByValue<sparse_solve_retval_base> Base;
+  typedef typename Base::Index Index;
+
+  sparse_solve_retval_base(const DecompositionType& dec, const Rhs& rhs)
+    : m_dec(dec), m_rhs(rhs)
+  {}
+
+  inline Index rows() const { return m_dec.cols(); }
+  inline Index cols() const { return m_rhs.cols(); }
+  inline const DecompositionType& dec() const { return m_dec; }
+  inline const RhsNestedCleaned& rhs() const { return m_rhs; }
+
+  template<typename Dest> inline void evalTo(Dest& dst) const
+  {
+    static_cast<const sparse_solve_retval<DecompositionType,Rhs>*>(this)->evalTo(dst);
+  }
+
+  protected:
+    const DecompositionType& m_dec;
+    const typename Rhs::Nested m_rhs;
+};
+
+#define EIGEN_MAKE_SPARSE_SOLVE_HELPERS(DecompositionType,Rhs) \
+  typedef typename DecompositionType::MatrixType MatrixType; \
+  typedef typename MatrixType::Scalar Scalar; \
+  typedef typename MatrixType::RealScalar RealScalar; \
+  typedef typename MatrixType::Index Index; \
+  typedef Eigen::internal::sparse_solve_retval_base<DecompositionType,Rhs> Base; \
+  using Base::dec; \
+  using Base::rhs; \
+  using Base::rows; \
+  using Base::cols; \
+  sparse_solve_retval(const DecompositionType& dec, const Rhs& rhs) \
+    : Base(dec, rhs) {}
+    
+} // namepsace internal
+
+#endif // EIGEN_SPARSE_SOLVE_H
diff --git a/unsupported/Eigen/src/SparseExtra/SparseLDLT.h b/unsupported/Eigen/src/SparseExtra/SparseLDLTLegacy.h
index 837d70295..14a89dcf9 100644
--- a/unsupported/Eigen/src/SparseExtra/SparseLDLT.h
+++ b/unsupported/Eigen/src/SparseExtra/SparseLDLTLegacy.h
@@ -60,8 +60,8 @@ LDL License:
     and a notice that the code was modified is included.
  */
 
-#ifndef EIGEN_SPARSELDLT_H
-#define EIGEN_SPARSELDLT_H
+#ifndef EIGEN_SPARSELDLT_LEGACY_H
+#define EIGEN_SPARSELDLT_LEGACY_H
 
 /** \ingroup Sparse_Module
   *
@@ -187,6 +187,8 @@ class SparseLDLT
     VectorXi m_parent; // elimination tree
     VectorXi m_nonZerosPerCol;
 //     VectorXi m_w; // workspace
+    PermutationMatrix<Dynamic> m_P;
+    PermutationMatrix<Dynamic> m_Pinv;
     RealScalar m_precision;
     int m_flags;
     mutable int m_status;
@@ -248,15 +250,22 @@ void SparseLDLT<_MatrixType,Backend>::_symbolic(const _MatrixType& a)
   const Index* Ap = a._outerIndexPtr();
   const Index* Ai = a._innerIndexPtr();
   Index* Lp = m_matrix._outerIndexPtr();
+  
   const Index* P = 0;
   Index* Pinv = 0;
 
-  if (P)
+  if(P)
   {
-    /* If P is present then compute Pinv, the inverse of P */
-    for (Index k = 0; k < size; ++k)
-      Pinv[P[k]] = k;
+    m_P.indices()     = VectorXi::Map(P,size);
+    m_Pinv = m_P.inverse();
+    Pinv = m_Pinv.indices().data();
   }
+  else
+  {
+    m_P.resize(0);
+    m_Pinv.resize(0);
+  }
+
   for (Index k = 0; k < size; ++k)
   {
     /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
@@ -311,9 +320,16 @@ bool SparseLDLT<_MatrixType,Backend>::_numeric(const _MatrixType& a)
   Scalar * y = ei_aligned_stack_new(Scalar, size);
   Index * pattern = ei_aligned_stack_new(Index, size);
   Index * tags = ei_aligned_stack_new(Index, size);
-
-  const Index* P = 0;
-  const Index* Pinv = 0;
+  
+  Index* P = 0;
+  Index* Pinv = 0;
+  
+  if(m_P.size()==size)
+  {
+    P = m_P.indices().data();
+    Pinv = m_Pinv.indices().data();
+  }
+  
   bool ok = true;
 
   for (Index k = 0; k < size; ++k)
@@ -383,6 +399,9 @@ bool SparseLDLT<_MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) cons
   eigen_assert(m_matrix.rows()==b.rows());
   if (!m_succeeded)
     return false;
+  
+  if(m_P.size()>0)
+    b = m_Pinv * b;
 
   if (m_matrix.nonZeros()>0) // otherwise L==I
     m_matrix.template triangularView<UnitLower>().solveInPlace(b);
@@ -390,7 +409,10 @@ bool SparseLDLT<_MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) cons
   if (m_matrix.nonZeros()>0) // otherwise L==I
     m_matrix.adjoint().template triangularView<UnitUpper>().solveInPlace(b);
 
+  if(m_P.size()>0)
+    b = m_P * b;
+  
   return true;
 }
 
-#endif // EIGEN_SPARSELDLT_H
+#endif // EIGEN_SPARSELDLT_LEGACY_H
diff --git a/unsupported/test/sparse_ldlt.cpp b/unsupported/test/sparse_ldlt.cpp
index 275839670..035e21123 100644
--- a/unsupported/test/sparse_ldlt.cpp
+++ b/unsupported/test/sparse_ldlt.cpp
@@ -31,6 +31,8 @@
 
 template<typename Scalar> void sparse_ldlt(int rows, int cols)
 {
+  static bool odd = true;
+  odd = !odd;
   double density = std::max(8./(rows*cols), 0.01);
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
   typedef Matrix<Scalar,Dynamic,1> DenseVector;
@@ -42,41 +44,126 @@ template<typename Scalar> void sparse_ldlt(int rows, int cols)
   DenseVector refX(cols), x(cols);
 
   initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0);
-  for(int i=0; i<rows; ++i)
-    m2.coeffRef(i,i) = refMat2(i,i) = internal::abs(internal::real(refMat2(i,i)));
-
-  refX = refMat2.template selfadjointView<Upper>().ldlt().solve(b);
+  
+  SparseMatrix<Scalar> m3 = m2 * m2.adjoint(), m3_lo(rows,rows), m3_up(rows,rows);
+  DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
+  
+  refX = refMat3.template selfadjointView<Upper>().ldlt().solve(b);
   typedef SparseMatrix<Scalar,Upper|SelfAdjoint> SparseSelfAdjointMatrix;
   x = b;
-  SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
+  SparseLDLT<SparseSelfAdjointMatrix> ldlt(m3);
   if (ldlt.succeeded())
     ldlt.solveInPlace(x);
   else
     std::cerr << "warning LDLT failed\n";
 
-  VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
+//   VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
   VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default");
+  
+  
+  
+  // new Simplicial LLT
+  
+  
+  // new API
+  {
+    SparseMatrix<Scalar> m2(rows, cols);
+    DenseMatrix refMat2(rows, cols);
 
-#ifdef EIGEN_CHOLMOD_SUPPORT
-  x = b;
-  SparseLDLT<SparseSelfAdjointMatrix, Cholmod> ldlt2(m2);
-  if (ldlt2.succeeded())
-    ldlt2.solveInPlace(x);
-  else
-    std::cerr << "warning LDLT failed\n";
+    DenseVector b = DenseVector::Random(cols);
+    DenseVector ref_x(cols), x(cols);
+    DenseMatrix B = DenseMatrix::Random(rows,cols);
+    DenseMatrix ref_X(rows,cols), X(rows,cols);
 
-  VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
-  VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solveInPlace");
+    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, 0, 0);
 
+    for(int i=0; i<rows; ++i)
+      m2.coeffRef(i,i) = refMat2(i,i) = internal::abs(internal::real(refMat2(i,i)));
+    
+    
+    SparseMatrix<Scalar> m3 = m2 * m2.adjoint(), m3_lo(rows,rows), m3_up(rows,rows);
+    DenseMatrix refMat3 = refMat2 * refMat2.adjoint();
+    
+    m3_lo.template selfadjointView<Lower>().rankUpdate(m2,0);
+    m3_up.template selfadjointView<Upper>().rankUpdate(m2,0);
+    
+    // with a single vector as the rhs
+    ref_x = refMat3.template selfadjointView<Lower>().llt().solve(b);
 
-  SparseLDLT<SparseSelfAdjointMatrix, Cholmod> ldlt3(m2);
-  if (ldlt3.succeeded())
-    x = ldlt3.solve(b);
-  else
-    std::cerr << "warning LDLT failed\n";
+    x = SimplicialCholesky<SparseMatrix<Scalar>, Lower>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(b);
+    VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, lower, single dense rhs");
+    
+    x = SimplicialCholesky<SparseMatrix<Scalar>, Upper>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(b);
+    VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, upper, single dense rhs");
+    
+//     x = SimplicialCholesky<SparseMatrix<Scalar>, Lower>(m3_lo).solve(b);
+//     VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, lower only, single dense rhs");
+    
+//     x = SimplicialCholesky<SparseMatrix<Scalar>, Upper>(m3_up).solve(b);
+//     VERIFY(ref_x.isApprox(x,test_precision<Scalar>()) && "SimplicialCholesky: solve, upper only, single dense rhs");
+    
+    
+    // with multiple rhs
+    ref_X = refMat3.template selfadjointView<Lower>().llt().solve(B);
+
+    X = SimplicialCholesky<SparseMatrix<Scalar>, Lower>()/*.setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt)*/.compute(m3).solve(B);
+    VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, lower, multiple dense rhs");
+    
+//     X = SimplicialCholesky<SparseMatrix<Scalar>, Upper>().setMode(odd ? SimplicialCholeskyLLt : SimplicialCholeskyLDLt).compute(m3).solve(B);
+//     VERIFY(ref_X.isApprox(X,test_precision<Scalar>()) && "SimplicialCholesky: solve, full storage, upper, multiple dense rhs");
+    
+    
+//     // with a sparse rhs
+//     SparseMatrix<Scalar> spB(rows,cols), spX(rows,cols);
+//     B.diagonal().array() += 1;
+//     spB = B.sparseView(0.5,1);
+//     
+//     ref_X = refMat3.template selfadjointView<Lower>().llt().solve(DenseMatrix(spB));
 
-  VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
-  VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solve");
+//     spX = SimplicialCholesky<SparseMatrix<Scalar>, Lower>(m3).solve(spB);
+//     VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
+//     
+//     spX = SimplicialCholesky<SparseMatrix<Scalar>, Upper>(m3).solve(spB);
+//     VERIFY(ref_X.isApprox(spX.toDense(),test_precision<Scalar>()) && "LLT: cholmod solve, multiple sparse rhs");
+  }
+  
+  
+  
+//   for(int i=0; i<rows; ++i)
+//     m2.coeffRef(i,i) = refMat2(i,i) = internal::abs(internal::real(refMat2(i,i)));
+// 
+//   refX = refMat2.template selfadjointView<Upper>().ldlt().solve(b);
+//   typedef SparseMatrix<Scalar,Upper|SelfAdjoint> SparseSelfAdjointMatrix;
+//   x = b;
+//   SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
+//   if (ldlt.succeeded())
+//     ldlt.solveInPlace(x);
+//   else
+//     std::cerr << "warning LDLT failed\n";
+// 
+//   VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
+//   VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default");
+
+#ifdef EIGEN_CHOLMOD_SUPPORT
+//   x = b;
+//   SparseLDLT<SparseSelfAdjointMatrix, Cholmod> ldlt2(m2);
+//   if (ldlt2.succeeded())
+//     ldlt2.solveInPlace(x);
+//   else
+//     std::cerr << "warning LDLT failed\n";
+// 
+//   VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
+//   VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solveInPlace");
+// 
+// 
+//   SparseLDLT<SparseSelfAdjointMatrix, Cholmod> ldlt3(m2);
+//   if (ldlt3.succeeded())
+//     x = ldlt3.solve(b);
+//   else
+//     std::cerr << "warning LDLT failed\n";
+// 
+//   VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
+//   VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: cholmod solve");
 
 #endif
 }