move sparse solvers from unsupported/ to main Eigen/ and remove the "not stable yet" warning

9 files changed, 0 insertions, 4350 deletions

diff --git a/unsupported/Eigen/src/SparseExtra/Amd.h b/unsupported/Eigen/src/SparseExtra/Amd.h
deleted file mode 100644
index 3cf8bd1e1..000000000
--- a/unsupported/Eigen/src/SparseExtra/Amd.h
+++ /dev/null
@@ -1,448 +0,0 @@
-// 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;
-}
-
-
-/** \internal
-  * Approximate minimum degree ordering algorithm.
-  * \returns the permutation P reducing the fill-in of the input matrix \a C
-  * The input matrix \a C must be a selfadjoint compressed column major SparseMatrix object. Both the upper and lower parts have to be stored, but the diagonal entries are optional.
-  * On exit the values of C are destroyed */
-template<typename Scalar, typename Index>
-void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic,Dynamic,Index>& perm)
-{
-  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;
-  
-  Index n = C.cols();
-  dense = std::max<Index> (16, 10 * sqrt ((double) n));   /* find dense threshold */
-  dense = std::min<Index> (n-2, dense);
-  
-  Index cnz = C.nonZeros();
-  perm.resize(n+1);
-  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    = perm.indices().data();                              /* 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<Index>(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<Index>(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<Index>(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<Index>(i, k, head, next, perm.indices().data(), w);
-  }
-  
-  perm.indices().conservativeResize(n);
-
-  delete[] W;
-}
-
-} // namespace internal
-
-#endif // EIGEN_SPARSE_AMD_H
diff --git a/unsupported/Eigen/src/SparseExtra/CholmodSupport.h b/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
deleted file mode 100644
index 3e502c0aa..000000000
--- a/unsupported/Eigen/src/SparseExtra/CholmodSupport.h
+++ /dev/null
@@ -1,399 +0,0 @@
-// 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/>.
-
-#ifndef EIGEN_CHOLMODSUPPORT_H
-#define EIGEN_CHOLMODSUPPORT_H
-
-namespace internal {
-
-template<typename Scalar, typename CholmodType>
-void cholmod_configure_matrix(CholmodType& mat)
-{
-  if (internal::is_same<Scalar,float>::value)
-  {
-    mat.xtype = CHOLMOD_REAL;
-    mat.dtype = CHOLMOD_SINGLE;
-  }
-  else if (internal::is_same<Scalar,double>::value)
-  {
-    mat.xtype = CHOLMOD_REAL;
-    mat.dtype = CHOLMOD_DOUBLE;
-  }
-  else if (internal::is_same<Scalar,std::complex<float> >::value)
-  {
-    mat.xtype = CHOLMOD_COMPLEX;
-    mat.dtype = CHOLMOD_SINGLE;
-  }
-  else if (internal::is_same<Scalar,std::complex<double> >::value)
-  {
-    mat.xtype = CHOLMOD_COMPLEX;
-    mat.dtype = CHOLMOD_DOUBLE;
-  }
-  else
-  {
-    eigen_assert(false && "Scalar type not supported by CHOLMOD");
-  }
-}
-
-} // namespace internal
-
-/** Wraps the Eigen sparse matrix \a mat into a Cholmod sparse matrix object.
-  * Note that the data are shared.
-  */
-template<typename _Scalar, int _Options, typename _Index>
-cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
-{
-  typedef SparseMatrix<_Scalar,_Options,_Index> MatrixType;
-  cholmod_sparse res;
-  res.nzmax   = mat.nonZeros();
-  res.nrow    = mat.rows();;
-  res.ncol    = mat.cols();
-  res.p       = mat._outerIndexPtr();
-  res.i       = mat._innerIndexPtr();
-  res.x       = mat._valuePtr();
-  res.sorted  = 1;
-  res.packed  = 1;
-  res.dtype   = 0;
-  res.stype   = -1;
-  
-  if (internal::is_same<_Index,int>::value)
-  {
-    res.itype = CHOLMOD_INT;
-  }
-  else
-  {
-    eigen_assert(false && "Index type different than int is not supported yet");
-  }
-
-  // setup res.xtype
-  internal::cholmod_configure_matrix<_Scalar>(res);
-  
-  res.stype = 0;
-  
-  return res;
-}
-
-template<typename _Scalar, int _Options, typename _Index>
-const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
-{
-  cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
-  return res;
-}
-
-/** Returns a view of the Eigen sparse matrix \a mat as Cholmod sparse matrix.
-  * The data are not copied but shared. */
-template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
-cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
-{
-  cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
-  
-  if(UpLo==Upper) res.stype =  1;
-  if(UpLo==Lower) res.stype = -1;
-
-  return res;
-}
-
-/** Returns a view of the Eigen \b dense matrix \a mat as Cholmod dense matrix.
-  * The data are not copied but shared. */
-template<typename Derived>
-cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
-{
-  EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
-  typedef typename Derived::Scalar Scalar;
-
-  cholmod_dense res;
-  res.nrow   = mat.rows();
-  res.ncol   = mat.cols();
-  res.nzmax  = res.nrow * res.ncol;
-  res.d      = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
-  res.x      = mat.derived().data();
-  res.z      = 0;
-
-  internal::cholmod_configure_matrix<Scalar>(res);
-
-  return res;
-}
-
-/** Returns a view of the Cholmod sparse matrix \a cm as an Eigen sparse matrix.
-  * The data are not copied but shared. */
-template<typename Scalar, int Flags, typename Index>
-MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm)
-{
-  return MappedSparseMatrix<Scalar,Flags,Index>
-         (cm.nrow, cm.ncol, reinterpret_cast<Index*>(cm.p)[cm.ncol],
-          reinterpret_cast<Index*>(cm.p), reinterpret_cast<Index*>(cm.i),reinterpret_cast<Scalar*>(cm.x) );
-}
-
-enum CholmodMode {
-  CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
-};
-
-/** \brief A Cholesky factorization and solver based on Cholmod
-  *
-  * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization
-  * using the Cholmod library. 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 CholmodDecomposition
-{
-  public:
-    typedef _MatrixType MatrixType;
-    enum { UpLo = _UpLo };
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
-    typedef MatrixType CholMatrixType;
-    typedef typename MatrixType::Index Index;
-
-  public:
-
-    CholmodDecomposition()
-      : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
-    {
-      cholmod_start(&m_cholmod);
-      setMode(CholmodLDLt);
-    }
-
-    CholmodDecomposition(const MatrixType& matrix)
-      : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
-    {
-      cholmod_start(&m_cholmod);
-      compute(matrix);
-    }
-
-    ~CholmodDecomposition()
-    {
-      if(m_cholmodFactor)
-        cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
-      cholmod_finish(&m_cholmod);
-    }
-    
-    inline Index cols() const { return m_cholmodFactor->n; }
-    inline Index rows() const { return m_cholmodFactor->n; }
-    
-    void setMode(CholmodMode mode)
-    {
-      switch(mode)
-      {
-        case CholmodAuto:
-          m_cholmod.final_asis = 1;
-          m_cholmod.supernodal = CHOLMOD_AUTO;
-          break;
-        case CholmodSimplicialLLt:
-          m_cholmod.final_asis = 0;
-          m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
-          m_cholmod.final_ll = 1;
-          break;
-        case CholmodSupernodalLLt:
-          m_cholmod.final_asis = 1;
-          m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
-          break;
-        case CholmodLDLt:
-          m_cholmod.final_asis = 1;
-          m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
-          break;
-        default:
-          break;
-      }
-    }
-    
-    /** \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 */
-    void compute(const MatrixType& matrix)
-    {
-      analyzePattern(matrix);
-      factorize(matrix);
-    }
-    
-    /** \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<CholmodDecomposition, Rhs>
-    solve(const MatrixBase<Rhs>& b) const
-    {
-      eigen_assert(m_isInitialized && "LLT is not initialized.");
-      eigen_assert(rows()==b.rows()
-                && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
-      return internal::solve_retval<CholmodDecomposition, 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 && "LLT is not initialized.");
-      eigen_assert(rows()==b.rows()
-                && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
-      return internal::sparse_solve_retval<CholmodDecomposition, 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& matrix)
-    {
-      if(m_cholmodFactor)
-      {
-        cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
-        m_cholmodFactor = 0;
-      }
-      cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
-      m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
-      
-      this->m_isInitialized = true;
-      this->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& matrix)
-    {
-      eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
-      cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
-      cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
-      
-      this->m_info = Success;
-      m_factorizationIsOk = true;
-    }
-    
-    /** Returns a reference to the Cholmod's configuration structure to get a full control over the performed operations.
-     *  See the Cholmod user guide for details. */
-    cholmod_common& cholmod() { return m_cholmod; }
-    
-    #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()");
-      const Index size = m_cholmodFactor->n;
-      eigen_assert(size==b.rows());
-
-      // note: cd stands for Cholmod Dense
-      cholmod_dense b_cd = viewAsCholmod(b.const_cast_derived());
-      cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
-      if(!x_cd)
-      {
-        this->m_info = NumericalIssue;
-      }
-      // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
-      dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
-      cholmod_free_dense(&x_cd, &m_cholmod);
-    }
-    
-    /** \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
-    {
-      eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
-      const Index size = m_cholmodFactor->n;
-      eigen_assert(size==b.rows());
-
-      // note: cs stands for Cholmod Sparse
-      cholmod_sparse b_cs = viewAsCholmod(b);
-      cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
-      if(!x_cs)
-      {
-        this->m_info = NumericalIssue;
-      }
-      // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
-      dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
-      cholmod_free_sparse(&x_cs, &m_cholmod);
-    }
-    #endif // EIGEN_PARSED_BY_DOXYGEN
-    
-    template<typename Stream>
-    void dumpMemory(Stream& s)
-    {}
-
-  protected:
-    mutable cholmod_common m_cholmod;
-    cholmod_factor* m_cholmodFactor;
-    mutable ComputationInfo m_info;
-    bool m_isInitialized;
-    int m_factorizationIsOk;
-    int m_analysisIsOk;
-};
-
-namespace internal {
-  
-template<typename _MatrixType, int _UpLo, typename Rhs>
-struct solve_retval<CholmodDecomposition<_MatrixType,_UpLo>, Rhs>
-  : solve_retval_base<CholmodDecomposition<_MatrixType,_UpLo>, Rhs>
-{
-  typedef CholmodDecomposition<_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<CholmodDecomposition<_MatrixType,_UpLo>, Rhs>
-  : sparse_solve_retval_base<CholmodDecomposition<_MatrixType,_UpLo>, Rhs>
-{
-  typedef CholmodDecomposition<_MatrixType,_UpLo> Dec;
-  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
-
-  template<typename Dest> void evalTo(Dest& dst) const
-  {
-    dec()._solve(rhs(),dst);
-  }
-};
-
-}
-
-#endif // EIGEN_CHOLMODSUPPORT_H
diff --git a/unsupported/Eigen/src/SparseExtra/CholmodSupportLegacy.h b/unsupported/Eigen/src/SparseExtra/CholmodSupportLegacy.h
deleted file mode 100644
index 33af6a176..000000000
--- a/unsupported/Eigen/src/SparseExtra/CholmodSupportLegacy.h
+++ /dev/null
@@ -1,520 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008-2009 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_CHOLMODSUPPORT_LEGACY_H
-#define EIGEN_CHOLMODSUPPORT_LEGACY_H
-
-namespace internal {
-
-template<typename Scalar, typename CholmodType>
-void cholmod_configure_matrix_legacy(CholmodType& mat)
-{
-  if (internal::is_same<Scalar,float>::value)
-  {
-    mat.xtype = CHOLMOD_REAL;
-    mat.dtype = CHOLMOD_SINGLE;
-  }
-  else if (internal::is_same<Scalar,double>::value)
-  {
-    mat.xtype = CHOLMOD_REAL;
-    mat.dtype = CHOLMOD_DOUBLE;
-  }
-  else if (internal::is_same<Scalar,std::complex<float> >::value)
-  {
-    mat.xtype = CHOLMOD_COMPLEX;
-    mat.dtype = CHOLMOD_SINGLE;
-  }
-  else if (internal::is_same<Scalar,std::complex<double> >::value)
-  {
-    mat.xtype = CHOLMOD_COMPLEX;
-    mat.dtype = CHOLMOD_DOUBLE;
-  }
-  else
-  {
-    eigen_assert(false && "Scalar type not supported by CHOLMOD");
-  }
-}
-
-template<typename _MatrixType>
-cholmod_sparse cholmod_map_eigen_to_sparse(_MatrixType& mat)
-{
-  typedef typename _MatrixType::Scalar Scalar;
-  cholmod_sparse res;
-  res.nzmax   = mat.nonZeros();
-  res.nrow    = mat.rows();;
-  res.ncol    = mat.cols();
-  res.p       = mat._outerIndexPtr();
-  res.i       = mat._innerIndexPtr();
-  res.x       = mat._valuePtr();
-  res.xtype   = CHOLMOD_REAL;
-  res.itype   = CHOLMOD_INT;
-  res.sorted  = 1;
-  res.packed  = 1;
-  res.dtype   = 0;
-  res.stype   = -1;
-
-  internal::cholmod_configure_matrix_legacy<Scalar>(res);
-
-
-  if (_MatrixType::Flags & SelfAdjoint)
-  {
-    if (_MatrixType::Flags & Upper)
-      res.stype = 1;
-    else if (_MatrixType::Flags & Lower)
-      res.stype = -1;
-    else
-      res.stype = 0;
-  }
-  else
-    res.stype = -1; // by default we consider the lower part
-
-  return res;
-}
-
-template<typename Derived>
-cholmod_dense cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat)
-{
-  EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
-  typedef typename Derived::Scalar Scalar;
-
-  cholmod_dense res;
-  res.nrow   = mat.rows();
-  res.ncol   = mat.cols();
-  res.nzmax  = res.nrow * res.ncol;
-  res.d      = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
-  res.x      = mat.derived().data();
-  res.z      = 0;
-
-  internal::cholmod_configure_matrix_legacy<Scalar>(res);
-
-  return res;
-}
-
-template<typename Scalar, int Flags, typename Index>
-MappedSparseMatrix<Scalar,Flags,Index> map_cholmod_sparse_to_eigen(cholmod_sparse& cm)
-{
-  return MappedSparseMatrix<Scalar,Flags,Index>
-         (cm.nrow, cm.ncol, reinterpret_cast<Index*>(cm.p)[cm.ncol],
-          reinterpret_cast<Index*>(cm.p), reinterpret_cast<Index*>(cm.i),reinterpret_cast<Scalar*>(cm.x) );
-}
-
-} // namespace internal
-
-/** \deprecated use class SimplicialLDLT, or class SimplicialLLT, class ConjugateGradient */
-template<typename _MatrixType>
-class SparseLLT<_MatrixType, Cholmod> : public SparseLLT<_MatrixType>
-{
-  protected:
-    typedef SparseLLT<_MatrixType> Base;
-    typedef typename Base::Scalar Scalar;
-    typedef typename Base::RealScalar RealScalar;
-    typedef typename Base::CholMatrixType CholMatrixType;
-    using Base::MatrixLIsDirty;
-    using Base::SupernodalFactorIsDirty;
-    using Base::m_flags;
-    using Base::m_matrix;
-    using Base::m_status;
-
-  public:
-    typedef _MatrixType MatrixType;
-    typedef typename MatrixType::Index Index;
-
-    /** \deprecated the entire class is deprecated */
-    EIGEN_DEPRECATED SparseLLT(int flags = 0)
-      : Base(flags), m_cholmodFactor(0)
-    {
-      cholmod_start(&m_cholmod);
-    }
-
-    /** \deprecated the entire class is deprecated */
-    EIGEN_DEPRECATED SparseLLT(const MatrixType& matrix, int flags = 0)
-      : Base(flags), m_cholmodFactor(0)
-    {
-      cholmod_start(&m_cholmod);
-      compute(matrix);
-    }
-
-    ~SparseLLT()
-    {
-      if (m_cholmodFactor)
-        cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
-      cholmod_finish(&m_cholmod);
-    }
-
-    inline const CholMatrixType& matrixL() const;
-
-    template<typename Derived>
-    bool solveInPlace(MatrixBase<Derived> &b) const;
-
-    template<typename Rhs>
-    inline const internal::solve_retval<SparseLLT<MatrixType, Cholmod>, Rhs>
-    solve(const MatrixBase<Rhs>& b) const
-    {
-      eigen_assert(true && "SparseLLT is not initialized.");
-      return internal::solve_retval<SparseLLT<MatrixType, Cholmod>, Rhs>(*this, b.derived());
-    }
-
-    void compute(const MatrixType& matrix);
-
-    inline Index cols() const { return m_matrix.cols(); }
-    inline Index rows() const { return m_matrix.rows(); }
-
-    inline const cholmod_factor* cholmodFactor() const
-    { return m_cholmodFactor; }
-
-    inline cholmod_common* cholmodCommon() const
-    { return &m_cholmod; }
-
-    bool succeeded() const;
-
-  protected:
-    mutable cholmod_common m_cholmod;
-    cholmod_factor* m_cholmodFactor;
-};
-
-
-namespace internal {
-
-template<typename _MatrixType, typename Rhs>
-  struct solve_retval<SparseLLT<_MatrixType, Cholmod>, Rhs>
-  : solve_retval_base<SparseLLT<_MatrixType, Cholmod>, Rhs>
-{
-  typedef SparseLLT<_MatrixType, Cholmod> SpLLTDecType;
-  EIGEN_MAKE_SOLVE_HELPERS(SpLLTDecType,Rhs)
-
-  template<typename Dest> void evalTo(Dest& dst) const
-  {
-    //Index size = dec().cholmodFactor()->n;
-    eigen_assert((Index)dec().cholmodFactor()->n==rhs().rows());
-    
-    cholmod_factor* cholmodFactor = const_cast<cholmod_factor*>(dec().cholmodFactor());
-    cholmod_common* cholmodCommon = const_cast<cholmod_common*>(dec().cholmodCommon());
-    // this uses Eigen's triangular sparse solver
-    // if (m_status & MatrixLIsDirty)
-    //   matrixL();
-    // Base::solveInPlace(b);
-    // as long as our own triangular sparse solver is not fully optimal,
-    // let's use CHOLMOD's one:
-    cholmod_dense cdb = internal::cholmod_map_eigen_to_dense(rhs().const_cast_derived());
-    cholmod_dense* x = cholmod_solve(CHOLMOD_A, cholmodFactor, &cdb, cholmodCommon);
-
-    dst = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x), rhs().rows());  
-
-    cholmod_free_dense(&x, cholmodCommon);
-
-  }
-    
-};
-
-} // namespace internal
-
-
-
-template<typename _MatrixType>
-void SparseLLT<_MatrixType,Cholmod>::compute(const _MatrixType& a)
-{
-  if (m_cholmodFactor)
-  {
-    cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
-    m_cholmodFactor = 0;
-  }
-
-  cholmod_sparse A = internal::cholmod_map_eigen_to_sparse(const_cast<_MatrixType&>(a));
-//   m_cholmod.supernodal = CHOLMOD_AUTO;
-  // TODO
-//   if (m_flags&IncompleteFactorization)
-//   {
-//     m_cholmod.nmethods = 1;
-//     m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
-//     m_cholmod.postorder = 0;
-//   }
-//   else
-//   {
-//     m_cholmod.nmethods = 1;
-//     m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
-//     m_cholmod.postorder = 0;
-//   }
-//   m_cholmod.final_ll = 1;
-  m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
-  cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
-
-  this->m_status = (this->m_status & ~Base::SupernodalFactorIsDirty) | Base::MatrixLIsDirty;
-}
-
-
-// TODO
-template<typename _MatrixType>
-bool SparseLLT<_MatrixType,Cholmod>::succeeded() const
-{ return true; }
-
-
-
-template<typename _MatrixType>
-inline const typename SparseLLT<_MatrixType,Cholmod>::CholMatrixType&
-SparseLLT<_MatrixType,Cholmod>::matrixL() const
-{
-  if (this->m_status & Base::MatrixLIsDirty)
-  {
-    eigen_assert(!(this->m_status & Base::SupernodalFactorIsDirty));
-
-    cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod);
-    const_cast<typename Base::CholMatrixType&>(this->m_matrix) = 
-      internal::map_cholmod_sparse_to_eigen<Scalar,ColMajor,Index>(*cmRes);
-    free(cmRes);
-
-    this->m_status = (this->m_status & ~Base::MatrixLIsDirty);
-  }
-  return this->m_matrix;
-}
-
-
-
-
-template<typename _MatrixType>
-template<typename Derived>
-bool SparseLLT<_MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
-{
-  //Index size = m_cholmodFactor->n;
-  eigen_assert((Index)m_cholmodFactor->n==b.rows());
-
-  // this uses Eigen's triangular sparse solver
-  //   if (m_status & MatrixLIsDirty)
-  //     matrixL();
-  //   Base::solveInPlace(b);
-  // as long as our own triangular sparse solver is not fully optimal,
-  // let's use CHOLMOD's one:
-  cholmod_dense cdb = internal::cholmod_map_eigen_to_dense(b);
-
-  cholmod_dense* x = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &cdb, &m_cholmod);
-  eigen_assert(x && "Eigen: cholmod_solve failed.");
-
-  b = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x),b.rows());
-  cholmod_free_dense(&x, &m_cholmod);
-  return true;
-}
-
-
-
-
-
-
-
-
-
-
-
-template<typename _MatrixType>
-class SparseLDLT<_MatrixType,Cholmod> : public SparseLDLT<_MatrixType>
-{
-  protected:
-    typedef SparseLDLT<_MatrixType> Base;
-    typedef typename Base::Scalar Scalar;
-    typedef typename Base::RealScalar RealScalar;
-    using Base::MatrixLIsDirty;
-    using Base::SupernodalFactorIsDirty;
-    using Base::m_flags;
-    using Base::m_matrix;
-    using Base::m_status;
-
-  public:
-    typedef _MatrixType MatrixType;
-    typedef typename MatrixType::Index Index;
-
-    SparseLDLT(int flags = 0)
-      : Base(flags), m_cholmodFactor(0)
-    {
-      cholmod_start(&m_cholmod);
-    }
-
-    SparseLDLT(const _MatrixType& matrix, int flags = 0)
-      : Base(flags), m_cholmodFactor(0)
-    {
-      cholmod_start(&m_cholmod);
-      compute(matrix);
-    }
-
-    ~SparseLDLT()
-    {
-      if (m_cholmodFactor)
-        cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
-      cholmod_finish(&m_cholmod);
-    }
-
-    inline const typename Base::CholMatrixType& matrixL(void) const;
-
-    template<typename Derived>
-    void solveInPlace(MatrixBase<Derived> &b) const;
-
-    template<typename Rhs>
-    inline const internal::solve_retval<SparseLDLT<MatrixType, Cholmod>, Rhs>
-    solve(const MatrixBase<Rhs>& b) const
-    {
-      eigen_assert(true && "SparseLDLT is not initialized.");
-      return internal::solve_retval<SparseLDLT<MatrixType, Cholmod>, Rhs>(*this, b.derived());
-    }
-
-    void compute(const _MatrixType& matrix);
-
-    inline Index cols() const { return m_matrix.cols(); }
-    inline Index rows() const { return m_matrix.rows(); }
-
-    inline const cholmod_factor* cholmodFactor() const
-    { return m_cholmodFactor; }
-
-    inline cholmod_common* cholmodCommon() const
-    { return &m_cholmod; }
-
-    bool succeeded() const;
-
-  protected:
-    mutable cholmod_common m_cholmod;
-    cholmod_factor* m_cholmodFactor;
-};
-
-
-
-namespace internal {
-
-template<typename _MatrixType, typename Rhs>
-  struct solve_retval<SparseLDLT<_MatrixType, Cholmod>, Rhs>
-  : solve_retval_base<SparseLDLT<_MatrixType, Cholmod>, Rhs>
-{
-  typedef SparseLDLT<_MatrixType, Cholmod> SpLDLTDecType;
-  EIGEN_MAKE_SOLVE_HELPERS(SpLDLTDecType,Rhs)
-
-  template<typename Dest> void evalTo(Dest& dst) const
-  {
-    //Index size = dec().cholmodFactor()->n;
-    eigen_assert((Index)dec().cholmodFactor()->n==rhs().rows());
-    
-    cholmod_factor* cholmodFactor = const_cast<cholmod_factor*>(dec().cholmodFactor());
-    cholmod_common* cholmodCommon = const_cast<cholmod_common*>(dec().cholmodCommon());
-    // this uses Eigen's triangular sparse solver
-    // if (m_status & MatrixLIsDirty)
-    //   matrixL();
-    // Base::solveInPlace(b);
-    // as long as our own triangular sparse solver is not fully optimal,
-    // let's use CHOLMOD's one:
-    cholmod_dense cdb = internal::cholmod_map_eigen_to_dense(rhs().const_cast_derived());
-    cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, cholmodFactor, &cdb, cholmodCommon);
-
-    dst = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x), rhs().rows());  
-    cholmod_free_dense(&x, cholmodCommon);
-
-  }
-    
-};
-
-
-} // namespace internal
-
-template<typename _MatrixType>
-void SparseLDLT<_MatrixType,Cholmod>::compute(const _MatrixType& a)
-{
-  if (m_cholmodFactor)
-  {
-    cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
-    m_cholmodFactor = 0;
-  }
-
-  cholmod_sparse A = internal::cholmod_map_eigen_to_sparse(const_cast<_MatrixType&>(a));
- 
-  //m_cholmod.supernodal = CHOLMOD_AUTO;
-  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
-  //m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
-  // TODO
-  if (this->m_flags & IncompleteFactorization)
-  {
-    m_cholmod.nmethods = 1;
-    //m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
-    m_cholmod.method[0].ordering = CHOLMOD_COLAMD;
-    m_cholmod.postorder = 1;
-  }
-  else
-  {
-    m_cholmod.nmethods = 1;
-    m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
-    m_cholmod.postorder = 0;
-  }
-  m_cholmod.final_ll = 0;
-  m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
-  cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
-  
-  this->m_status = (this->m_status & ~Base::SupernodalFactorIsDirty) | Base::MatrixLIsDirty;
-}
-
-
-// TODO
-template<typename _MatrixType>
-bool SparseLDLT<_MatrixType,Cholmod>::succeeded() const
-{ return true; }
-
-
-template<typename _MatrixType>
-inline const typename SparseLDLT<_MatrixType>::CholMatrixType&
-SparseLDLT<_MatrixType,Cholmod>::matrixL() const
-{
-  if (this->m_status & Base::MatrixLIsDirty)
-  {
-    eigen_assert(!(this->m_status & Base::SupernodalFactorIsDirty));
-
-    cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod);
-    const_cast<typename Base::CholMatrixType&>(this->m_matrix) = MappedSparseMatrix<Scalar>(*cmRes);
-    free(cmRes);
-
-    this->m_status = (this->m_status & ~Base::MatrixLIsDirty);
-  }
-  return this->m_matrix;
-}
-
-
-
-
-
-
-template<typename _MatrixType>
-template<typename Derived>
-void SparseLDLT<_MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
-{
-  //Index size = m_cholmodFactor->n;
-  eigen_assert((Index)m_cholmodFactor->n == b.rows());
-
-  // this uses Eigen's triangular sparse solver
-  //   if (m_status & MatrixLIsDirty)
-  //     matrixL();
-  //   Base::solveInPlace(b);
-  // as long as our own triangular sparse solver is not fully optimal,
-  // let's use CHOLMOD's one:
-  cholmod_dense cdb = internal::cholmod_map_eigen_to_dense(b);
-  cholmod_dense* x = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &cdb, &m_cholmod);
-  b = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x),b.rows());
-  cholmod_free_dense(&x, &m_cholmod);
-}
-
-
-
-
-
-
-#endif // EIGEN_CHOLMODSUPPORT_LEGACY_H
diff --git a/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h b/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h
deleted file mode 100644
index 2147af258..000000000
--- a/unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h
+++ /dev/null
@@ -1,802 +0,0 @@
-// 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 factorizations
-  *
-  * These classes provide LL^T and LDL^T Cholesky factorizations of sparse matrices that are
-  * selfadjoint and positive definite. The factorization allows for solving A.X = B where
-  * 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 Derived>
-class SimplicialCholeskyBase
-{
-  public:
-    typedef typename internal::traits<Derived>::MatrixType MatrixType;
-    enum { UpLo = internal::traits<Derived>::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,Dynamic,1> VectorType;
-
-  public:
-
-    SimplicialCholeskyBase()
-      : m_info(Success), m_isInitialized(false)
-    {}
-
-    SimplicialCholeskyBase(const MatrixType& matrix)
-      : m_info(Success), m_isInitialized(false)
-    {
-      compute(matrix);
-    }
-
-    ~SimplicialCholeskyBase()
-    {
-    }
-
-    Derived& derived() { return *static_cast<Derived*>(this); }
-    const Derived& derived() const { return *static_cast<const Derived*>(this); }
-    
-    inline Index cols() const { return m_matrix.cols(); }
-    inline Index rows() const { return m_matrix.rows(); }
-    
-    /** \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 */
-    Derived& compute(const MatrixType& matrix)
-    {
-      derived().analyzePattern(matrix);
-      derived().factorize(matrix);
-      return 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::solve_retval<SimplicialCholeskyBase, Rhs>
-    solve(const MatrixBase<Rhs>& b) const
-    {
-      eigen_assert(m_isInitialized && "Simplicial LLt or LDLt is not initialized.");
-      eigen_assert(rows()==b.rows()
-                && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
-      return internal::solve_retval<SimplicialCholeskyBase, 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<SimplicialCholeskyBase, Rhs>
-    solve(const SparseMatrixBase<Rhs>& b) const
-    {
-      eigen_assert(m_isInitialized && "Simplicial LLt or LDLt 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<SimplicialCholeskyBase, Rhs>(*this, b.derived());
-    }
-    
-    /** \returns the permutation P
-      * \sa permutationPinv() */
-    const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const
-    { return m_P; }
-    
-    /** \returns the inverse P^-1 of the permutation P
-      * \sa permutationP() */
-    const PermutationMatrix<Dynamic,Dynamic,Index>& permutationPinv() const
-    { return m_Pinv; }
-
-#ifndef EIGEN_PARSED_BY_DOXYGEN
-    /** \internal */
-    template<typename Stream>
-    void dumpMemory(Stream& s)
-    {
-      int total = 0;
-      s << "  L:        " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
-      s << "  diag:     " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
-      s << "  tree:     " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
-      s << "  nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
-      s << "  perm:     " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
-      s << "  perm^-1:  " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
-      s << "  TOTAL:    " << (total>> 20) << "Mb" << "\n";
-    }
-
-    /** \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_matrix.nonZeros()>0) // otherwise L==I
-        derived().matrixL().solveInPlace(dest);
-
-      if(m_diag.size()>0)
-        dest = m_diag.asDiagonal().inverse() * dest;
-
-      if (m_matrix.nonZeros()>0) // otherwise I==I
-        derived().matrixU().solveInPlace(dest);
-
-      if(m_P.size()>0)
-        dest = m_P * dest;
-    }
-
-    /** \internal */
-    template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
-    void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &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());
-      
-      // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix.
-      static const int NbColsAtOnce = 4;
-      int rhsCols = b.cols();
-      int size = b.rows();
-      Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols);
-      for(int k=0; k<rhsCols; k+=NbColsAtOnce)
-      {
-        int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce);
-        tmp.leftCols(actualCols) = b.middleCols(k,actualCols);
-        tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols));
-        dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView();
-      }
-    }
-
-#endif // EIGEN_PARSED_BY_DOXYGEN
-
-  protected:
-
-    template<bool DoLDLt>
-    void factorize(const MatrixType& a);
-
-    void analyzePattern(const MatrixType& a, bool doLDLt);
-
-    /** keeps off-diagonal entries; drops diagonal entries */
-    struct keep_diag {
-      inline bool operator() (const Index& row, const Index& col, const Scalar&) const
-      {
-        return row!=col;
-      }
-    };
-
-    mutable ComputationInfo m_info;
-    bool m_isInitialized;
-    bool m_factorizationIsOk;
-    bool m_analysisIsOk;
-    
-    CholMatrixType m_matrix;
-    VectorType m_diag;                                // the diagonal coefficients (LDLt mode)
-    VectorXi m_parent;                                // elimination tree
-    VectorXi m_nonZerosPerCol;
-    PermutationMatrix<Dynamic,Dynamic,Index> m_P;     // the permutation
-    PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv;  // the inverse permutation
-};
-
-template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLt;
-template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLt;
-template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky;
-
-namespace internal {
-
-template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLt<_MatrixType,_UpLo> >
-{
-  typedef _MatrixType MatrixType;
-  enum { UpLo = _UpLo };
-  typedef typename MatrixType::Scalar                             Scalar;
-  typedef typename MatrixType::Index                              Index;
-  typedef SparseMatrix<Scalar, ColMajor, Index>          CholMatrixType;
-  typedef SparseTriangularView<CholMatrixType, Eigen::Lower>   MatrixL;
-  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper>   MatrixU;
-  inline static MatrixL getL(const MatrixType& m) { return m; }
-  inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
-};
-
-//template<typename _MatrixType> struct traits<SimplicialLLt<_MatrixType,Upper> >
-//{
-//  typedef _MatrixType MatrixType;
-//  enum { UpLo = Upper };
-//  typedef typename MatrixType::Scalar                             Scalar;
-//  typedef typename MatrixType::Index                              Index;
-//  typedef SparseMatrix<Scalar, ColMajor, Index>          CholMatrixType;
-//  typedef TriangularView<CholMatrixType, Eigen::Lower>   MatrixL;
-//  typedef TriangularView<CholMatrixType, Eigen::Upper>   MatrixU;
-//  inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
-//  inline static MatrixU getU(const MatrixType& m) { return m; }
-//};
-
-template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLt<_MatrixType,_UpLo> >
-{
-  typedef _MatrixType MatrixType;
-  enum { UpLo = _UpLo };
-  typedef typename MatrixType::Scalar                               Scalar;
-  typedef typename MatrixType::Index                                Index;
-  typedef SparseMatrix<Scalar, ColMajor, Index>            CholMatrixType;
-  typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
-  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
-  inline static MatrixL getL(const MatrixType& m) { return m; }
-  inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
-};
-
-//template<typename _MatrixType> struct traits<SimplicialLDLt<_MatrixType,Upper> >
-//{
-//  typedef _MatrixType MatrixType;
-//  enum { UpLo = Upper };
-//  typedef typename MatrixType::Scalar                               Scalar;
-//  typedef typename MatrixType::Index                                Index;
-//  typedef SparseMatrix<Scalar, ColMajor, Index>            CholMatrixType;
-//  typedef TriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
-//  typedef TriangularView<CholMatrixType, Eigen::UnitUpper> MatrixU;
-//  inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
-//  inline static MatrixU getU(const MatrixType& m) { return m; }
-//};
-
-template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
-{
-  typedef _MatrixType MatrixType;
-  enum { UpLo = _UpLo };
-};
-
-}
-
-/** \class SimplicialLLt
-  * \brief A direct sparse LLt Cholesky factorizations
-  *
-  * This class provides a LL^T Cholesky factorizations of sparse matrices that are
-  * selfadjoint and positive definite. The factorization allows for solving A.X = B where
-  * 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.
-  *
-  * \sa class SimplicialLDLt
-  */
-template<typename _MatrixType, int _UpLo>
-    class SimplicialLLt : public SimplicialCholeskyBase<SimplicialLLt<_MatrixType,_UpLo> >
-{
-public:
-    typedef _MatrixType MatrixType;
-    enum { UpLo = _UpLo };
-    typedef SimplicialCholeskyBase<SimplicialLLt> Base;
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
-    typedef typename MatrixType::Index Index;
-    typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
-    typedef Matrix<Scalar,Dynamic,1> VectorType;
-    typedef internal::traits<SimplicialLLt> Traits;
-    typedef typename Traits::MatrixL  MatrixL;
-    typedef typename Traits::MatrixU  MatrixU;
-public:
-    SimplicialLLt() : Base() {}
-    SimplicialLLt(const MatrixType& matrix)
-        : Base(matrix) {}
-
-    inline const MatrixL matrixL() const {
-        eigen_assert(Base::m_factorizationIsOk && "Simplicial LLt not factorized");
-        return Traits::getL(Base::m_matrix);
-    }
-
-    inline const MatrixU matrixU() const {
-        eigen_assert(Base::m_factorizationIsOk && "Simplicial LLt not factorized");
-        return Traits::getU(Base::m_matrix);
-    }
-
-    /** 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)
-    {
-      Base::analyzePattern(a, 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)
-    {
-      Base::template factorize<false>(a);
-    }
-
-    Scalar determinant() const
-    {
-      Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
-      return internal::abs2(detL);
-    }
-};
-
-/** \class SimplicialLDLt
-  * \brief A direct sparse LDLt Cholesky factorizations without square root.
-  *
-  * This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are
-  * selfadjoint and positive definite. The factorization allows for solving A.X = B where
-  * 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.
-  *
-  * \sa class SimplicialLLt
-  */
-template<typename _MatrixType, int _UpLo>
-    class SimplicialLDLt : public SimplicialCholeskyBase<SimplicialLDLt<_MatrixType,_UpLo> >
-{
-public:
-    typedef _MatrixType MatrixType;
-    enum { UpLo = _UpLo };
-    typedef SimplicialCholeskyBase<SimplicialLDLt> Base;
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
-    typedef typename MatrixType::Index Index;
-    typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
-    typedef Matrix<Scalar,Dynamic,1> VectorType;
-    typedef internal::traits<SimplicialLDLt> Traits;
-    typedef typename Traits::MatrixL  MatrixL;
-    typedef typename Traits::MatrixU  MatrixU;
-public:
-    SimplicialLDLt() : Base() {}
-    SimplicialLDLt(const MatrixType& matrix)
-        : Base(matrix) {}
-
-    inline const VectorType vectorD() const {
-        eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLt not factorized");
-        return Base::m_diag;
-    }
-    inline const MatrixL matrixL() const {
-        eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLt not factorized");
-        return Traits::getL(Base::m_matrix);
-    }
-
-    inline const MatrixU matrixU() const {
-        eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLt not factorized");
-        return Traits::getU(Base::m_matrix);
-    }
-
-    /** 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)
-    {
-      Base::analyzePattern(a, true);
-    }
-
-    /** 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)
-    {
-      Base::template factorize<true>(a);
-    }
-
-    Scalar determinant() const
-    {
-      return Base::m_diag.prod();
-    }
-};
-
-/** \class SimplicialCholesky
-  * \deprecated
-  * \sa class SimplicialLDLt, class SimplicialLLt
-  */
-template<typename _MatrixType, int _UpLo>
-    class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> >
-{
-public:
-    typedef _MatrixType MatrixType;
-    enum { UpLo = _UpLo };
-    typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
-    typedef typename MatrixType::Index Index;
-    typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
-    typedef Matrix<Scalar,Dynamic,1> VectorType;
-    typedef internal::traits<SimplicialCholesky> Traits;
-    typedef internal::traits<SimplicialLDLt<MatrixType,UpLo> > LDLtTraits;
-    typedef internal::traits<SimplicialLLt<MatrixType,UpLo>  > LLtTraits;
-  public:
-    SimplicialCholesky() : Base(), m_LDLt(true) {}
-    SimplicialCholesky(const MatrixType& matrix)
-      : Base(), m_LDLt(true)
-    {
-      Base::compute(matrix);
-    }
-
-    SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
-    {
-      switch(mode)
-      {
-      case SimplicialCholeskyLLt:
-        m_LDLt = false;
-        break;
-      case SimplicialCholeskyLDLt:
-        m_LDLt = true;
-        break;
-      default:
-        break;
-      }
-
-      return *this;
-    }
-
-    inline const VectorType vectorD() const {
-        eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
-        return Base::m_diag;
-    }
-    inline const CholMatrixType rawMatrix() const {
-        eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
-        return Base::m_matrix;
-    }
-
-    /** 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)
-    {
-      Base::analyzePattern(a, m_LDLt);
-    }
-
-    /** 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)
-    {
-      if(m_LDLt)
-        Base::template factorize<true>(a);
-      else
-        Base::template factorize<false>(a);
-    }
-
-    /** \internal */
-    template<typename Rhs,typename Dest>
-    void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
-    {
-      eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
-      eigen_assert(Base::m_matrix.rows()==b.rows());
-
-      if(Base::m_info!=Success)
-        return;
-
-      if(Base::m_P.size()>0)
-        dest = Base::m_Pinv * b;
-      else
-        dest = b;
-
-      if(Base::m_matrix.nonZeros()>0) // otherwise L==I
-      {
-        if(m_LDLt)
-          LDLtTraits::getL(Base::m_matrix).solveInPlace(dest);
-        else
-          LLtTraits::getL(Base::m_matrix).solveInPlace(dest);
-      }
-
-      if(Base::m_diag.size()>0)
-        dest = Base::m_diag.asDiagonal().inverse() * dest;
-
-      if (Base::m_matrix.nonZeros()>0) // otherwise I==I
-      {
-        if(m_LDLt)
-          LDLtTraits::getU(Base::m_matrix).solveInPlace(dest);
-        else
-          LLtTraits::getU(Base::m_matrix).solveInPlace(dest);
-      }
-
-      if(Base::m_P.size()>0)
-        dest = Base::m_P * dest;
-    }
-    
-    Scalar determinant() const
-    {
-      if(m_LDLt)
-      {
-        return Base::m_diag.prod();
-      }
-      else
-      {
-        Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
-        return internal::abs2(detL);
-      }
-    }
-    
-  protected:
-    bool m_LDLt;
-};
-
-template<typename Derived>
-void SimplicialCholeskyBase<Derived>::analyzePattern(const MatrixType& a, bool doLDLt)
-{
-  eigen_assert(a.rows()==a.cols());
-  const Index size = a.rows();
-  m_matrix.resize(size, size);
-  m_parent.resize(size);
-  m_nonZerosPerCol.resize(size);
-  
-  ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
-  
-  // TODO allows to configure the permutation
-  {
-    CholMatrixType C;
-    C = a.template selfadjointView<UpLo>();
-    // remove diagonal entries:
-    C.prune(keep_diag());
-    internal::minimum_degree_ordering(C, m_P);
-  }
-  
-  if(m_P.size()>0)
-    m_Pinv  = m_P.inverse();
-  else
-    m_Pinv.resize(0);
-  
-  SparseMatrix<Scalar,ColMajor,Index> ap(size,size);
-  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
-  
-  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 */
-    for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it)
-    {
-      Index i = 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 */
-        }
-      }
-    }
-  }
-  
-  /* 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] + (doLDLt ? 0 : 1);
-
-  m_matrix.resizeNonZeros(Lp[size]);
-  
-  m_isInitialized     = true;
-  m_info              = Success;
-  m_analysisIsOk      = true;
-  m_factorizationIsOk = false;
-}
-
-
-template<typename Derived>
-template<bool DoLDLt>
-void SimplicialCholeskyBase<Derived>::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();
-
-  ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
-  ei_declare_aligned_stack_constructed_variable(Index,  pattern, size, 0);
-  ei_declare_aligned_stack_constructed_variable(Index,  tags, size, 0);
-
-  SparseMatrix<Scalar,ColMajor,Index> ap(size,size);
-  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
-  
-  bool ok = true;
-  m_diag.resize(DoLDLt ? 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
-    for(typename MatrixType::InnerIterator it(ap,k); it; ++it)
-    {
-      Index i = 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)
-    {
-      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(DoLDLt)
-        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] + (DoLDLt ? 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(DoLDLt)
-      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;
-    }
-  }
-
-  m_info = ok ? Success : NumericalIssue;
-  m_factorizationIsOk = true;
-}
-
-namespace internal {
-  
-template<typename Derived, typename Rhs>
-struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
-  : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
-{
-  typedef SimplicialCholeskyBase<Derived> Dec;
-  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
-
-  template<typename Dest> void evalTo(Dest& dst) const
-  {
-    dec().derived()._solve(rhs(),dst);
-  }
-};
-
-template<typename Derived, typename Rhs>
-struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
-  : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
-{
-  typedef SimplicialCholeskyBase<Derived> Dec;
-  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
-
-  template<typename Dest> void evalTo(Dest& dst) const
-  {
-    dec().derived()._solve_sparse(rhs(),dst);
-  }
-};
-
-}
-
-#endif // EIGEN_SIMPLICIAL_CHOLESKY_H
diff --git a/unsupported/Eigen/src/SparseExtra/Solve.h b/unsupported/Eigen/src/SparseExtra/Solve.h
deleted file mode 100644
index 5b6c859ae..000000000
--- a/unsupported/Eigen/src/SparseExtra/Solve.h
+++ /dev/null
@@ -1,122 +0,0 @@
-// 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) {}
-
-
-
-template<typename DecompositionType, typename Rhs, typename Guess> struct solve_retval_with_guess;
-
-template<typename DecompositionType, typename Rhs, typename Guess>
-struct traits<solve_retval_with_guess<DecompositionType, Rhs, Guess> >
-{
-  typedef typename DecompositionType::MatrixType MatrixType;
-  typedef Matrix<typename Rhs::Scalar,
-                 MatrixType::ColsAtCompileTime,
-                 Rhs::ColsAtCompileTime,
-                 Rhs::PlainObject::Options,
-                 MatrixType::MaxColsAtCompileTime,
-                 Rhs::MaxColsAtCompileTime> ReturnType;
-};
-
-template<typename DecompositionType, typename Rhs, typename Guess> struct solve_retval_with_guess
- : public ReturnByValue<solve_retval_with_guess<DecompositionType, Rhs, Guess> >
-{
-  typedef typename DecompositionType::Index Index;
-
-  solve_retval_with_guess(const DecompositionType& dec, const Rhs& rhs, const Guess& guess)
-    : m_dec(dec), m_rhs(rhs), m_guess(guess)
-  {}
-
-  inline Index rows() const { return m_dec.cols(); }
-  inline Index cols() const { return m_rhs.cols(); }
-
-  template<typename Dest> inline void evalTo(Dest& dst) const
-  {
-    dst = m_guess;
-    m_dec._solveWithGuess(m_rhs,dst);
-  }
-
-  protected:
-    const DecompositionType& m_dec;
-    const typename Rhs::Nested m_rhs;
-    const typename Guess::Nested m_guess;
-};
-
-} // namepsace internal
-
-#endif // EIGEN_SPARSE_SOLVE_H
diff --git a/unsupported/Eigen/src/SparseExtra/SuperLUSupport.h b/unsupported/Eigen/src/SparseExtra/SuperLUSupport.h
deleted file mode 100644
index e485a9f50..000000000
--- a/unsupported/Eigen/src/SparseExtra/SuperLUSupport.h
+++ /dev/null
@@ -1,989 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008-2011 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_SUPERLUSUPPORT_H
-#define EIGEN_SUPERLUSUPPORT_H
-
-#define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE)                                                              \
-    extern "C" {                                                                                          \
-      typedef struct { FLOATTYPE for_lu; FLOATTYPE total_needed; int expansions; } PREFIX##mem_usage_t;   \
-      extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,                  \
-                                char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,           \
-                                void *, int, SuperMatrix *, SuperMatrix *,                                \
-                                FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *,                       \
-                                PREFIX##mem_usage_t *, SuperLUStat_t *, int *);                           \
-    }                                                                                                     \
-    inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,                                \
-         int *perm_c, int *perm_r, int *etree, char *equed,                                               \
-         FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                                      \
-         SuperMatrix *U, void *work, int lwork,                                                           \
-         SuperMatrix *B, SuperMatrix *X,                                                                  \
-         FLOATTYPE *recip_pivot_growth,                                                                   \
-         FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr,                                              \
-         SuperLUStat_t *stats, int *info, KEYTYPE) {                                                      \
-    PREFIX##mem_usage_t mem_usage;                                                                        \
-    PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L,                                      \
-         U, work, lwork, B, X, recip_pivot_growth, rcond,                                                 \
-         ferr, berr, &mem_usage, stats, info);                                                            \
-    return mem_usage.for_lu; /* bytes used by the factor storage */                                       \
-  }
-
-DECL_GSSVX(s,float,float)
-DECL_GSSVX(c,float,std::complex<float>)
-DECL_GSSVX(d,double,double)
-DECL_GSSVX(z,double,std::complex<double>)
-
-#ifdef MILU_ALPHA
-#define EIGEN_SUPERLU_HAS_ILU
-#endif
-
-#ifdef EIGEN_SUPERLU_HAS_ILU
-
-// similarly for the incomplete factorization using gsisx
-#define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE)                                                    \
-    extern "C" {                                                                                \
-      extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *,        \
-                         char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,        \
-                         void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *,   \
-                         PREFIX##mem_usage_t *, SuperLUStat_t *, int *);                        \
-    }                                                                                           \
-    inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A,                      \
-         int *perm_c, int *perm_r, int *etree, char *equed,                                     \
-         FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                            \
-         SuperMatrix *U, void *work, int lwork,                                                 \
-         SuperMatrix *B, SuperMatrix *X,                                                        \
-         FLOATTYPE *recip_pivot_growth,                                                         \
-         FLOATTYPE *rcond,                                                                      \
-         SuperLUStat_t *stats, int *info, KEYTYPE) {                                            \
-    PREFIX##mem_usage_t mem_usage;                                                              \
-    PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L,                            \
-         U, work, lwork, B, X, recip_pivot_growth, rcond,                                       \
-         &mem_usage, stats, info);                                                              \
-    return mem_usage.for_lu; /* bytes used by the factor storage */                             \
-  }
-
-DECL_GSISX(s,float,float)
-DECL_GSISX(c,float,std::complex<float>)
-DECL_GSISX(d,double,double)
-DECL_GSISX(z,double,std::complex<double>)
-
-#endif
-
-template<typename MatrixType>
-struct SluMatrixMapHelper;
-
-/** \internal
-  *
-  * A wrapper class for SuperLU matrices. It supports only compressed sparse matrices
-  * and dense matrices. Supernodal and other fancy format are not supported by this wrapper.
-  *
-  * This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure.
-  */
-struct SluMatrix : SuperMatrix
-{
-  SluMatrix()
-  {
-    Store = &storage;
-  }
-
-  SluMatrix(const SluMatrix& other)
-    : SuperMatrix(other)
-  {
-    Store = &storage;
-    storage = other.storage;
-  }
-
-  SluMatrix& operator=(const SluMatrix& other)
-  {
-    SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
-    Store = &storage;
-    storage = other.storage;
-    return *this;
-  }
-
-  struct
-  {
-    union {int nnz;int lda;};
-    void *values;
-    int *innerInd;
-    int *outerInd;
-  } storage;
-
-  void setStorageType(Stype_t t)
-  {
-    Stype = t;
-    if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
-      Store = &storage;
-    else
-    {
-      eigen_assert(false && "storage type not supported");
-      Store = 0;
-    }
-  }
-
-  template<typename Scalar>
-  void setScalarType()
-  {
-    if (internal::is_same<Scalar,float>::value)
-      Dtype = SLU_S;
-    else if (internal::is_same<Scalar,double>::value)
-      Dtype = SLU_D;
-    else if (internal::is_same<Scalar,std::complex<float> >::value)
-      Dtype = SLU_C;
-    else if (internal::is_same<Scalar,std::complex<double> >::value)
-      Dtype = SLU_Z;
-    else
-    {
-      eigen_assert(false && "Scalar type not supported by SuperLU");
-    }
-  }
-
-  template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
-  static SluMatrix Map(Matrix<Scalar,Rows,Cols,Options,MRows,MCols>& mat)
-  {
-    typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
-    eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
-    SluMatrix res;
-    res.setStorageType(SLU_DN);
-    res.setScalarType<Scalar>();
-    res.Mtype     = SLU_GE;
-
-    res.nrow      = mat.rows();
-    res.ncol      = mat.cols();
-
-    res.storage.lda       = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride();
-    res.storage.values    = mat.data();
-    return res;
-  }
-
-  template<typename MatrixType>
-  static SluMatrix Map(SparseMatrixBase<MatrixType>& mat)
-  {
-    SluMatrix res;
-    if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
-    {
-      res.setStorageType(SLU_NR);
-      res.nrow      = mat.cols();
-      res.ncol      = mat.rows();
-    }
-    else
-    {
-      res.setStorageType(SLU_NC);
-      res.nrow      = mat.rows();
-      res.ncol      = mat.cols();
-    }
-
-    res.Mtype       = SLU_GE;
-
-    res.storage.nnz       = mat.nonZeros();
-    res.storage.values    = mat.derived()._valuePtr();
-    res.storage.innerInd  = mat.derived()._innerIndexPtr();
-    res.storage.outerInd  = mat.derived()._outerIndexPtr();
-
-    res.setScalarType<typename MatrixType::Scalar>();
-
-    // FIXME the following is not very accurate
-    if (MatrixType::Flags & Upper)
-      res.Mtype = SLU_TRU;
-    if (MatrixType::Flags & Lower)
-      res.Mtype = SLU_TRL;
-
-    eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
-
-    return res;
-  }
-};
-
-template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
-struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
-{
-  typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
-  static void run(MatrixType& mat, SluMatrix& res)
-  {
-    eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
-    res.setStorageType(SLU_DN);
-    res.setScalarType<Scalar>();
-    res.Mtype     = SLU_GE;
-
-    res.nrow      = mat.rows();
-    res.ncol      = mat.cols();
-
-    res.storage.lda       = mat.outerStride();
-    res.storage.values    = mat.data();
-  }
-};
-
-template<typename Derived>
-struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
-{
-  typedef Derived MatrixType;
-  static void run(MatrixType& mat, SluMatrix& res)
-  {
-    if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
-    {
-      res.setStorageType(SLU_NR);
-      res.nrow      = mat.cols();
-      res.ncol      = mat.rows();
-    }
-    else
-    {
-      res.setStorageType(SLU_NC);
-      res.nrow      = mat.rows();
-      res.ncol      = mat.cols();
-    }
-
-    res.Mtype       = SLU_GE;
-
-    res.storage.nnz       = mat.nonZeros();
-    res.storage.values    = mat._valuePtr();
-    res.storage.innerInd  = mat._innerIndexPtr();
-    res.storage.outerInd  = mat._outerIndexPtr();
-
-    res.setScalarType<typename MatrixType::Scalar>();
-
-    // FIXME the following is not very accurate
-    if (MatrixType::Flags & Upper)
-      res.Mtype = SLU_TRU;
-    if (MatrixType::Flags & Lower)
-      res.Mtype = SLU_TRL;
-
-    eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
-  }
-};
-
-namespace internal {
-
-template<typename MatrixType>
-SluMatrix asSluMatrix(MatrixType& mat)
-{
-  return SluMatrix::Map(mat);
-}
-
-/** View a Super LU matrix as an Eigen expression */
-template<typename Scalar, int Flags, typename Index>
-MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
-{
-  eigen_assert((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR
-         || (Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC);
-
-  Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
-
-  return MappedSparseMatrix<Scalar,Flags,Index>(
-    sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
-    sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
-}
-
-} // end namespace internal
-
-
-template<typename _MatrixType, typename Derived>
-class SuperLUBase
-{
-  public:
-    typedef _MatrixType MatrixType;
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
-    typedef typename MatrixType::Index Index;
-    typedef Matrix<Scalar,Dynamic,1> Vector;
-    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
-    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;    
-    typedef SparseMatrix<Scalar> LUMatrixType;
-
-  public:
-
-    SuperLUBase() {}
-
-    ~SuperLUBase()
-    {
-      clearFactors();
-    }
-    
-    Derived& derived() { return *static_cast<Derived*>(this); }
-    const Derived& derived() const { return *static_cast<const Derived*>(this); }
-    
-    inline Index rows() const { return m_matrix.rows(); }
-    inline Index cols() const { return m_matrix.cols(); }
-    
-    /** \returns a reference to the Super LU option object to configure the  Super LU algorithms. */
-    inline superlu_options_t& options() { return m_sluOptions; }
-    
-    /** \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 */
-    void compute(const MatrixType& matrix)
-    {
-      derived().analyzePattern(matrix);
-      derived().factorize(matrix);
-    }
-    
-    /** \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<SuperLUBase, Rhs> solve(const MatrixBase<Rhs>& b) const
-    {
-      eigen_assert(m_isInitialized && "SuperLU is not initialized.");
-      eigen_assert(rows()==b.rows()
-                && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
-      return internal::solve_retval<SuperLUBase, 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<SuperLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
-//     {
-//       eigen_assert(m_isInitialized && "SuperLU is not initialized.");
-//       eigen_assert(rows()==b.rows()
-//                 && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
-//       return internal::sparse_solve_retval<SuperLU, 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& /*matrix*/)
-    {
-      m_isInitialized = true;
-      m_info = Success;
-      m_analysisIsOk = true;
-      m_factorizationIsOk = false;
-    }
-    
-    template<typename Stream>
-    void dumpMemory(Stream& s)
-    {}
-    
-  protected:
-    
-    void initFactorization(const MatrixType& a)
-    {
-      const int size = a.rows();
-      m_matrix = a;
-
-      m_sluA = internal::asSluMatrix(m_matrix);
-      clearFactors();
-
-      m_p.resize(size);
-      m_q.resize(size);
-      m_sluRscale.resize(size);
-      m_sluCscale.resize(size);
-      m_sluEtree.resize(size);
-
-      // set empty B and X
-      m_sluB.setStorageType(SLU_DN);
-      m_sluB.setScalarType<Scalar>();
-      m_sluB.Mtype          = SLU_GE;
-      m_sluB.storage.values = 0;
-      m_sluB.nrow           = 0;
-      m_sluB.ncol           = 0;
-      m_sluB.storage.lda    = size;
-      m_sluX                = m_sluB;
-      
-      m_extractedDataAreDirty = true;
-    }
-    
-    void init()
-    {
-      m_info = InvalidInput;
-      m_isInitialized = false;
-      m_sluL.Store = 0;
-      m_sluU.Store = 0;
-    }
-    
-    void extractData() const;
-
-    void clearFactors()
-    {
-      if(m_sluL.Store)
-        Destroy_SuperNode_Matrix(&m_sluL);
-      if(m_sluU.Store)
-        Destroy_CompCol_Matrix(&m_sluU);
-
-      m_sluL.Store = 0;
-      m_sluU.Store = 0;
-
-      memset(&m_sluL,0,sizeof m_sluL);
-      memset(&m_sluU,0,sizeof m_sluU);
-    }
-
-    // cached data to reduce reallocation, etc.
-    mutable LUMatrixType m_l;
-    mutable LUMatrixType m_u;
-    mutable IntColVectorType m_p;
-    mutable IntRowVectorType m_q;
-
-    mutable LUMatrixType m_matrix;  // copy of the factorized matrix
-    mutable SluMatrix m_sluA;
-    mutable SuperMatrix m_sluL, m_sluU;
-    mutable SluMatrix m_sluB, m_sluX;
-    mutable SuperLUStat_t m_sluStat;
-    mutable superlu_options_t m_sluOptions;
-    mutable std::vector<int> m_sluEtree;
-    mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale;
-    mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr;
-    mutable char m_sluEqued;
-
-    mutable ComputationInfo m_info;
-    bool m_isInitialized;
-    int m_factorizationIsOk;
-    int m_analysisIsOk;
-    mutable bool m_extractedDataAreDirty;
-};
-
-
-template<typename _MatrixType>
-class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
-{
-  public:
-    typedef SuperLUBase<_MatrixType,SuperLU> Base;
-    typedef _MatrixType MatrixType;
-    typedef typename Base::Scalar Scalar;
-    typedef typename Base::RealScalar RealScalar;
-    typedef typename Base::Index Index;
-    typedef typename Base::IntRowVectorType IntRowVectorType;
-    typedef typename Base::IntColVectorType IntColVectorType;    
-    typedef typename Base::LUMatrixType LUMatrixType;
-    typedef TriangularView<LUMatrixType, Lower|UnitDiag>  LMatrixType;
-    typedef TriangularView<LUMatrixType,  Upper>           UMatrixType;
-
-  public:
-
-    SuperLU() : Base() { init(); }
-
-    SuperLU(const MatrixType& matrix) : Base()
-    {
-      Base::init();
-      compute(matrix);
-    }
-
-    ~SuperLU()
-    {
-    }
-    
-    /** 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& matrix)
-    {
-      init();
-      Base::analyzePattern(matrix);
-    }
-    
-    /** 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& matrix);
-    
-    #ifndef EIGEN_PARSED_BY_DOXYGEN
-    /** \internal */
-    template<typename Rhs,typename Dest>
-    void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
-    #endif // EIGEN_PARSED_BY_DOXYGEN
-    
-    inline const LMatrixType& matrixL() const
-    {
-      if (m_extractedDataAreDirty) this->extractData();
-      return m_l;
-    }
-
-    inline const UMatrixType& matrixU() const
-    {
-      if (m_extractedDataAreDirty) this->extractData();
-      return m_u;
-    }
-
-    inline const IntColVectorType& permutationP() const
-    {
-      if (m_extractedDataAreDirty) this->extractData();
-      return m_p;
-    }
-
-    inline const IntRowVectorType& permutationQ() const
-    {
-      if (m_extractedDataAreDirty) this->extractData();
-      return m_q;
-    }
-    
-    Scalar determinant() const;
-    
-  protected:
-    
-    using Base::m_matrix;
-    using Base::m_sluOptions;
-    using Base::m_sluA;
-    using Base::m_sluB;
-    using Base::m_sluX;
-    using Base::m_p;
-    using Base::m_q;
-    using Base::m_sluEtree;
-    using Base::m_sluEqued;
-    using Base::m_sluRscale;
-    using Base::m_sluCscale;
-    using Base::m_sluL;
-    using Base::m_sluU;
-    using Base::m_sluStat;
-    using Base::m_sluFerr;
-    using Base::m_sluBerr;
-    using Base::m_l;
-    using Base::m_u;
-    
-    using Base::m_analysisIsOk;
-    using Base::m_factorizationIsOk;
-    using Base::m_extractedDataAreDirty;
-    using Base::m_isInitialized;
-    using Base::m_info;
-    
-    void init()
-    {
-      Base::init();
-      
-      set_default_options(&this->m_sluOptions);
-      m_sluOptions.PrintStat        = NO;
-      m_sluOptions.ConditionNumber  = NO;
-      m_sluOptions.Trans            = NOTRANS;
-      m_sluOptions.ColPerm          = COLAMD;
-    }
-};
-
-template<typename MatrixType>
-void SuperLU<MatrixType>::factorize(const MatrixType& a)
-{
-  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
-  if(!m_analysisIsOk)
-  {
-    m_info = InvalidInput;
-    return;
-  }
-  
-  initFactorization(a);
-  
-  int info = 0;
-  RealScalar recip_pivot_growth, rcond;
-  RealScalar ferr, berr;
-
-  StatInit(&m_sluStat);
-  SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
-                &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
-                &m_sluL, &m_sluU,
-                NULL, 0,
-                &m_sluB, &m_sluX,
-                &recip_pivot_growth, &rcond,
-                &ferr, &berr,
-                &m_sluStat, &info, Scalar());
-  StatFree(&m_sluStat);
-
-  m_extractedDataAreDirty = true;
-
-  // FIXME how to better check for errors ???
-  m_info = info == 0 ? Success : NumericalIssue;
-  m_factorizationIsOk = true;
-}
-
-template<typename MatrixType>
-template<typename Rhs,typename Dest>
-void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
-{
-  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
-
-  const int size = m_matrix.rows();
-  const int rhsCols = b.cols();
-  eigen_assert(size==b.rows());
-
-  m_sluOptions.Trans = NOTRANS;
-  m_sluOptions.Fact = FACTORED;
-  m_sluOptions.IterRefine = NOREFINE;
-  
-
-  m_sluFerr.resize(rhsCols);
-  m_sluBerr.resize(rhsCols);
-  m_sluB = SluMatrix::Map(b.const_cast_derived());
-  m_sluX = SluMatrix::Map(x.derived());
-  
-  typename Rhs::PlainObject b_cpy;
-  if(m_sluEqued!='N')
-  {
-    b_cpy = b;
-    m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());  
-  }
-
-  StatInit(&m_sluStat);
-  int info = 0;
-  RealScalar recip_pivot_growth, rcond;
-  SuperLU_gssvx(&m_sluOptions, &m_sluA,
-                m_q.data(), m_p.data(),
-                &m_sluEtree[0], &m_sluEqued,
-                &m_sluRscale[0], &m_sluCscale[0],
-                &m_sluL, &m_sluU,
-                NULL, 0,
-                &m_sluB, &m_sluX,
-                &recip_pivot_growth, &rcond,
-                &m_sluFerr[0], &m_sluBerr[0],
-                &m_sluStat, &info, Scalar());
-  StatFree(&m_sluStat);
-  m_info = info==0 ? Success : NumericalIssue;
-}
-
-// the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
-//
-//  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
-//
-//  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
-//  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
-//
-template<typename MatrixType, typename Derived>
-void SuperLUBase<MatrixType,Derived>::extractData() const
-{
-  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
-  if (m_extractedDataAreDirty)
-  {
-    int         upper;
-    int         fsupc, istart, nsupr;
-    int         lastl = 0, lastu = 0;
-    SCformat    *Lstore = static_cast<SCformat*>(m_sluL.Store);
-    NCformat    *Ustore = static_cast<NCformat*>(m_sluU.Store);
-    Scalar      *SNptr;
-
-    const int size = m_matrix.rows();
-    m_l.resize(size,size);
-    m_l.resizeNonZeros(Lstore->nnz);
-    m_u.resize(size,size);
-    m_u.resizeNonZeros(Ustore->nnz);
-
-    int* Lcol = m_l._outerIndexPtr();
-    int* Lrow = m_l._innerIndexPtr();
-    Scalar* Lval = m_l._valuePtr();
-
-    int* Ucol = m_u._outerIndexPtr();
-    int* Urow = m_u._innerIndexPtr();
-    Scalar* Uval = m_u._valuePtr();
-
-    Ucol[0] = 0;
-    Ucol[0] = 0;
-
-    /* for each supernode */
-    for (int k = 0; k <= Lstore->nsuper; ++k)
-    {
-      fsupc   = L_FST_SUPC(k);
-      istart  = L_SUB_START(fsupc);
-      nsupr   = L_SUB_START(fsupc+1) - istart;
-      upper   = 1;
-
-      /* for each column in the supernode */
-      for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
-      {
-        SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
-
-        /* Extract U */
-        for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
-        {
-          Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
-          /* Matlab doesn't like explicit zero. */
-          if (Uval[lastu] != 0.0)
-            Urow[lastu++] = U_SUB(i);
-        }
-        for (int i = 0; i < upper; ++i)
-        {
-          /* upper triangle in the supernode */
-          Uval[lastu] = SNptr[i];
-          /* Matlab doesn't like explicit zero. */
-          if (Uval[lastu] != 0.0)
-            Urow[lastu++] = L_SUB(istart+i);
-        }
-        Ucol[j+1] = lastu;
-
-        /* Extract L */
-        Lval[lastl] = 1.0; /* unit diagonal */
-        Lrow[lastl++] = L_SUB(istart + upper - 1);
-        for (int i = upper; i < nsupr; ++i)
-        {
-          Lval[lastl] = SNptr[i];
-          /* Matlab doesn't like explicit zero. */
-          if (Lval[lastl] != 0.0)
-            Lrow[lastl++] = L_SUB(istart+i);
-        }
-        Lcol[j+1] = lastl;
-
-        ++upper;
-      } /* for j ... */
-
-    } /* for k ... */
-
-    // squeeze the matrices :
-    m_l.resizeNonZeros(lastl);
-    m_u.resizeNonZeros(lastu);
-
-    m_extractedDataAreDirty = false;
-  }
-}
-
-template<typename MatrixType>
-typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
-{
-  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for computing the determinant, you must first call either compute() or analyzePattern()/factorize()");
-  
-  if (m_extractedDataAreDirty)
-    this->extractData();
-
-  Scalar det = Scalar(1);
-  for (int j=0; j<m_u.cols(); ++j)
-  {
-    if (m_u._outerIndexPtr()[j+1]-m_u._outerIndexPtr()[j] > 0)
-    {
-      int lastId = m_u._outerIndexPtr()[j+1]-1;
-      eigen_assert(m_u._innerIndexPtr()[lastId]<=j);
-      if (m_u._innerIndexPtr()[lastId]==j)
-        det *= m_u._valuePtr()[lastId];
-    }
-  }
-  if(m_sluEqued!='N')
-    return det/m_sluRscale.prod()/m_sluCscale.prod();
-  else
-    return det;
-}
-
-#ifdef EIGEN_SUPERLU_HAS_ILU
-template<typename _MatrixType>
-class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
-{
-  public:
-    typedef SuperLUBase<_MatrixType,SuperILU> Base;
-    typedef _MatrixType MatrixType;
-    typedef typename Base::Scalar Scalar;
-    typedef typename Base::RealScalar RealScalar;
-    typedef typename Base::Index Index;
-
-  public:
-
-    SuperILU() : Base() { init(); }
-
-    SuperILU(const MatrixType& matrix) : Base()
-    {
-      init();
-      compute(matrix);
-    }
-
-    ~SuperILU()
-    {
-    }
-    
-    /** 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& matrix)
-    {
-      Base::analyzePattern(matrix);
-    }
-    
-    /** 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& matrix);
-    
-    #ifndef EIGEN_PARSED_BY_DOXYGEN
-    /** \internal */
-    template<typename Rhs,typename Dest>
-    void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
-    #endif // EIGEN_PARSED_BY_DOXYGEN
-    
-  protected:
-    
-    using Base::m_matrix;
-    using Base::m_sluOptions;
-    using Base::m_sluA;
-    using Base::m_sluB;
-    using Base::m_sluX;
-    using Base::m_p;
-    using Base::m_q;
-    using Base::m_sluEtree;
-    using Base::m_sluEqued;
-    using Base::m_sluRscale;
-    using Base::m_sluCscale;
-    using Base::m_sluL;
-    using Base::m_sluU;
-    using Base::m_sluStat;
-    using Base::m_sluFerr;
-    using Base::m_sluBerr;
-    using Base::m_l;
-    using Base::m_u;
-    
-    using Base::m_analysisIsOk;
-    using Base::m_factorizationIsOk;
-    using Base::m_extractedDataAreDirty;
-    using Base::m_isInitialized;
-    using Base::m_info;
-
-    void init()
-    {
-      Base::init();
-      
-      ilu_set_default_options(&m_sluOptions);
-      m_sluOptions.PrintStat        = NO;
-      m_sluOptions.ConditionNumber  = NO;
-      m_sluOptions.Trans            = NOTRANS;
-      m_sluOptions.ColPerm          = MMD_AT_PLUS_A;
-      
-      // no attempt to preserve column sum
-      m_sluOptions.ILU_MILU = SILU;
-      // only basic ILU(k) support -- no direct control over memory consumption
-      // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
-      // and set ILU_FillFactor to max memory growth
-      m_sluOptions.ILU_DropRule = DROP_BASIC;
-      m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
-    }
-};
-
-template<typename MatrixType>
-void SuperILU<MatrixType>::factorize(const MatrixType& a)
-{
-  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
-  if(!m_analysisIsOk)
-  {
-    m_info = InvalidInput;
-    return;
-  }
-  
-  this->initFactorization(a);
-
-  int info = 0;
-  RealScalar recip_pivot_growth, rcond;
-
-  StatInit(&m_sluStat);
-  SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
-                &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
-                &m_sluL, &m_sluU,
-                NULL, 0,
-                &m_sluB, &m_sluX,
-                &recip_pivot_growth, &rcond,
-                &m_sluStat, &info, Scalar());
-  StatFree(&m_sluStat);
-
-  // FIXME how to better check for errors ???
-  m_info = info == 0 ? Success : NumericalIssue;
-  m_factorizationIsOk = true;
-}
-
-template<typename MatrixType>
-template<typename Rhs,typename Dest>
-void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
-{
-  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
-
-  const int size = m_matrix.rows();
-  const int rhsCols = b.cols();
-  eigen_assert(size==b.rows());
-
-  m_sluOptions.Trans = NOTRANS;
-  m_sluOptions.Fact = FACTORED;
-  m_sluOptions.IterRefine = NOREFINE;
-
-  m_sluFerr.resize(rhsCols);
-  m_sluBerr.resize(rhsCols);
-  m_sluB = SluMatrix::Map(b.const_cast_derived());
-  m_sluX = SluMatrix::Map(x.derived());
-
-  typename Rhs::PlainObject b_cpy;
-  if(m_sluEqued!='N')
-  {
-    b_cpy = b;
-    m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());  
-  }
-  
-  int info = 0;
-  RealScalar recip_pivot_growth, rcond;
-
-  StatInit(&m_sluStat);
-  SuperLU_gsisx(&m_sluOptions, &m_sluA,
-                m_q.data(), m_p.data(),
-                &m_sluEtree[0], &m_sluEqued,
-                &m_sluRscale[0], &m_sluCscale[0],
-                &m_sluL, &m_sluU,
-                NULL, 0,
-                &m_sluB, &m_sluX,
-                &recip_pivot_growth, &rcond,
-                &m_sluStat, &info, Scalar());
-  StatFree(&m_sluStat);
-
-  m_info = info==0 ? Success : NumericalIssue;
-}
-#endif
-
-namespace internal {
-  
-template<typename _MatrixType, typename Derived, typename Rhs>
-struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
-  : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
-{
-  typedef SuperLUBase<_MatrixType,Derived> Dec;
-  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
-
-  template<typename Dest> void evalTo(Dest& dst) const
-  {
-    dec().derived()._solve(rhs(),dst);
-  }
-};
-
-template<typename _MatrixType, typename Derived, typename Rhs>
-struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
-  : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
-{
-  typedef SuperLUBase<_MatrixType,Derived> Dec;
-  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
-
-  template<typename Dest> void evalTo(Dest& dst) const
-  {
-    dec().derived()._solve(rhs(),dst);
-  }
-};
-
-}
-
-#endif // EIGEN_SUPERLUSUPPORT_H
diff --git a/unsupported/Eigen/src/SparseExtra/SuperLUSupportLegacy.h b/unsupported/Eigen/src/SparseExtra/SuperLUSupportLegacy.h
deleted file mode 100644
index e85d8d36e..000000000
--- a/unsupported/Eigen/src/SparseExtra/SuperLUSupportLegacy.h
+++ /dev/null
@@ -1,407 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008-2009 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_SUPERLUSUPPORT_LEGACY_H
-#define EIGEN_SUPERLUSUPPORT_LEGACY_H
-
-/** \deprecated use class BiCGSTAB, class SuperLU, or class UmfPackLU */
-template<typename MatrixType>
-class SparseLU<MatrixType,SuperLULegacy> : public SparseLU<MatrixType>
-{
-  protected:
-    typedef SparseLU<MatrixType> Base;
-    typedef typename Base::Scalar Scalar;
-    typedef typename Base::RealScalar RealScalar;
-    typedef Matrix<Scalar,Dynamic,1> Vector;
-    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
-    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
-    typedef SparseMatrix<Scalar,Lower|UnitDiag> LMatrixType;
-    typedef SparseMatrix<Scalar,Upper> UMatrixType;
-    using Base::m_flags;
-    using Base::m_status;
-
-  public:
-
-    /** \deprecated the entire class is deprecated */
-  EIGEN_DEPRECATED SparseLU(int flags = NaturalOrdering)
-      : Base(flags)
-    {
-    }
-
-    /** \deprecated the entire class is deprecated */
-    EIGEN_DEPRECATED SparseLU(const MatrixType& matrix, int flags = NaturalOrdering)
-      : Base(flags)
-    {
-      compute(matrix);
-    }
-
-    ~SparseLU()
-    {
-      Destroy_SuperNode_Matrix(&m_sluL);
-      Destroy_CompCol_Matrix(&m_sluU);
-    }
-
-    inline const LMatrixType& matrixL() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_l;
-    }
-
-    inline const UMatrixType& matrixU() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_u;
-    }
-
-    inline const IntColVectorType& permutationP() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_p;
-    }
-
-    inline const IntRowVectorType& permutationQ() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_q;
-    }
-
-    Scalar determinant() const;
-
-    template<typename BDerived, typename XDerived>
-    bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed = SvNoTrans) const;
-
-    void compute(const MatrixType& matrix);
-
-  protected:
-
-    void extractData() const;
-
-  protected:
-    // cached data to reduce reallocation, etc.
-    mutable LMatrixType m_l;
-    mutable UMatrixType m_u;
-    mutable IntColVectorType m_p;
-    mutable IntRowVectorType m_q;
-
-    mutable SparseMatrix<Scalar> m_matrix;
-    mutable SluMatrix m_sluA;
-    mutable SuperMatrix m_sluL, m_sluU;
-    mutable SluMatrix m_sluB, m_sluX;
-    mutable SuperLUStat_t m_sluStat;
-    mutable superlu_options_t m_sluOptions;
-    mutable std::vector<int> m_sluEtree;
-    mutable std::vector<RealScalar> m_sluRscale, m_sluCscale;
-    mutable std::vector<RealScalar> m_sluFerr, m_sluBerr;
-    mutable char m_sluEqued;
-    mutable bool m_extractedDataAreDirty;
-};
-
-template<typename MatrixType>
-void SparseLU<MatrixType,SuperLULegacy>::compute(const MatrixType& a)
-{
-  const int size = a.rows();
-  m_matrix = a;
-
-  set_default_options(&m_sluOptions);
-  m_sluOptions.ColPerm = NATURAL;
-  m_sluOptions.PrintStat = NO;
-  m_sluOptions.ConditionNumber = NO;
-  m_sluOptions.Trans = NOTRANS;
-  // m_sluOptions.Equil = NO;
-
-  switch (Base::orderingMethod())
-  {
-      case NaturalOrdering          : m_sluOptions.ColPerm = NATURAL; break;
-      case MinimumDegree_AT_PLUS_A  : m_sluOptions.ColPerm = MMD_AT_PLUS_A; break;
-      case MinimumDegree_ATA        : m_sluOptions.ColPerm = MMD_ATA; break;
-      case ColApproxMinimumDegree   : m_sluOptions.ColPerm = COLAMD; break;
-      default:
-        //std::cerr << "Eigen: ordering method \"" << Base::orderingMethod() << "\" not supported by the SuperLU backend\n";
-        m_sluOptions.ColPerm = NATURAL;
-  };
-
-  m_sluA = internal::asSluMatrix(m_matrix);
-  memset(&m_sluL,0,sizeof m_sluL);
-  memset(&m_sluU,0,sizeof m_sluU);
-  m_sluEqued = 'N';
-  int info = 0;
-
-  m_p.resize(size);
-  m_q.resize(size);
-  m_sluRscale.resize(size);
-  m_sluCscale.resize(size);
-  m_sluEtree.resize(size);
-
-  RealScalar recip_pivot_gross, rcond;
-  RealScalar ferr, berr;
-
-  // set empty B and X
-  m_sluB.setStorageType(SLU_DN);
-  m_sluB.setScalarType<Scalar>();
-  m_sluB.Mtype = SLU_GE;
-  m_sluB.storage.values = 0;
-  m_sluB.nrow = m_sluB.ncol = 0;
-  m_sluB.storage.lda = size;
-  m_sluX = m_sluB;
-
-  StatInit(&m_sluStat);
-  if (m_flags&IncompleteFactorization)
-  {
-    #ifdef EIGEN_SUPERLU_HAS_ILU
-    ilu_set_default_options(&m_sluOptions);
-
-    // no attempt to preserve column sum
-    m_sluOptions.ILU_MILU = SILU;
-
-    // only basic ILU(k) support -- no direct control over memory consumption
-    // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
-    // and set ILU_FillFactor to max memory growth
-    m_sluOptions.ILU_DropRule = DROP_BASIC;
-    m_sluOptions.ILU_DropTol = Base::m_precision;
-
-    SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
-      &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
-      &m_sluL, &m_sluU,
-      NULL, 0,
-      &m_sluB, &m_sluX,
-      &recip_pivot_gross, &rcond,
-      &m_sluStat, &info, Scalar());
-    #else
-    //std::cerr << "Incomplete factorization is only available in SuperLU v4\n";
-    Base::m_succeeded = false;
-    return;
-    #endif
-  }
-  else
-  {
-    SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
-      &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
-      &m_sluL, &m_sluU,
-      NULL, 0,
-      &m_sluB, &m_sluX,
-      &recip_pivot_gross, &rcond,
-      &ferr, &berr,
-      &m_sluStat, &info, Scalar());
-  }
-  StatFree(&m_sluStat);
-
-  m_extractedDataAreDirty = true;
-
-  // FIXME how to better check for errors ???
-  Base::m_succeeded = (info == 0);
-}
-
-template<typename MatrixType>
-template<typename BDerived,typename XDerived>
-bool SparseLU<MatrixType,SuperLULegacy>::solve(const MatrixBase<BDerived> &b,
-                        MatrixBase<XDerived> *x, const int transposed) const
-{
-  const int size = m_matrix.rows();
-  const int rhsCols = b.cols();
-  eigen_assert(size==b.rows());
-
-  switch (transposed) {
-      case SvNoTrans    :  m_sluOptions.Trans = NOTRANS; break;
-      case SvTranspose  :  m_sluOptions.Trans = TRANS;   break;
-      case SvAdjoint    :  m_sluOptions.Trans = CONJ;    break;
-      default:
-        //std::cerr << "Eigen: transposition  option \"" << transposed << "\" not supported by the SuperLU backend\n";
-        m_sluOptions.Trans = NOTRANS;
-  }
-
-  m_sluOptions.Fact = FACTORED;
-  m_sluOptions.IterRefine = NOREFINE;
-
-  m_sluFerr.resize(rhsCols);
-  m_sluBerr.resize(rhsCols);
-  m_sluB = SluMatrix::Map(b.const_cast_derived());
-  m_sluX = SluMatrix::Map(x->derived());
-  
-  typename BDerived::PlainObject b_cpy;
-  if(m_sluEqued!='N')
-  {
-    b_cpy = b;
-    m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());  
-  }
-
-  StatInit(&m_sluStat);
-  int info = 0;
-  RealScalar recip_pivot_gross, rcond;
-
-  if (m_flags&IncompleteFactorization)
-  {
-    #ifdef EIGEN_SUPERLU_HAS_ILU
-    SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
-      &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
-      &m_sluL, &m_sluU,
-      NULL, 0,
-      &m_sluB, &m_sluX,
-      &recip_pivot_gross, &rcond,
-      &m_sluStat, &info, Scalar());
-    #else
-    //std::cerr << "Incomplete factorization is only available in SuperLU v4\n";
-    return false;
-    #endif
-  }
-  else
-  {
-    SuperLU_gssvx(
-      &m_sluOptions, &m_sluA,
-      m_q.data(), m_p.data(),
-      &m_sluEtree[0], &m_sluEqued,
-      &m_sluRscale[0], &m_sluCscale[0],
-      &m_sluL, &m_sluU,
-      NULL, 0,
-      &m_sluB, &m_sluX,
-      &recip_pivot_gross, &rcond,
-      &m_sluFerr[0], &m_sluBerr[0],
-      &m_sluStat, &info, Scalar());
-  }
-  StatFree(&m_sluStat);
-
-  // reset to previous state
-  m_sluOptions.Trans = NOTRANS;
-  return info==0;
-}
-
-//
-// the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
-//
-//  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
-//
-//  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
-//  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
-//
-template<typename MatrixType>
-void SparseLU<MatrixType,SuperLULegacy>::extractData() const
-{
-  if (m_extractedDataAreDirty)
-  {
-    int         upper;
-    int         fsupc, istart, nsupr;
-    int         lastl = 0, lastu = 0;
-    SCformat    *Lstore = static_cast<SCformat*>(m_sluL.Store);
-    NCformat    *Ustore = static_cast<NCformat*>(m_sluU.Store);
-    Scalar      *SNptr;
-
-    const int size = m_matrix.rows();
-    m_l.resize(size,size);
-    m_l.resizeNonZeros(Lstore->nnz);
-    m_u.resize(size,size);
-    m_u.resizeNonZeros(Ustore->nnz);
-
-    int* Lcol = m_l._outerIndexPtr();
-    int* Lrow = m_l._innerIndexPtr();
-    Scalar* Lval = m_l._valuePtr();
-
-    int* Ucol = m_u._outerIndexPtr();
-    int* Urow = m_u._innerIndexPtr();
-    Scalar* Uval = m_u._valuePtr();
-
-    Ucol[0] = 0;
-    Ucol[0] = 0;
-
-    /* for each supernode */
-    for (int k = 0; k <= Lstore->nsuper; ++k)
-    {
-      fsupc   = L_FST_SUPC(k);
-      istart  = L_SUB_START(fsupc);
-      nsupr   = L_SUB_START(fsupc+1) - istart;
-      upper   = 1;
-
-      /* for each column in the supernode */
-      for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
-      {
-        SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
-
-        /* Extract U */
-        for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
-        {
-          Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
-          /* Matlab doesn't like explicit zero. */
-          if (Uval[lastu] != 0.0)
-            Urow[lastu++] = U_SUB(i);
-        }
-        for (int i = 0; i < upper; ++i)
-        {
-          /* upper triangle in the supernode */
-          Uval[lastu] = SNptr[i];
-          /* Matlab doesn't like explicit zero. */
-          if (Uval[lastu] != 0.0)
-            Urow[lastu++] = L_SUB(istart+i);
-        }
-        Ucol[j+1] = lastu;
-
-        /* Extract L */
-        Lval[lastl] = 1.0; /* unit diagonal */
-        Lrow[lastl++] = L_SUB(istart + upper - 1);
-        for (int i = upper; i < nsupr; ++i)
-        {
-          Lval[lastl] = SNptr[i];
-          /* Matlab doesn't like explicit zero. */
-          if (Lval[lastl] != 0.0)
-            Lrow[lastl++] = L_SUB(istart+i);
-        }
-        Lcol[j+1] = lastl;
-
-        ++upper;
-      } /* for j ... */
-
-    } /* for k ... */
-
-    // squeeze the matrices :
-    m_l.resizeNonZeros(lastl);
-    m_u.resizeNonZeros(lastu);
-
-    m_extractedDataAreDirty = false;
-  }
-}
-
-template<typename MatrixType>
-typename SparseLU<MatrixType,SuperLULegacy>::Scalar SparseLU<MatrixType,SuperLULegacy>::determinant() const
-{
-  assert((!NumTraits<Scalar>::IsComplex) && "This function is not implemented for complex yet");
-  if (m_extractedDataAreDirty)
-    extractData();
-
-  // TODO this code could be moved to the default/base backend
-  // FIXME perhaps we have to take into account the scale factors m_sluRscale and m_sluCscale ???
-  Scalar det = Scalar(1);
-  for (int j=0; j<m_u.cols(); ++j)
-  {
-    if (m_u._outerIndexPtr()[j+1]-m_u._outerIndexPtr()[j] > 0)
-    {
-      int lastId = m_u._outerIndexPtr()[j+1]-1;
-      eigen_assert(m_u._innerIndexPtr()[lastId]<=j);
-      if (m_u._innerIndexPtr()[lastId]==j)
-      {
-        det *= m_u._valuePtr()[lastId];
-      }
-    }
-//       std::cout << m_sluRscale[j] << " " << m_sluCscale[j] << "   \n";
-  }
-  return det;
-}
-
-#endif // EIGEN_SUPERLUSUPPORT_LEGACY_H
diff --git a/unsupported/Eigen/src/SparseExtra/UmfPackSupport.h b/unsupported/Eigen/src/SparseExtra/UmfPackSupport.h
deleted file mode 100644
index e41de8337..000000000
--- a/unsupported/Eigen/src/SparseExtra/UmfPackSupport.h
+++ /dev/null
@@ -1,406 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008-2011 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_UMFPACKSUPPORT_H
-#define EIGEN_UMFPACKSUPPORT_H
-
-/* TODO extract L, extract U, compute det, etc... */
-
-// generic double/complex<double> wrapper functions:
-
-inline void umfpack_free_numeric(void **Numeric, double)
-{ umfpack_di_free_numeric(Numeric); *Numeric = 0; }
-
-inline void umfpack_free_numeric(void **Numeric, std::complex<double>)
-{ umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
-
-inline void umfpack_free_symbolic(void **Symbolic, double)
-{ umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
-
-inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>)
-{ umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
-
-inline int umfpack_symbolic(int n_row,int n_col,
-                            const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
-                            const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
-{
-  return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
-}
-
-inline int umfpack_symbolic(int n_row,int n_col,
-                            const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic,
-                            const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
-{
-  return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&internal::real_ref(Ax[0]),0,Symbolic,Control,Info);
-}
-
-inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[],
-                            void *Symbolic, void **Numeric,
-                            const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
-{
-  return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
-}
-
-inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[],
-                            void *Symbolic, void **Numeric,
-                            const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
-{
-  return umfpack_zi_numeric(Ap,Ai,&internal::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
-}
-
-inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[],
-                          double X[], const double B[], void *Numeric,
-                          const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
-{
-  return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
-}
-
-inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[],
-                          std::complex<double> X[], const std::complex<double> B[], void *Numeric,
-                          const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
-{
-  return umfpack_zi_solve(sys,Ap,Ai,&internal::real_ref(Ax[0]),0,&internal::real_ref(X[0]),0,&internal::real_ref(B[0]),0,Numeric,Control,Info);
-}
-
-inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
-{
-  return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
-}
-
-inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>)
-{
-  return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
-}
-
-inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[],
-                               int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
-{
-  return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
-}
-
-inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[],
-                               int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric)
-{
-  double& lx0_real = internal::real_ref(Lx[0]);
-  double& ux0_real = internal::real_ref(Ux[0]);
-  double& dx0_real = internal::real_ref(Dx[0]);
-  return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
-                                Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
-}
-
-inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
-{
-  return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
-}
-
-inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
-{
-  double& mx_real = internal::real_ref(*Mx);
-  return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
-}
-
-/** \brief A sparse LU factorization and solver based on UmfPack
-  *
-  * This class allows to solve for A.X = B sparse linear problems via a LU factorization
-  * using the UmfPack library. The sparse matrix A must be column-major, squared and full rank.
-  * 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<>
-  *
-  */
-template<typename _MatrixType>
-class UmfPackLU
-{
-  public:
-    typedef _MatrixType MatrixType;
-    typedef typename MatrixType::Scalar Scalar;
-    typedef typename MatrixType::RealScalar RealScalar;
-    typedef typename MatrixType::Index Index;
-    typedef Matrix<Scalar,Dynamic,1> Vector;
-    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
-    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
-    typedef SparseMatrix<Scalar> LUMatrixType;
-
-  public:
-
-    UmfPackLU() { init(); }
-
-    UmfPackLU(const MatrixType& matrix)
-    {
-      init();
-      compute(matrix);
-    }
-
-    ~UmfPackLU()
-    {
-      if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
-      if(m_numeric)  umfpack_free_numeric(&m_numeric,Scalar());
-    }
-
-    inline Index rows() const { return m_matrixRef->rows(); }
-    inline Index cols() const { return m_matrixRef->cols(); }
-
-    /** \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;
-    }
-
-    inline const LUMatrixType& matrixL() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_l;
-    }
-
-    inline const LUMatrixType& matrixU() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_u;
-    }
-
-    inline const IntColVectorType& permutationP() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_p;
-    }
-
-    inline const IntRowVectorType& permutationQ() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_q;
-    }
-
-    /** Computes the sparse Cholesky decomposition of \a matrix */
-    void compute(const MatrixType& matrix)
-    {
-      analyzePattern(matrix);
-      factorize(matrix);
-    }
-
-    /** \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<UmfPackLU, Rhs> solve(const MatrixBase<Rhs>& b) const
-    {
-      eigen_assert(m_isInitialized && "UmfPAckLU is not initialized.");
-      eigen_assert(rows()==b.rows()
-                && "UmfPAckLU::solve(): invalid number of rows of the right hand side matrix b");
-      return internal::solve_retval<UmfPackLU, 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<UmfPAckLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
-//     {
-//       eigen_assert(m_isInitialized && "UmfPAckLU is not initialized.");
-//       eigen_assert(rows()==b.rows()
-//                 && "UmfPAckLU::solve(): invalid number of rows of the right hand side matrix b");
-//       return internal::sparse_solve_retval<UmfPAckLU, 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& matrix)
-    {
-      eigen_assert((MatrixType::Flags&RowMajorBit)==0 && "UmfPackLU: Row major matrices are not supported yet");
-
-      if(m_symbolic)
-        umfpack_free_symbolic(&m_symbolic,Scalar());
-      if(m_numeric)
-        umfpack_free_numeric(&m_numeric,Scalar());
-
-      int errorCode = 0;
-      errorCode = umfpack_symbolic(matrix.rows(), matrix.cols(), matrix._outerIndexPtr(), matrix._innerIndexPtr(), matrix._valuePtr(),
-                                   &m_symbolic, 0, 0);
-
-      m_isInitialized = true;
-      m_info = errorCode ? InvalidInput : 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& matrix)
-    {
-      eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
-      if(m_numeric)
-        umfpack_free_numeric(&m_numeric,Scalar());
-
-      m_matrixRef = &matrix;
-
-      int errorCode;
-      errorCode = umfpack_numeric(matrix._outerIndexPtr(), matrix._innerIndexPtr(), matrix._valuePtr(),
-                                  m_symbolic, &m_numeric, 0, 0);
-
-      m_info = errorCode ? NumericalIssue : Success;
-      m_factorizationIsOk = true;
-    }
-
-    #ifndef EIGEN_PARSED_BY_DOXYGEN
-    /** \internal */
-    template<typename BDerived,typename XDerived>
-    bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
-    #endif
-
-    Scalar determinant() const;
-
-    void extractData() const;
-
-  protected:
-
-
-    void init()
-    {
-      m_info = InvalidInput;
-      m_isInitialized = false;
-      m_numeric = 0;
-      m_symbolic = 0;
-    }
-
-    // cached data to reduce reallocation, etc.
-    mutable LUMatrixType m_l;
-    mutable LUMatrixType m_u;
-    mutable IntColVectorType m_p;
-    mutable IntRowVectorType m_q;
-
-    const MatrixType* m_matrixRef;
-    void* m_numeric;
-    void* m_symbolic;
-
-    mutable ComputationInfo m_info;
-    bool m_isInitialized;
-    int m_factorizationIsOk;
-    int m_analysisIsOk;
-    mutable bool m_extractedDataAreDirty;
-};
-
-
-template<typename MatrixType>
-void UmfPackLU<MatrixType>::extractData() const
-{
-  if (m_extractedDataAreDirty)
-  {
-    // get size of the data
-    int lnz, unz, rows, cols, nz_udiag;
-    umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
-
-    // allocate data
-    m_l.resize(rows,(std::min)(rows,cols));
-    m_l.resizeNonZeros(lnz);
-
-    m_u.resize((std::min)(rows,cols),cols);
-    m_u.resizeNonZeros(unz);
-
-    m_p.resize(rows);
-    m_q.resize(cols);
-
-    // extract
-    umfpack_get_numeric(m_l._outerIndexPtr(), m_l._innerIndexPtr(), m_l._valuePtr(),
-                        m_u._outerIndexPtr(), m_u._innerIndexPtr(), m_u._valuePtr(),
-                        m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
-
-    m_extractedDataAreDirty = false;
-  }
-}
-
-template<typename MatrixType>
-typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const
-{
-  Scalar det;
-  umfpack_get_determinant(&det, 0, m_numeric, 0);
-  return det;
-}
-
-template<typename MatrixType>
-template<typename BDerived,typename XDerived>
-bool UmfPackLU<MatrixType>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const
-{
-  const int rhsCols = b.cols();
-  eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet");
-  eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet");
-
-  int errorCode;
-  for (int j=0; j<rhsCols; ++j)
-  {
-    errorCode = umfpack_solve(UMFPACK_A,
-        m_matrixRef->_outerIndexPtr(), m_matrixRef->_innerIndexPtr(), m_matrixRef->_valuePtr(),
-        &x.col(j).coeffRef(0), &b.const_cast_derived().col(j).coeffRef(0), m_numeric, 0, 0);
-    if (errorCode!=0)
-      return false;
-  }
-
-  return true;
-}
-
-
-namespace internal {
-
-template<typename _MatrixType, typename Rhs>
-struct solve_retval<UmfPackLU<_MatrixType>, Rhs>
-  : solve_retval_base<UmfPackLU<_MatrixType>, Rhs>
-{
-  typedef UmfPackLU<_MatrixType> Dec;
-  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
-
-  template<typename Dest> void evalTo(Dest& dst) const
-  {
-    dec()._solve(rhs(),dst);
-  }
-};
-
-template<typename _MatrixType, typename Rhs>
-struct sparse_solve_retval<UmfPackLU<_MatrixType>, Rhs>
-  : sparse_solve_retval_base<UmfPackLU<_MatrixType>, Rhs>
-{
-  typedef UmfPackLU<_MatrixType> Dec;
-  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
-
-  template<typename Dest> void evalTo(Dest& dst) const
-  {
-    dec()._solve(rhs(),dst);
-  }
-};
-
-}
-
-#endif // EIGEN_UMFPACKSUPPORT_H
diff --git a/unsupported/Eigen/src/SparseExtra/UmfPackSupportLegacy.h b/unsupported/Eigen/src/SparseExtra/UmfPackSupportLegacy.h
deleted file mode 100644
index 3d30e1ed1..000000000
--- a/unsupported/Eigen/src/SparseExtra/UmfPackSupportLegacy.h
+++ /dev/null
@@ -1,257 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008-2009 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_UMFPACKSUPPORT_LEGACY_H
-#define EIGEN_UMFPACKSUPPORT_LEGACY_H
-
-/** \deprecated use class BiCGSTAB, class SuperLU, or class UmfPackLU */
-template<typename _MatrixType>
-class SparseLU<_MatrixType,UmfPack> : public SparseLU<_MatrixType>
-{
-  protected:
-    typedef SparseLU<_MatrixType> Base;
-    typedef typename Base::Scalar Scalar;
-    typedef typename Base::RealScalar RealScalar;
-    typedef Matrix<Scalar,Dynamic,1> Vector;
-    typedef Matrix<int, 1, _MatrixType::ColsAtCompileTime> IntRowVectorType;
-    typedef Matrix<int, _MatrixType::RowsAtCompileTime, 1> IntColVectorType;
-    typedef SparseMatrix<Scalar,Lower|UnitDiag> LMatrixType;
-    typedef SparseMatrix<Scalar,Upper> UMatrixType;
-    using Base::m_flags;
-    using Base::m_status;
-
-  public:
-    typedef _MatrixType MatrixType;
-    typedef typename MatrixType::Index Index;
-
-    /** \deprecated the entire class is deprecated */
-    EIGEN_DEPRECATED SparseLU(int flags = NaturalOrdering)
-      : Base(flags), m_numeric(0)
-    {
-    }
-
-    /** \deprecated the entire class is deprecated */
-    EIGEN_DEPRECATED SparseLU(const MatrixType& matrix, int flags = NaturalOrdering)
-      : Base(flags), m_numeric(0)
-    {
-      compute(matrix);
-    }
-
-    ~SparseLU()
-    {
-      if (m_numeric)
-        umfpack_free_numeric(&m_numeric,Scalar());
-    }
-
-    inline const LMatrixType& matrixL() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_l;
-    }
-
-    inline const UMatrixType& matrixU() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_u;
-    }
-
-    inline const IntColVectorType& permutationP() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_p;
-    }
-
-    inline const IntRowVectorType& permutationQ() const
-    {
-      if (m_extractedDataAreDirty) extractData();
-      return m_q;
-    }
-
-    Scalar determinant() const;
-
-    template<typename BDerived, typename XDerived>
-    bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x) const;
-
-    template<typename Rhs>
-      inline const internal::solve_retval<SparseLU<MatrixType, UmfPack>, Rhs>
-    solve(const MatrixBase<Rhs>& b) const
-    {
-      eigen_assert(true && "SparseLU is not initialized.");
-      return internal::solve_retval<SparseLU<MatrixType, UmfPack>, Rhs>(*this, b.derived());
-    }
-
-    void compute(const MatrixType& matrix);
-
-    inline Index cols() const { return m_matrixRef->cols(); }
-    inline Index rows() const { return m_matrixRef->rows(); }
-
-    inline const MatrixType& matrixLU() const
-    {
-      //eigen_assert(m_isInitialized && "LU is not initialized.");
-      return *m_matrixRef;
-    }
-
-    const void* numeric() const
-    {
-      return m_numeric;
-    }
-
-  protected:
-
-    void extractData() const;
-  
-  protected:
-    // cached data:
-    void* m_numeric;
-    const MatrixType* m_matrixRef;
-    mutable LMatrixType m_l;
-    mutable UMatrixType m_u;
-    mutable IntColVectorType m_p;
-    mutable IntRowVectorType m_q;
-    mutable bool m_extractedDataAreDirty;
-};
-
-namespace internal {
-
-template<typename _MatrixType, typename Rhs>
-  struct solve_retval<SparseLU<_MatrixType, UmfPack>, Rhs>
-  : solve_retval_base<SparseLU<_MatrixType, UmfPack>, Rhs>
-{
-  typedef SparseLU<_MatrixType, UmfPack> SpLUDecType;
-  EIGEN_MAKE_SOLVE_HELPERS(SpLUDecType,Rhs)
-
-  template<typename Dest> void evalTo(Dest& dst) const
-  {
-    const int rhsCols = rhs().cols();
-
-    eigen_assert((Rhs::Flags&RowMajorBit)==0 && "UmfPack backend does not support non col-major rhs yet");
-    eigen_assert((Dest::Flags&RowMajorBit)==0 && "UmfPack backend does not support non col-major result yet");
-
-    void* numeric = const_cast<void*>(dec().numeric());
-
-    EIGEN_UNUSED int errorCode = 0;
-    for (int j=0; j<rhsCols; ++j)
-    {
-      errorCode = umfpack_solve(UMFPACK_A,
-                                dec().matrixLU()._outerIndexPtr(), dec().matrixLU()._innerIndexPtr(), dec().matrixLU()._valuePtr(),
-                                &dst.col(j).coeffRef(0), &rhs().const_cast_derived().col(j).coeffRef(0), numeric, 0, 0);
-      eigen_assert(!errorCode && "UmfPack could not solve the system.");
-    }
-  }
-    
-};
-
-} // end namespace internal
-
-template<typename MatrixType>
-void SparseLU<MatrixType,UmfPack>::compute(const MatrixType& a)
-{
-  typedef typename MatrixType::Index Index;
-  const Index rows = a.rows();
-  const Index cols = a.cols();
-  eigen_assert((MatrixType::Flags&RowMajorBit)==0 && "Row major matrices are not supported yet");
-
-  m_matrixRef = &a;
-
-  if (m_numeric)
-    umfpack_free_numeric(&m_numeric,Scalar());
-
-  void* symbolic;
-  int errorCode = 0;
-  errorCode = umfpack_symbolic(rows, cols, a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(),
-                                  &symbolic, 0, 0);
-  if (errorCode==0)
-    errorCode = umfpack_numeric(a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(),
-                                   symbolic, &m_numeric, 0, 0);
-
-  umfpack_free_symbolic(&symbolic,Scalar());
-
-  m_extractedDataAreDirty = true;
-
-  Base::m_succeeded = (errorCode==0);
-}
-
-template<typename MatrixType>
-void SparseLU<MatrixType,UmfPack>::extractData() const
-{
-  if (m_extractedDataAreDirty)
-  {
-    // get size of the data
-    int lnz, unz, rows, cols, nz_udiag;
-    umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
-
-    // allocate data
-    m_l.resize(rows,(std::min)(rows,cols));
-    m_l.resizeNonZeros(lnz);
-    
-    m_u.resize((std::min)(rows,cols),cols);
-    m_u.resizeNonZeros(unz);
-
-    m_p.resize(rows);
-    m_q.resize(cols);
-
-    // extract
-    umfpack_get_numeric(m_l._outerIndexPtr(), m_l._innerIndexPtr(), m_l._valuePtr(),
-                        m_u._outerIndexPtr(), m_u._innerIndexPtr(), m_u._valuePtr(),
-                        m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
-    
-    m_extractedDataAreDirty = false;
-  }
-}
-
-template<typename MatrixType>
-typename SparseLU<MatrixType,UmfPack>::Scalar SparseLU<MatrixType,UmfPack>::determinant() const
-{
-  Scalar det;
-  umfpack_get_determinant(&det, 0, m_numeric, 0);
-  return det;
-}
-
-template<typename MatrixType>
-template<typename BDerived,typename XDerived>
-bool SparseLU<MatrixType,UmfPack>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> *x) const
-{
-  //const int size = m_matrix.rows();
-  const int rhsCols = b.cols();
-//   eigen_assert(size==b.rows());
-  eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPack backend does not support non col-major rhs yet");
-  eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPack backend does not support non col-major result yet");
-
-  int errorCode;
-  for (int j=0; j<rhsCols; ++j)
-  {
-    errorCode = umfpack_solve(UMFPACK_A,
-        m_matrixRef->_outerIndexPtr(), m_matrixRef->_innerIndexPtr(), m_matrixRef->_valuePtr(),
-        &x->col(j).coeffRef(0), &b.const_cast_derived().col(j).coeffRef(0), m_numeric, 0, 0);
-    if (errorCode!=0)
-      return false;
-  }
-//   errorCode = umfpack_di_solve(UMFPACK_A,
-//       m_matrixRef._outerIndexPtr(), m_matrixRef._innerIndexPtr(), m_matrixRef._valuePtr(),
-//       x->derived().data(), b.derived().data(), m_numeric, 0, 0);
-
-  return true;
-}
-
-#endif // EIGEN_UMFPACKSUPPORT_H