mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-18 14:34:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

extend sparse product unit tests

Browse Source
This commit is contained in:
Gael Guennebaud 2011-03-22 11:58:22 +01:00
parent 5fda8cdfb3
commit 22c7609d72
4 changed files with 108 additions and 60 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
test/sparse.h
View File

@ -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();

1
test/sparse_basic.cpp
View File

@ -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

114
test/sparse_product.cpp
View File

@ -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);
// 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 @@ 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> >()) );
}
}

2
test/sparse_vector.cpp
View File

@ -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