sparse matrix must be row major like Scilab 5. 30/15230/13
Cedric Delamarre [Wed, 17 Sep 2014 07:56:59 +0000 (09:56 +0200)]
Change-Id: I39ecfa0f0684020cc14fe7dd07b4adc2ea04042e

scilab/modules/ast/includes/types/matrixiterator.hxx
scilab/modules/ast/includes/types/sparse.hxx
scilab/modules/ast/includes/types/sparseOp.hxx
scilab/modules/ast/includes/types/type_traits.hxx
scilab/modules/ast/src/cpp/types/sparse.cpp
scilab/modules/sparse/sci_gateway/cpp/sci_sp2adj.cpp

index 4485e9b..3f4de4b 100644 (file)
@@ -84,25 +84,25 @@ template<> std::complex<double> get(types::Sparse SPARSE_CONST& s, int r, int c)
     return s.get(r, c);
 }
 
-template<> double get(Eigen::SparseMatrix<double, 0, int> SPARSE_CONST&s, int r, int c)
+template<> double get(types::Sparse::RealSparse_t SPARSE_CONST&s, int r, int c)
 {
     return s.coeff(r, c);
 }
-template<> std::complex<double> get(Eigen::SparseMatrix<double, 0, int> SPARSE_CONST&s, int r, int c)
+template<> std::complex<double> get(types::Sparse::RealSparse_t SPARSE_CONST&s, int r, int c)
 {
     return std::complex<double>(s.coeff(r, c), 0.);
 }
 
-template<> bool get(Eigen::SparseMatrix<bool> SPARSE_CONST& d, int r, int c)
+template<> bool get(types::SparseBool::BoolSparse_t SPARSE_CONST& d, int r, int c)
 {
     return d.coeff(r, c);
 }
 
-template<> double get(Eigen::SparseMatrix<std::complex<double>, 0, int> SPARSE_CONST&s, int r, int c)
+template<> double get(types::Sparse::CplxSparse_t SPARSE_CONST&s, int r, int c)
 {
     return s.coeff(r, c).real();
 }
-template<> std::complex<double> get(Eigen::SparseMatrix<std::complex<double>, 0, int> SPARSE_CONST&s, int r, int c)
+template<> std::complex<double> get(types::Sparse::CplxSparse_t SPARSE_CONST&s, int r, int c)
 {
     return s.coeff(r, c);
 }
@@ -160,7 +160,7 @@ template<> bool set(types::SparseBool & d, int r, int c, int v)
 /*
  * TODO report possible bug in Eigen when inserting 0. invalidates Eigen::InnerIterator
  */
