extend sparse product unit tests

4 files changed, 109 insertions, 61 deletions

diff --git a/test/sparse.h b/test/sparse.h
index 949a597fc..530ae30bc 100644
--- a/test/sparse.h
+++ b/test/sparse.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
+// 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
@@ -58,30 +58,35 @@ enum {
  * \param zeroCoords and nonzeroCoords allows to get the coordinate lists of the non zero,
  *        and zero coefficients respectively.
  */
-template<typename Scalar> void
+template<typename Scalar,int Opt1,int Opt2> void
 initSparse(double density,
-           Matrix<Scalar,Dynamic,Dynamic>& refMat,
-           SparseMatrix<Scalar>& sparseMat,
+           Matrix<Scalar,Dynamic,Dynamic,Opt1>& refMat,
+           SparseMatrix<Scalar,Opt2>& sparseMat,
            int flags = 0,
            std::vector<Vector2i>* zeroCoords = 0,
            std::vector<Vector2i>* nonzeroCoords = 0)
 {
+  enum { IsRowMajor = SparseMatrix<Scalar,Opt2>::IsRowMajor };
   sparseMat.setZero();
   sparseMat.reserve(int(refMat.rows()*refMat.cols()*density));
-  for(int j=0; j<refMat.cols(); j++)
+  
+  for(int j=0; j<sparseMat.outerSize(); j++)
   {
     sparseMat.startVec(j);
-    for(int i=0; i<refMat.rows(); i++)
+    for(int i=0; i<sparseMat.innerSize(); i++)
     {
+      int ai(i), aj(j);
+      if(IsRowMajor)
+        std::swap(ai,aj);
       Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
       if ((flags&ForceNonZeroDiag) && (i==j))
       {
         v = internal::random<Scalar>()*Scalar(3.);
         v = v*v + Scalar(5.);
       }
-      if ((flags & MakeLowerTriangular) && j>i)
+      if ((flags & MakeLowerTriangular) && aj>ai)
         v = Scalar(0);
-      else if ((flags & MakeUpperTriangular) && j<i)
+      else if ((flags & MakeUpperTriangular) && aj<ai)
         v = Scalar(0);
 
       if ((flags&ForceRealDiag) && (i==j))
@@ -91,42 +96,46 @@ initSparse(double density,
       {
         sparseMat.insertBackByOuterInner(j,i) = v;
         if (nonzeroCoords)
-          nonzeroCoords->push_back(Vector2i(i,j));
+          nonzeroCoords->push_back(Vector2i(ai,aj));
       }
       else if (zeroCoords)
       {
-        zeroCoords->push_back(Vector2i(i,j));
+        zeroCoords->push_back(Vector2i(ai,aj));
       }
-      refMat(i,j) = v;
+      refMat(ai,aj) = v;
     }
   }
   sparseMat.finalize();
 }
 
-template<typename Scalar> void
+template<typename Scalar,int Opt1,int Opt2> void
 initSparse(double density,
-           Matrix<Scalar,Dynamic,Dynamic>& refMat,
-           DynamicSparseMatrix<Scalar>& sparseMat,
+           Matrix<Scalar,Dynamic,Dynamic, Opt1>& refMat,
+           DynamicSparseMatrix<Scalar, Opt2>& sparseMat,
            int flags = 0,
            std::vector<Vector2i>* zeroCoords = 0,
            std::vector<Vector2i>* nonzeroCoords = 0)
 {
+  enum { IsRowMajor = DynamicSparseMatrix<Scalar,Opt2>::IsRowMajor };
   sparseMat.setZero();
   sparseMat.reserve(int(refMat.rows()*refMat.cols()*density));
-  for(int j=0; j<refMat.cols(); j++)
+  for(int j=0; j<sparseMat.outerSize(); j++)
   {
     sparseMat.startVec(j); // not needed for DynamicSparseMatrix
-    for(int i=0; i<refMat.rows(); i++)
+    for(int i=0; i<sparseMat.innerSize(); i++)
     {
+      int ai(i), aj(j);
+      if(IsRowMajor)
+        std::swap(ai,aj);
       Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
       if ((flags&ForceNonZeroDiag) && (i==j))
       {
         v = internal::random<Scalar>()*Scalar(3.);
         v = v*v + Scalar(5.);
       }
-      if ((flags & MakeLowerTriangular) && j>i)
+      if ((flags & MakeLowerTriangular) && aj>ai)
         v = Scalar(0);
-      else if ((flags & MakeUpperTriangular) && j<i)
+      else if ((flags & MakeUpperTriangular) && aj<ai)
         v = Scalar(0);
 
       if ((flags&ForceRealDiag) && (i==j))
@@ -136,13 +145,13 @@ initSparse(double density,
       {
         sparseMat.insertBackByOuterInner(j,i) = v;
         if (nonzeroCoords)
-          nonzeroCoords->push_back(Vector2i(i,j));
+          nonzeroCoords->push_back(Vector2i(ai,aj));
       }
       else if (zeroCoords)
       {
-        zeroCoords->push_back(Vector2i(i,j));
+        zeroCoords->push_back(Vector2i(ai,aj));
       }
-      refMat(i,j) = v;
+      refMat(ai,aj) = v;
     }
   }
   sparseMat.finalize();
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index 9d79ca740..7910bbf8f 100644
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -1,6 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
+// Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
 // Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp
index 90ec3781e..1d2183bc3 100644
--- a/test/sparse_product.cpp
+++ b/test/sparse_product.cpp
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
+// 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
@@ -24,11 +24,37 @@
 
 #include "sparse.h"
 
-template<typename SparseMatrixType> void sparse_product(const SparseMatrixType& ref)
+template<typename SparseMatrixType, typename DenseMatrix, bool IsRowMajor=SparseMatrixType::IsRowMajor> struct test_outer;
+
+template<typename SparseMatrixType, typename DenseMatrix> struct test_outer<SparseMatrixType,DenseMatrix,false> {
+  static void run(SparseMatrixType& m2, SparseMatrixType& m4, DenseMatrix& refMat2, DenseMatrix& refMat4) {
+    int c  = internal::random(0,m2.cols()-1);
+    int c1 = internal::random(0,m2.cols()-1);
+    VERIFY_IS_APPROX(m4=m2.col(c)*refMat2.col(c1).transpose(), refMat4=refMat2.col(c)*refMat2.col(c1).transpose());
+    VERIFY_IS_APPROX(m4=refMat2.col(c1)*m2.col(c).transpose(), refMat4=refMat2.col(c1)*refMat2.col(c).transpose());
+  }
+};
+
+template<typename SparseMatrixType, typename DenseMatrix> struct test_outer<SparseMatrixType,DenseMatrix,true> {
+  static void run(SparseMatrixType& m2, SparseMatrixType& m4, DenseMatrix& refMat2, DenseMatrix& refMat4) {
+    int r  = internal::random(0,m2.rows()-1);
+    int c1 = internal::random(0,m2.cols()-1);
+    VERIFY_IS_APPROX(m4=m2.row(r).transpose()*refMat2.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*refMat2.col(c1).transpose());
+    VERIFY_IS_APPROX(m4=refMat2.col(c1)*m2.row(r), refMat4=refMat2.col(c1)*refMat2.row(r));
+  }
+};
+
+// (m2,m4,refMat2,refMat4,dv1);
+//     VERIFY_IS_APPROX(m4=m2.innerVector(c)*dv1.transpose(), refMat4=refMat2.colVector(c)*dv1.transpose());
+//     VERIFY_IS_APPROX(m4=dv1*mcm.col(c).transpose(), refMat4=dv1*refMat2.col(c).transpose());
+
+template<typename SparseMatrixType> void sparse_product()
 {
   typedef typename SparseMatrixType::Index Index;
-  const Index rows = ref.rows();
-  const Index cols = ref.cols();
+  Index n = 100;
+  const Index rows  = internal::random<int>(1,n);
+  const Index cols  = internal::random<int>(1,n);
+  const Index depth = internal::random<int>(1,n);
   typedef typename SparseMatrixType::Scalar Scalar;
   enum { Flags = SparseMatrixType::Flags };
 
@@ -41,25 +67,37 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
 
   // test matrix-matrix product
   {
-    DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
-    DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows);
-    DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows);
-    DenseMatrix refMat5 = DenseMatrix::Random(rows, rows);
+    DenseMatrix refMat2  = DenseMatrix::Zero(rows, depth);
+    DenseMatrix refMat2t = DenseMatrix::Zero(depth, rows);
+    DenseMatrix refMat3  = DenseMatrix::Zero(depth, cols);
+    DenseMatrix refMat3t = DenseMatrix::Zero(cols, depth);
+    DenseMatrix refMat4  = DenseMatrix::Zero(rows, cols);
+    DenseMatrix refMat4t = DenseMatrix::Zero(cols, rows);
+    DenseMatrix refMat5  = DenseMatrix::Random(depth, cols);
+    DenseMatrix refMat6  = DenseMatrix::Random(rows, rows);
     DenseMatrix dm4 = DenseMatrix::Zero(rows, rows);
-    DenseVector dv1 = DenseVector::Random(rows);
-    SparseMatrixType m2(rows, rows);
-    SparseMatrixType m3(rows, rows);
-    SparseMatrixType m4(rows, rows);
-    initSparse<Scalar>(density, refMat2, m2);
-    initSparse<Scalar>(density, refMat3, m3);
-    initSparse<Scalar>(density, refMat4, m4);
-
-    int c = internal::random<int>(0,rows-1);
+//     DenseVector dv1 = DenseVector::Random(rows);
+    SparseMatrixType m2 (rows, depth);
+    SparseMatrixType m2t(depth, rows);
+    SparseMatrixType m3 (depth, cols);
+    SparseMatrixType m3t(cols, depth);
+    SparseMatrixType m4 (rows, cols);
+    SparseMatrixType m4t(cols, rows);
+    SparseMatrixType m6(rows, rows);
+    initSparse(density, refMat2,  m2);
+    initSparse(density, refMat2t, m2t);
+    initSparse(density, refMat3,  m3);
+    initSparse(density, refMat3t, m3t);
+    initSparse(density, refMat4,  m4);
+    initSparse(density, refMat4t, m4t);
+    initSparse(density, refMat6, m6);
+
+//     int c = internal::random<int>(0,depth-1);
 
     VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3);
-    VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
-    VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
-    VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose());
+    VERIFY_IS_APPROX(m4=m2t.transpose()*m3, refMat4=refMat2t.transpose()*refMat3);
+    VERIFY_IS_APPROX(m4=m2t.transpose()*m3t.transpose(), refMat4=refMat2t.transpose()*refMat3t.transpose());
+    VERIFY_IS_APPROX(m4=m2*m3t.transpose(), refMat4=refMat2*refMat3t.transpose());
 
     VERIFY_IS_APPROX(m4 = m2*m3/s1, refMat4 = refMat2*refMat3/s1);
     VERIFY_IS_APPROX(m4 = m2*m3*s1, refMat4 = refMat2*refMat3*s1);
@@ -67,24 +105,23 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
 
     // sparse * dense
     VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3);
-    VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose());
-    VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3, refMat4=refMat2.transpose()*refMat3);
-    VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
+    VERIFY_IS_APPROX(dm4=m2*refMat3t.transpose(), refMat4=refMat2*refMat3t.transpose());
+    VERIFY_IS_APPROX(dm4=m2t.transpose()*refMat3, refMat4=refMat2t.transpose()*refMat3);
+    VERIFY_IS_APPROX(dm4=m2t.transpose()*refMat3t.transpose(), refMat4=refMat2t.transpose()*refMat3t.transpose());
 
     VERIFY_IS_APPROX(dm4=m2*(refMat3+refMat3), refMat4=refMat2*(refMat3+refMat3));
