extend sparse product unit tests
diff --git a/test/sparse.h b/test/sparse.h
index 949a597..530ae30 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 @@
* \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 @@
{
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 @@
{
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 9d79ca7..7910bbf 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 90ec378..1d2183b 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 @@
// 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);
+// 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,rows-1);
+// 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 @@
// 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 @@
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 @@
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 be85740..b324991 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