-template<> bool set(Eigen::SparseMatrix<double, 0, int>& s, int r, int c, double v)
+template<> bool set(types::Sparse::RealSparse_t& s, int r, int c, double v)
 {
     if (v != 0.)
     {
@@ -169,7 +169,7 @@ template<> bool set(Eigen::SparseMatrix<double, 0, int>& s, int r, int c, double
     return true;
 }
 
-template<> bool set(Eigen::SparseMatrix<double, 0, int>& s, int r, int c, std::complex<double> v)
+template<> bool set(types::Sparse::RealSparse_t& s, int r, int c, std::complex<double> v)
 {
     if ( v.real() != 0.)
     {
@@ -178,7 +178,7 @@ template<> bool set(Eigen::SparseMatrix<double, 0, int>& s, int r, int c, std::c
     return  true;
 }
 // should we make this a compile error ?
-template<> bool set(Eigen::SparseMatrix<std::complex<double>, 0, int>& s, int r, int c, double v)
+template<> bool set(types::Sparse::CplxSparse_t& s, int r, int c, double v)
 {
     if (v != 0.)
     {
@@ -191,7 +191,7 @@ namespace
 {
 std::complex<double> const cplxZero(0., 0.);
 }
-template<> bool set(Eigen::SparseMatrix<std::complex<double>, 0, int>& s, int r, int c, std::complex<double> v)
+template<> bool set(types::Sparse::CplxSparse_t& s, int r, int c, std::complex<double> v)
 {
     if (v != cplxZero)
     {
@@ -200,7 +200,7 @@ template<> bool set(Eigen::SparseMatrix<std::complex<double>, 0, int>& s, int r,
     return true;
 }
 
-template<> bool set(Eigen::SparseMatrix<bool>& s, int r, int c, bool v)
+template<> bool set(types::SparseBool::BoolSparse_t& s, int r, int c, bool v)
 {
     if (v)
     {
index d350494..1f5cfef 100644 (file)
@@ -495,8 +495,8 @@ struct EXTERN_AST Sparse : GenericType
 
     SparseBool* newLesserThan(Sparse const&o);
 
-    typedef Eigen::SparseMatrix<double >   RealSparse_t;
-    typedef Eigen::SparseMatrix<std::complex<double > >    CplxSparse_t;
+    typedef Eigen::SparseMatrix<double, Eigen::RowMajor>                 RealSparse_t;
+    typedef Eigen::SparseMatrix<std::complex<double>, Eigen::RowMajor>   CplxSparse_t;
     /**
        One and only one of the args should be 0.
        @param realSp ptr to an Eigen sparse matrix of double values
@@ -709,7 +709,7 @@ struct EXTERN_AST SparseBool : GenericType
     SparseBool* newLogicalOr(SparseBool const&o) const;
     SparseBool* newLogicalAnd(SparseBool const&o) const;
 
-    typedef Eigen::SparseMatrix<bool> BoolSparse_t;
+    typedef Eigen::SparseMatrix<bool, Eigen::RowMajor> BoolSparse_t;
     SparseBool(BoolSparse_t* o);
     BoolSparse_t* matrixBool;
 
index 81eb6e7..53373d9 100644 (file)
@@ -158,11 +158,11 @@ template<typename Scalar> struct OperatorTraits<std::logical_or<Scalar> > : With
  * @return ptr to the new Eigen::Sparse result containing the cwise binary op with the scalar.
  */
 template<typename Sp1, typename Scalar2, typename Op>
-Eigen::SparseMatrix<typename Op::result_type>* scalarOp(Eigen::EigenBase<Sp1> const& op1, Scalar2 op2, Op op)
+Eigen::SparseMatrix<typename Op::result_type, Eigen::RowMajor>* scalarOp(Eigen::EigenBase<Sp1> const& op1, Scalar2 op2, Op op)
 {
     typedef typename Eigen::internal::traits<Sp1>::Scalar Scalar1;
     typedef typename Op::result_type result_scalar;
-    typedef Eigen::SparseMatrix<result_scalar> result_t;
+    typedef Eigen::SparseMatrix<result_scalar, Eigen::RowMajor> result_t;
     result_t* res;
     if (op(Scalar1(), op2) == result_scalar())
     {
@@ -188,11 +188,11 @@ Eigen::SparseMatrix<typename Op::result_type>* scalarOp(Eigen::EigenBase<Sp1> co
  * @return ptr to the new Eigen::Sparse result containing the cwise binary op with the scalar.
  */
 template<typename Scalar1, typename Sp2, typename Op>
-Eigen::SparseMatrix<typename Op::result_type>* scalarOp(Scalar1 op1, Eigen::EigenBase<Sp2> const& op2, Op op)
+Eigen::SparseMatrix<typename Op::result_type, Eigen::RowMajor>* scalarOp(Scalar1 op1, Eigen::EigenBase<Sp2> const& op2, Op op)
 {
     typedef typename Eigen::internal::traits<Sp2>::Scalar Scalar2;
     typedef typename Op::result_type result_scalar;
-    typedef Eigen::SparseMatrix<result_scalar> result_t;
+    typedef Eigen::SparseMatrix<result_scalar, Eigen::RowMajor> result_t;
     result_t* res;
     if (op(op1, Scalar2()) == result_scalar())
     {
@@ -217,7 +217,7 @@ Eigen::SparseMatrix<typename Op::result_type>* scalarOp(Scalar1 op1, Eigen::Eige
  * @return ptr to the new Eigen::Sparse result containing the cwise binary op with the two matrices.
  */
 template< typename Sp1, typename Sp2 , typename Op>
-Eigen::SparseMatrix<typename Op::result_type>* cwiseOp(Eigen::EigenBase<Sp1> const& op1, Eigen::EigenBase<Sp2> const& op2, Op op)
+Eigen::SparseMatrix<typename Op::result_type, Eigen::RowMajor>* cwiseOp(Eigen::EigenBase<Sp1> const& op1, Eigen::EigenBase<Sp2> const& op2, Op op)
 {
     if (op1.rows() == 1 && op1.cols() == 1)
     {
@@ -239,10 +239,10 @@ Eigen::SparseMatrix<typename Op::result_type>* cwiseOp(Eigen::EigenBase<Sp1> con
  * @return ptr to the new Eigen::Sparse result containing the cwise binary op between the two matrices
  */
 template< typename Sp1, typename Sp2 , typename Op>
-Eigen::SparseMatrix<typename Op::result_type>* cwiseOp(Eigen::EigenBase<Sp1> const& op1, Eigen::EigenBase<Sp2> const& op2, Op op, FullTraversal /*unused*/)
+Eigen::SparseMatrix<typename Op::result_type, Eigen::RowMajor>* cwiseOp(Eigen::EigenBase<Sp1> const& op1, Eigen::EigenBase<Sp2> const& op2, Op op, FullTraversal /*unused*/)
 {
     typedef typename Op::result_type result_scalar;
-    typedef Eigen::SparseMatrix<result_scalar> result_t;
+    typedef Eigen::SparseMatrix<result_scalar, Eigen::RowMajor> result_t;
     typedef typename Eigen::internal::traits<Sp1>::Scalar Scalar1;
     typedef typename Eigen::internal::traits<Sp2>::Scalar Scalar2;
     // TODO remove dense temp when Eigen provides sparse full traversal API
@@ -261,10 +261,10 @@ Eigen::SparseMatrix<typename Op::result_type>* cwiseOp(Eigen::EigenBase<Sp1> con
  * @return ptr to the new Eigen::Sparse result containing the cwise binary op with the scalar.
  */
 template< typename Sp1, typename Sp2 , typename Op>
-Eigen::SparseMatrix<typename Op::result_type>* cwiseOp(Eigen::EigenBase<Sp1> const& op1, Eigen::EigenBase<Sp2> const& op2, Op op, UnionTraversal /*unused*/)
+Eigen::SparseMatrix<typename Op::result_type, Eigen::RowMajor>* cwiseOp(Eigen::EigenBase<Sp1> const& op1, Eigen::EigenBase<Sp2> const& op2, Op op, UnionTraversal /*unused*/)
 {
     typedef typename Op::result_type result_scalar;
-    typedef Eigen::SparseMatrix<result_scalar> result_t;
+    typedef Eigen::SparseMatrix<result_scalar, Eigen::RowMajor> result_t;
     result_t* res(new result_t(op1.derived().binaryExpr(op2.derived(), op)));
     res->prune(&keepForSparse<result_scalar>);
     return res;
index faedfde..f60176f 100644 (file)
@@ -77,7 +77,7 @@ struct type_traits
     }
 
     template<typename T>
-    inline static void neg(const int r, const int c, const T * const in, Eigen::SparseMatrix<bool> * const out)
+    inline static void neg(const int r, const int c, const T * const in, Eigen::SparseMatrix<bool, Eigen::RowMajor> * const out)
     {
         for (int i = 0; i < r; i++)
         {
index 16ce5d4..7e5b929 100644 (file)
@@ -134,35 +134,22 @@ template<typename T> std::wstring toString(T const& m, int precision)
     }
     ostr << L" sparse matrix\n\n";
 
-    const typename Eigen::internal::traits<T>::Index* pInner = m.innerIndexPtr();
-    const typename Eigen::internal::traits<T>::Index* pOuter = m.outerIndexPtr();
+    const typename Eigen::internal::traits<T>::Index* pIColPos      = m.innerIndexPtr();
+    const typename Eigen::internal::traits<T>::Index* pINbItemByRow = m.outerIndexPtr();
 
-    int iRow = 0;
-    int iCol = 0;
+    int iPos = 0;
 
-    for (size_t j = 0 ; j < m.rows() ; j++)
+    for (size_t j = 1 ; j < m.rows() + 1 ; j++)
     {
-        iRow = (int)j;
-        for (size_t i = 0 ; i < m.nonZeros() ; i++)
+        for (size_t i = pINbItemByRow[j - 1] ; i < pINbItemByRow[j] ; i++)
         {
-            if (pInner[i] == j)
-            {
-                //good row
-                for (size_t k = 0 ; k < m.outerSize() + 1; k++)
-                {
-                    if (pOuter[k] > i)
-                    {
-                        iCol = (int)k;
-                        break;
-                    }
-                }
+            ostr << L"(";
+            addUnsignedIntValue<unsigned long long>(&ostr, (int)j, iWidthRows);
+            ostr << L",";
+            addUnsignedIntValue<unsigned long long>(&ostr, pIColPos[iPos] + 1, iWidthCols);
+            ostr << L")\t" << p(m.valuePtr()[iPos]) << std::endl;
 
-                ostr << L"(";
-                addUnsignedIntValue<unsigned long long>(&ostr, iRow + 1, iWidthRows);
-                ostr << L",";
-                addUnsignedIntValue<unsigned long long>(&ostr, iCol, iWidthCols);
-                ostr << L")\t" << p(m.valuePtr()[i]) << std::endl;
-            }
+            iPos++;
         }
     }
 
@@ -1867,7 +1854,12 @@ int* Sparse::getColPos(int* _piColPos)
         mycopy_n(matrixReal->innerIndexPtr(), nonZeros(), _piColPos);
     }
 
-    std::transform(_piColPos, _piColPos + nonZeros(), _piColPos, std::bind2nd(std::plus<double>(), 1));
+    // std::transform(_piColPos, _piColPos + nonZeros(), _piColPos, std::bind2nd(std::plus<double>(), 1));
+    for (int i = 0; i < nonZeros(); i++)
+    {
+        _piColPos[i]++;
+    }
+
     return _piColPos;
 }
 
@@ -1881,10 +1873,10 @@ template<typename S> struct GetReal: std::unary_function<typename S::InnerIterat
         return it.value();
     }
 };
-template<> struct GetReal< Eigen::SparseMatrix<std::complex<double > > >
-        : std::unary_function<Eigen::SparseMatrix<std::complex<double > > ::InnerIterator, double>
+template<> struct GetReal< Eigen::SparseMatrix<std::complex<double >, Eigen::RowMajor > >
+        : std::unary_function<Sparse::CplxSparse_t::InnerIterator, double>
 {
-    double operator()( Eigen::SparseMatrix<std::complex<double > > ::InnerIterator it) const
+    double operator()( Sparse::CplxSparse_t::InnerIterator it) const
     {
         return it.value().real();
     }
@@ -3048,7 +3040,12 @@ int* SparseBool::getNbItemByRow(int* _piNbItemByRows)
 
 int* SparseBool::getColPos(int* _piColPos)
 {
-    mycopy_n(matrixBool->innerIndexPtr(), getRows(), _piColPos);
+    mycopy_n(matrixBool->innerIndexPtr(), nbTrue(), _piColPos);
+    for (int i = 0; i < nbTrue(); i++)
+    {
+        _piColPos[i]++;
+    }
+
     return _piColPos;
 }
 
index 6967401..edc596e 100644 (file)
@@ -45,7 +45,10 @@ Function::ReturnValue sci_sp2adj(typed_list &in, int nbRes, typed_list &out)
         return Function::Error;
     }
 
-    Sparse* SPARSE_CONST sp = in[0]->getAs<Sparse>();
+    InternalType* pIT = NULL;
+    Sparse* SPARSE_CONST spIn = in[0]->getAs<Sparse>();
+    spIn->transpose(pIT);
+    Sparse* sp = pIT->getAs<Sparse>();
     std::size_t const nonZeros = sp->nonZeros();
 
     types::Double* res = new Double(sp->getCols() + 1, 1);
@@ -75,5 +78,6 @@ Function::ReturnValue sci_sp2adj(typed_list &in, int nbRes, typed_list &out)
         out.push_back(res);
     }
 
+    delete pIT;
     return Function::OK;
 }