add missing files

5 files changed, 433 insertions, 0 deletions

diff --git a/unsupported/Eigen/src/SparseExtra/CMakeLists.txt b/unsupported/Eigen/src/SparseExtra/CMakeLists.txt
new file mode 100644
index 000000000..7ea32ca5e
--- /dev/null
+++ b/unsupported/Eigen/src/SparseExtra/CMakeLists.txt
@@ -0,0 +1,6 @@
+FILE(GLOB Eigen_SparseExtra_SRCS "*.h")
+
+INSTALL(FILES
+  ${Eigen_SparseExtra_SRCS}
+  DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/SparseExtra COMPONENT Devel
+  )
diff --git a/unsupported/test/sparse_extra.cpp b/unsupported/test/sparse_extra.cpp
new file mode 100644
index 000000000..fa6dc50f5
--- /dev/null
+++ b/unsupported/test/sparse_extra.cpp
@@ -0,0 +1,153 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2010 Gael Guennebaud <g.gael@free.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/>.
+
+#include "sparse.h"
+#include <Eigen/SparseExtra>
+
+template<typename SetterType,typename DenseType, typename Scalar, int Options>
+bool test_random_setter(SparseMatrix<Scalar,Options>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords)
+{
+  typedef SparseMatrix<Scalar,Options> SparseType;
+  {
+    sm.setZero();
+    SetterType w(sm);
+    std::vector<Vector2i> remaining = nonzeroCoords;
+    while(!remaining.empty())
+    {
+      int i = ei_random<int>(0,static_cast<int>(remaining.size())-1);
+      w(remaining[i].x(),remaining[i].y()) = ref.coeff(remaining[i].x(),remaining[i].y());
+      remaining[i] = remaining.back();
+      remaining.pop_back();
+    }
+  }
+  return sm.isApprox(ref);
+}
+
+template<typename SetterType,typename DenseType, typename T>
+bool test_random_setter(DynamicSparseMatrix<T>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords)
+{
+  sm.setZero();
+  std::vector<Vector2i> remaining = nonzeroCoords;
+  while(!remaining.empty())
+  {
+    int i = ei_random<int>(0,static_cast<int>(remaining.size())-1);
+    sm.coeffRef(remaining[i].x(),remaining[i].y()) = ref.coeff(remaining[i].x(),remaining[i].y());
+    remaining[i] = remaining.back();
+    remaining.pop_back();
+  }
+  return sm.isApprox(ref);
+}
+
+template<typename SparseMatrixType> void sparse_extra(const SparseMatrixType& ref)
+{
+  const int rows = ref.rows();
+  const int cols = ref.cols();
+  typedef typename SparseMatrixType::Scalar Scalar;
+  enum { Flags = SparseMatrixType::Flags };
+
+  double density = std::max(8./(rows*cols), 0.01);
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+  typedef Matrix<Scalar,Dynamic,1> DenseVector;
+  Scalar eps = 1e-6;
+
+  SparseMatrixType m(rows, cols);
+  DenseMatrix refMat = DenseMatrix::Zero(rows, cols);
+  DenseVector vec1 = DenseVector::Random(rows);
+
+  std::vector<Vector2i> zeroCoords;
+  std::vector<Vector2i> nonzeroCoords;
+  initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords);
+
+  if (zeroCoords.size()==0 || nonzeroCoords.size()==0)
+    return;
+
+  // test coeff and coeffRef
+  for (int i=0; i<(int)zeroCoords.size(); ++i)
+  {
+    VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps );
+    if(ei_is_same_type<SparseMatrixType,SparseMatrix<Scalar,Flags> >::ret)
+      VERIFY_RAISES_ASSERT( m.coeffRef(zeroCoords[0].x(),zeroCoords[0].y()) = 5 );
+  }
+  VERIFY_IS_APPROX(m, refMat);
+
+  m.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
+  refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
+
+  VERIFY_IS_APPROX(m, refMat);
+
+  // random setter
+//   {
+//     m.setZero();
+//     VERIFY_IS_NOT_APPROX(m, refMat);
+//     SparseSetter<SparseMatrixType, RandomAccessPattern> w(m);
+//     std::vector<Vector2i> remaining = nonzeroCoords;
+//     while(!remaining.empty())
+//     {
+//       int i = ei_random<int>(0,remaining.size()-1);
+//       w->coeffRef(remaining[i].x(),remaining[i].y()) = refMat.coeff(remaining[i].x(),remaining[i].y());
+//       remaining[i] = remaining.back();
+//       remaining.pop_back();
+//     }
+//   }
+//   VERIFY_IS_APPROX(m, refMat);
+
+    VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, StdMapTraits> >(m,refMat,nonzeroCoords) ));
+    #ifdef EIGEN_UNORDERED_MAP_SUPPORT
+    VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, StdUnorderedMapTraits> >(m,refMat,nonzeroCoords) ));
+    #endif
+    #ifdef _DENSE_HASH_MAP_H_
+    VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleDenseHashMapTraits> >(m,refMat,nonzeroCoords) ));
+    #endif
+    #ifdef _SPARSE_HASH_MAP_H_
+    VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleSparseHashMapTraits> >(m,refMat,nonzeroCoords) ));
+    #endif
+
+
+  // test RandomSetter
+  /*{
+    SparseMatrixType m1(rows,cols), m2(rows,cols);
+    DenseMatrix refM1 = DenseMatrix::Zero(rows, rows);
+    initSparse<Scalar>(density, refM1, m1);
+    {
+      Eigen::RandomSetter<SparseMatrixType > setter(m2);
+      for (int j=0; j<m1.outerSize(); ++j)
+        for (typename SparseMatrixType::InnerIterator i(m1,j); i; ++i)
+          setter(i.index(), j) = i.value();
+    }
+    VERIFY_IS_APPROX(m1, m2);
+  }*/
+
+
+}
+
+void test_sparse_extra()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST_1( sparse_extra(SparseMatrix<double>(8, 8)) );
+    CALL_SUBTEST_2( sparse_extra(SparseMatrix<std::complex<double> >(16, 16)) );
+    CALL_SUBTEST_1( sparse_extra(SparseMatrix<double>(33, 33)) );
+
+    CALL_SUBTEST_3( sparse_extra(DynamicSparseMatrix<double>(8, 8)) );
+  }
+}
diff --git a/unsupported/test/sparse_ldlt.cpp b/unsupported/test/sparse_ldlt.cpp
new file mode 100644
index 000000000..8671807c9
--- /dev/null
+++ b/unsupported/test/sparse_ldlt.cpp
@@ -0,0 +1,68 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2010 Gael Guennebaud <g.gael@free.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/>.
+
+#include "sparse.h"
+
+#ifdef EIGEN_TAUCS_SUPPORT
+#include <Eigen/TaucsSupport>
+#endif
+
+template<typename Scalar> void sparse_ldlt(int rows, int cols)
+{
+  double density = std::max(8./(rows*cols), 0.01);
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+  typedef Matrix<Scalar,Dynamic,1> DenseVector;
+
+  SparseMatrix<Scalar> m2(rows, cols);
+  DenseMatrix refMat2(rows, cols);
+
+  DenseVector b = DenseVector::Random(cols);
+  DenseVector refX(cols), x(cols);
+
+  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0);
+  for(int i=0; i<rows; ++i)
+    m2.coeffRef(i,i) = refMat2(i,i) = ei_abs(ei_real(refMat2(i,i)));
+
+  refX = refMat2.template selfadjointView<Upper>().ldlt().solve(b);
+  typedef SparseMatrix<Scalar,Upper|SelfAdjoint> SparseSelfAdjointMatrix;
+  x = b;
+  SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
+  if (ldlt.succeeded())
+    ldlt.solveInPlace(x);
+  else
+    std::cerr << "warning LDLT failed\n";
+
+  VERIFY_IS_APPROX(refMat2.template selfadjointView<Upper>() * x, b);
+  VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default");
+}
+
+void test_sparse_ldlt()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST_1(sparse_ldlt<double>(8, 8) );
+    int s = ei_random<int>(1,300);
+    CALL_SUBTEST_2(sparse_ldlt<std::complex<double> >(s,s) );
+    CALL_SUBTEST_1(sparse_ldlt<double>(s,s) );
+  }
+}
diff --git a/unsupported/test/sparse_llt.cpp b/unsupported/test/sparse_llt.cpp
new file mode 100644
index 000000000..97360f1ec
--- /dev/null
+++ b/unsupported/test/sparse_llt.cpp
@@ -0,0 +1,91 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2010 Gael Guennebaud <g.gael@free.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/>.
+
+#include "sparse.h"
+
+#ifdef EIGEN_CHOLMOD_SUPPORT
+#include <Eigen/CholmodSupport>
+#endif
+
+#ifdef EIGEN_TAUCS_SUPPORT
+#include <Eigen/TaucsSupport>
+#endif
+
+template<typename Scalar> void sparse_llt(int rows, int cols)
+{
+  double density = std::max(8./(rows*cols), 0.01);
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+  typedef Matrix<Scalar,Dynamic,1> DenseVector;
+
+    // TODO fix the issue with complex (see SparseLLT::solveInPlace)
+    SparseMatrix<Scalar> m2(rows, cols);
+    DenseMatrix refMat2(rows, cols);
+
+    DenseVector b = DenseVector::Random(cols);
+    DenseVector refX(cols), x(cols);
+
+    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, 0, 0);
+    for(int i=0; i<rows; ++i)
+      m2.coeffRef(i,i) = refMat2(i,i) = ei_abs(ei_real(refMat2(i,i)));
+
+    refX = refMat2.template selfadjointView<Lower>().llt().solve(b);
+    if (!NumTraits<Scalar>::IsComplex)
+    {
+      x = b;
+      SparseLLT<SparseMatrix<Scalar> > (m2).solveInPlace(x);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
+    }
+    #ifdef EIGEN_CHOLMOD_SUPPORT
+    x = b;
+    SparseLLT<SparseMatrix<Scalar> ,Cholmod>(m2).solveInPlace(x);
+    VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod");
+    #endif
+
+    #ifdef EIGEN_TAUCS_SUPPORT
+    // TODO fix TAUCS with complexes
+    if (!NumTraits<Scalar>::IsComplex)
+    {
+      x = b;
+//       SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,IncompleteFactorization).solveInPlace(x);
+//       VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
+
+      x = b;
+      SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
+      x = b;
+      SparseLLT<SparseMatrix<Scalar> ,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x);
+      VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)");
+    }
+    #endif
+}
+
+void test_sparse_llt()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST_1(sparse_llt<double>(8, 8) );
+    int s = ei_random<int>(1,300);
+    CALL_SUBTEST_2(sparse_llt<std::complex<double> >(s,s) );
+    CALL_SUBTEST_1(sparse_llt<double>(s,s) );
+  }
+}
diff --git a/unsupported/test/sparse_lu.cpp b/unsupported/test/sparse_lu.cpp
new file mode 100644
index 000000000..bdf79ca23
--- /dev/null
+++ b/unsupported/test/sparse_lu.cpp
@@ -0,0 +1,115 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2010 Gael Guennebaud <g.gael@free.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/>.
+
+#include "sparse.h"
+
+#ifdef EIGEN_UMFPACK_SUPPORT
+#include <Eigen/UmfPackSupport>
+#endif
+
+#ifdef EIGEN_SUPERLU_SUPPORT
+#include <Eigen/SuperLUSupport>
+#endif
+
+template<typename Scalar> void sparse_lu(int rows, int cols)
+{
+  double density = std::max(8./(rows*cols), 0.01);
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+  typedef Matrix<Scalar,Dynamic,1> DenseVector;
+
+  DenseVector vec1 = DenseVector::Random(rows);
+
+  std::vector<Vector2i> zeroCoords;
+  std::vector<Vector2i> nonzeroCoords;
+
+    static int count = 0;
+    SparseMatrix<Scalar> m2(rows, cols);
+    DenseMatrix refMat2(rows, cols);
+
+    DenseVector b = DenseVector::Random(cols);
+    DenseVector refX(cols), x(cols);
+
+    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords);
+
+    FullPivLU<DenseMatrix> refLu(refMat2);
+    refX = refLu.solve(b);
+    #if defined(EIGEN_SUPERLU_SUPPORT) || defined(EIGEN_UMFPACK_SUPPORT)
+    Scalar refDet = refLu.determinant();
+    #endif
+    x.setZero();
+    // // SparseLU<SparseMatrix<Scalar> > (m2).solve(b,&x);
+    // // VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: default");
+    #ifdef EIGEN_SUPERLU_SUPPORT
+    {
+      x.setZero();
+      SparseLU<SparseMatrix<Scalar>,SuperLU> slu(m2);
+      if (slu.succeeded())
+      {
+        if (slu.solve(b,&x)) {
+          VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: SuperLU");
+        }
+        // std::cerr << refDet << " == " << slu.determinant() << "\n";
+        if (slu.solve(b, &x, SvTranspose)) {
+          VERIFY(b.isApprox(m2.transpose() * x, test_precision<Scalar>()));
+        }
+
+        if (slu.solve(b, &x, SvAdjoint)) {
+         VERIFY(b.isApprox(m2.adjoint() * x, test_precision<Scalar>()));
+        }
+
+        if (count==0) {
+          VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex
+        }
+      }
+    }
+    #endif
+    #ifdef EIGEN_UMFPACK_SUPPORT
+    {
+      // check solve
+      x.setZero();
+      SparseLU<SparseMatrix<Scalar>,UmfPack> slu(m2);
+      if (slu.succeeded()) {
+        if (slu.solve(b,&x)) {
+          if (count==0) {
+            VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: umfpack");  // FIXME solve is not very stable for complex
+          }
+        }
+        VERIFY_IS_APPROX(refDet,slu.determinant());
+        // TODO check the extracted data
+        //std::cerr << slu.matrixL() << "\n";
+      }
+    }
+    #endif
+    count++;
+}
+
+void test_sparse_lu()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST_1(sparse_lu<double>(8, 8) );
+    int s = ei_random<int>(1,300);
+    CALL_SUBTEST_2(sparse_lu<std::complex<double> >(s,s) );
+    CALL_SUBTEST_1(sparse_lu<double>(s,s) );
+  }
+}