-    VERIFY_IS_APPROX(dm4=m2.transpose()*(refMat3+refMat5)*0.5, refMat4=refMat2.transpose()*(refMat3+refMat5)*0.5);
+    VERIFY_IS_APPROX(dm4=m2t.transpose()*(refMat3+refMat5)*0.5, refMat4=refMat2t.transpose()*(refMat3+refMat5)*0.5);
 
     // dense * sparse
     VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3);
-    VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose());
-    VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
-    VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
+    VERIFY_IS_APPROX(dm4=refMat2*m3t.transpose(), refMat4=refMat2*refMat3t.transpose());
+    VERIFY_IS_APPROX(dm4=refMat2t.transpose()*m3, refMat4=refMat2t.transpose()*refMat3);
+    VERIFY_IS_APPROX(dm4=refMat2t.transpose()*m3t.transpose(), refMat4=refMat2t.transpose()*refMat3t.transpose());
 
     // sparse * dense and dense * sparse outer product
-    VERIFY_IS_APPROX(m4=m2.col(c)*dv1.transpose(), refMat4=refMat2.col(c)*dv1.transpose());
-    VERIFY_IS_APPROX(m4=dv1*m2.col(c).transpose(), refMat4=dv1*refMat2.col(c).transpose());
+    test_outer<SparseMatrixType,DenseMatrix>::run(m2,m4,refMat2,refMat4);
 
