* sparse LU: add extraction of L,U,P, and Q, as well as determinant

             for both backends.
* extended a bit the sparse unit tests

5 files changed, 301 insertions, 47 deletions

diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
index 554d8bd3c..fda369499 100644
--- a/Eigen/src/LU/LU.h
+++ b/Eigen/src/LU/LU.h
@@ -190,10 +190,7 @@ template<typename MatrixType> class LU
       * \sa MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse()
       */
     template<typename OtherDerived, typename ResultType>
-    bool solve(
-    const MatrixBase<OtherDerived>& b,
-    ResultType *result
-    ) const;
+    bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
 
     /** \returns the determinant of the matrix of which
       * *this is the LU decomposition. It has only linear complexity
diff --git a/Eigen/src/Sparse/SparseLLT.h b/Eigen/src/Sparse/SparseLLT.h
index b7d4f5bbd..f9eb459dc 100644
--- a/Eigen/src/Sparse/SparseLLT.h
+++ b/Eigen/src/Sparse/SparseLLT.h
@@ -160,7 +160,7 @@ void SparseLLT<MatrixType,Backend>::compute(const MatrixType& a)
       {
         Scalar y = it.value();
         x -= ei_abs2(y);
-        ++it; // skip j-th element, and process remaing column coefficients
+        ++it; // skip j-th element, and process remaining column coefficients
         tempVector.restart();
         for (; it; ++it)
         {
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index 74e81f4bb..1346b310a 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -189,6 +189,10 @@ class SparseMatrix
         m_outerSize = outerSize;
       }
     }
+    void resizeNonZeros(int size)
+    {
+      m_data.resize(size);
+    }
 
     inline SparseMatrix()
       : m_outerSize(0), m_innerSize(0), m_outerIndex(0)
diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/Eigen/src/Sparse/SuperLUSupport.h
index 2ccfe2267..aeab835aa 100644
--- a/Eigen/src/Sparse/SuperLUSupport.h
+++ b/Eigen/src/Sparse/SuperLUSupport.h
@@ -211,6 +211,10 @@ class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType>
     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|UnitDiagBit> LMatrixType;
+    typedef SparseMatrix<Scalar,Upper> UMatrixType;
     using Base::m_flags;
     using Base::m_status;
 
@@ -231,23 +235,59 @@ class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType>
     {
     }
 
+    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;
 
     void compute(const MatrixType& matrix);
 
   protected:
-    // cached data to reduce reallocation:
+
+    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 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, m_sluPermR, m_sluPermC;
+    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>
@@ -261,6 +301,7 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
   m_sluOptions.PrintStat = NO;
   m_sluOptions.ConditionNumber = NO;
   m_sluOptions.Trans = NOTRANS;
+  // m_sluOptions.Equil = NO;
 
   switch (Base::orderingMethod())
   {
@@ -279,8 +320,8 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
   m_sluEqued = 'B';
   int info = 0;
 
-  m_sluPermR.resize(size);
-  m_sluPermC.resize(size);
+  m_p.resize(size);
+  m_q.resize(size);
   m_sluRscale.resize(size);
   m_sluCscale.resize(size);
   m_sluEtree.resize(size);
@@ -298,7 +339,7 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
   m_sluX = m_sluB;
 
   StatInit(&m_sluStat);
-  SuperLU_gssvx(&m_sluOptions, &m_sluA, &m_sluPermC[0], &m_sluPermR[0], &m_sluEtree[0],
+  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,
@@ -308,26 +349,12 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
           &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>
-// inline const MatrixType&
-// SparseLU<MatrixType,SuperLU>::matrixL() const
-// {
-//   ei_assert(false && "matrixL() is Not supported by the SuperLU backend");
-//   return m_matrix;
-// }
-//
-// template<typename MatrixType>
-// inline const MatrixType&
-// SparseLU<MatrixType,SuperLU>::matrixU() const
-// {
-//   ei_assert(false && "matrixU() is Not supported by the SuperLU backend");
-//   return m_matrix;
-// }
-
 template<typename MatrixType>
 template<typename BDerived,typename XDerived>
 bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> *x) const
@@ -349,7 +376,7 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBa
   RealScalar recip_pivot_gross, rcond;
   SuperLU_gssvx(
     &m_sluOptions, &m_sluA,
-    &m_sluPermC[0], &m_sluPermR[0],
+    m_q.data(), m_p.data(),
     &m_sluEtree[0], &m_sluEqued,
     &m_sluRscale[0], &m_sluCscale[0],
     &m_sluL, &m_sluU,
@@ -363,4 +390,122 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, MatrixBa
   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,SuperLU>::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,SuperLU>::Scalar SparseLU<MatrixType,SuperLU>::determinant() const
+{
+  if (m_extractedDataAreDirty)
+    extractData();
+
+  // TODO this code coule 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;
+      ei_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] << "   ";
+  }
+  return det;
+}
+
 #endif // EIGEN_SUPERLUSUPPORT_H
diff --git a/Eigen/src/Sparse/UmfPackSupport.h b/Eigen/src/Sparse/UmfPackSupport.h
index 751e390f6..00fa7e57d 100644
--- a/Eigen/src/Sparse/UmfPackSupport.h
+++ b/Eigen/src/Sparse/UmfPackSupport.h
@@ -25,7 +25,21 @@
 #ifndef EIGEN_UMFPACKSUPPORT_H
 #define EIGEN_UMFPACKSUPPORT_H
 
-/* TODO extract L, extrac U, compute det, etc... */
+/* 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); }
+
+inline void umfpack_free_numeric(void **Numeric, std::complex<double>)
+{ umfpack_zi_free_numeric(Numeric); }
+
+inline void umfpack_free_symbolic(void **Symbolic, double)
+{ umfpack_di_free_symbolic(Symbolic); }
+
+inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>)
+{ umfpack_zi_free_symbolic(Symbolic); }
 
 inline int umfpack_symbolic(int n_row,int n_col,
                             const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
@@ -69,6 +83,39 @@ inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::co
   return umfpack_zi_solve(sys,Ap,Ai,&Ax[0].real(),0,&X[0].real(),0,&B[0].real(),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)
+{
+  return umfpack_zi_get_numeric(Lp,Lj,Lx?&Lx[0].real():0,0,Up,Ui,Ux?&Ux[0].real():0,0,P,Q,
+                               Dx?&Dx[0].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])
+{
+  return umfpack_zi_get_determinant(&Mx->real(),0,Ex,NumericHandle,User_Info);
+}
+
 
 template<typename MatrixType>
 class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
@@ -78,6 +125,10 @@ class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
     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|UnitDiagBit> LMatrixType;
+    typedef SparseMatrix<Scalar,Upper> UMatrixType;
     using Base::m_flags;
     using Base::m_status;
 
@@ -97,18 +148,53 @@ class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
     ~SparseLU()
     {
       if (m_numeric)
-        umfpack_di_free_numeric(&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;
 
     void compute(const MatrixType& matrix);
 
   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;
 };
 
 template<typename MatrixType>
@@ -121,7 +207,7 @@ void SparseLU<MatrixType,UmfPack>::compute(const MatrixType& a)
   m_matrixRef = &a;
 
   if (m_numeric)
-    umfpack_di_free_numeric(&m_numeric);
+    umfpack_free_numeric(&m_numeric,Scalar());
 
   void* symbolic;
   int errorCode = 0;
@@ -131,26 +217,48 @@ void SparseLU<MatrixType,UmfPack>::compute(const MatrixType& a)
     errorCode = umfpack_numeric(a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(),
                                    symbolic, &m_numeric, 0, 0);
 
-  umfpack_di_free_symbolic(&symbolic);
+  umfpack_free_symbolic(&symbolic,Scalar());
+
+  m_extractedDataAreDirty = true;
 
   Base::m_succeeded = (errorCode==0);
 }
 
-// template<typename MatrixType>
-// inline const MatrixType&
-// SparseLU<MatrixType,SuperLU>::matrixL() const
-// {
-//   ei_assert(false && "matrixL() is Not supported by the SuperLU backend");
-//   return m_matrix;
-// }
-//
-// template<typename MatrixType>
-// inline const MatrixType&
-// SparseLU<MatrixType,SuperLU>::matrixU() const
-// {
-//   ei_assert(false && "matrixU() is Not supported by the SuperLU backend");
-//   return m_matrix;
-// }
+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>