-    VERIFY_IS_APPROX(m3=m3*m3, refMat3=refMat3*refMat3);
+    VERIFY_IS_APPROX(m6=m6*m6, refMat6=refMat6*refMat6);
   }
 
   // test matrix - diagonal product
@@ -116,18 +153,19 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
     do {
       initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular);
     } while (refUp.isZero());
-    refLo = refUp.transpose().conjugate();
-    mLo = mUp.transpose().conjugate();
+    refLo = refUp.adjoint();
+    mLo = mUp.adjoint();
     refS = refUp + refLo;
     refS.diagonal() *= 0.5;
     mS = mUp + mLo;
+    // TODO be able to address the diagonal....
     for (int k=0; k<mS.outerSize(); ++k)
       for (typename SparseMatrixType::InnerIterator it(mS,k); it; ++it)
         if (it.index() == k)
           it.valueRef() *= 0.5;
 
     VERIFY_IS_APPROX(refS.adjoint(), refS);
-    VERIFY_IS_APPROX(mS.transpose().conjugate(), mS);
+    VERIFY_IS_APPROX(mS.adjoint(), mS);
     VERIFY_IS_APPROX(mS, refS);
     VERIFY_IS_APPROX(x=mS*b, refX=refS*b);
 
@@ -162,12 +200,12 @@ template<typename SparseMatrixType, typename DenseMatrixType> void sparse_produc
 void test_sparse_product()
 {
   for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST_1( sparse_product(SparseMatrix<double>(8, 8)) );
-    CALL_SUBTEST_2( sparse_product(SparseMatrix<std::complex<double> >(16, 16)) );
-    CALL_SUBTEST_1( sparse_product(SparseMatrix<double>(33, 33)) );
-
-    CALL_SUBTEST_3( sparse_product(DynamicSparseMatrix<double>(8, 8)) );
-
+    CALL_SUBTEST_1( (sparse_product<SparseMatrix<double,ColMajor> >()) );
+    CALL_SUBTEST_1( (sparse_product<SparseMatrix<double,RowMajor> >()) );
+    CALL_SUBTEST_2( (sparse_product<SparseMatrix<std::complex<double>, ColMajor > >()) );
+    CALL_SUBTEST_2( (sparse_product<SparseMatrix<std::complex<double>, RowMajor > >()) );
+    CALL_SUBTEST_3( (sparse_product<DynamicSparseMatrix<float, ColMajor> >()) );
+    CALL_SUBTEST_3( (sparse_product<DynamicSparseMatrix<float, RowMajor> >()) );
     CALL_SUBTEST_4( (sparse_product_regression_test<SparseMatrix<double,RowMajor>, Matrix<double, Dynamic, Dynamic, RowMajor> >()) );
   }
 }
diff --git a/test/sparse_vector.cpp b/test/sparse_vector.cpp
index be85740c0..b3249915e 100644
--- a/test/sparse_vector.cpp
+++ b/test/sparse_vector.cpp
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
+// 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