diff options
author | Gael Guennebaud <g.gael@free.fr> | 2009-02-09 09:59:30 +0000 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2009-02-09 09:59:30 +0000 |
commit | a9688f0b717568713ff1acef3466a0da00a68980 (patch) | |
tree | 8eea0db4c744fdc20c110cf767c884e80f2b83e4 /test | |
parent | e0be0206222d3891c1511fc2ec4f4d41c3ccdf84 (diff) |
- add diagonal * sparse product as an expression
- split sparse_basic unit test
- various fixes in sparse module
Diffstat (limited to 'test')
-rw-r--r-- | test/CMakeLists.txt | 1 | ||||
-rw-r--r-- | test/sparse_basic.cpp | 65 | ||||
-rw-r--r-- | test/sparse_product.cpp | 132 | ||||
-rw-r--r-- | test/sparse_vector.cpp | 1 |
4 files changed, 136 insertions, 63 deletions
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 99a21835f..4d2ac7501 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -146,6 +146,7 @@ if(QT4_FOUND) endif(QT4_FOUND) ei_add_test(sparse_vector) ei_add_test(sparse_basic) +ei_add_test(sparse_product) ei_add_test(sparse_solvers " " "${SPARSE_LIBS}") ei_add_test(reverse) diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 23b8526b7..410ef96a6 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -238,6 +238,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m1+=m2, refM1+=refM2); VERIFY_IS_APPROX(m1-=m2, refM1-=refM2); + VERIFY_IS_APPROX(m1.col(0).dot(refM2.row(0)), refM1.col(0).dot(refM2.row(0))); + refM4.setRandom(); // sparse cwise* dense VERIFY_IS_APPROX(m3.cwise()*refM4, refM3.cwise()*refM4); @@ -281,69 +283,6 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); } - - // test matrix product - { - DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); - DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows); - DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows); - DenseMatrix dm4 = DenseMatrix::Zero(rows, 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); - 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()); - - // 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()); - - // 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()); - } - - // test self adjoint products - { - DenseMatrix b = DenseMatrix::Random(rows, rows); - DenseMatrix x = DenseMatrix::Random(rows, rows); - DenseMatrix refX = DenseMatrix::Random(rows, rows); - DenseMatrix refUp = DenseMatrix::Zero(rows, rows); - DenseMatrix refLo = DenseMatrix::Zero(rows, rows); - DenseMatrix refS = DenseMatrix::Zero(rows, rows); - SparseMatrixType mUp(rows, rows); - SparseMatrixType mLo(rows, rows); - SparseMatrixType mS(rows, rows); - do { - initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular); - } while (refUp.isZero()); - refLo = refUp.transpose().conjugate(); - mLo = mUp.transpose().conjugate(); - refS = refUp + refLo; - refS.diagonal() *= 0.5; - mS = mUp + mLo; - 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, refS); - VERIFY_IS_APPROX(x=mS*b, refX=refS*b); - VERIFY_IS_APPROX(x=mUp.template marked<UpperTriangular|SelfAdjoint>()*b, refX=refS*b); - VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b); - VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b); - } // test prune { diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp new file mode 100644 index 000000000..95b3546c1 --- /dev/null +++ b/test/sparse_product.cpp @@ -0,0 +1,132 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com> +// +// 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" + +template<typename SparseMatrixType> void sparse_product(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; + + // test matrix-matrix product + /* + { + DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); + DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows); + DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows); + DenseMatrix dm4 = DenseMatrix::Zero(rows, 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); + 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()); + + // 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()); + + // 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()); + } + */ + + // test matrix - diagonal product + { + DenseMatrix refM2 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM3 = DenseMatrix::Zero(rows, rows); + DiagonalMatrix<Scalar,Dynamic> d1(DenseVector::Random(rows)); + SparseMatrixType m2(rows, rows); + SparseMatrixType m3(rows, rows); + initSparse<Scalar>(density, refM2, m2); + initSparse<Scalar>(density, refM3, m3); +// std::cerr << "foo\n" << (m2*d1).toDense() << "\n\n" << refM2*d1 << "\n\n"; + VERIFY_IS_APPROX(m3=m2*d1, refM3=refM2*d1); + VERIFY_IS_APPROX(m3=m2.transpose()*d1, refM3=refM2.transpose()*d1); + VERIFY_IS_APPROX(m3=d1*m2, refM3=d1*refM2); +// std::cerr << "foo\n" << (d1*m2.transpose()).toDense() << "\n\n" << d1 * refM2.transpose() << "\n\n"; + VERIFY_IS_APPROX(m3=d1*m2.transpose(), refM3=d1 * refM2.transpose()); + } + + // test self adjoint products +// { +// DenseMatrix b = DenseMatrix::Random(rows, rows); +// DenseMatrix x = DenseMatrix::Random(rows, rows); +// DenseMatrix refX = DenseMatrix::Random(rows, rows); +// DenseMatrix refUp = DenseMatrix::Zero(rows, rows); +// DenseMatrix refLo = DenseMatrix::Zero(rows, rows); +// DenseMatrix refS = DenseMatrix::Zero(rows, rows); +// SparseMatrixType mUp(rows, rows); +// SparseMatrixType mLo(rows, rows); +// SparseMatrixType mS(rows, rows); +// do { +// initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular); +// } while (refUp.isZero()); +// refLo = refUp.transpose().conjugate(); +// mLo = mUp.transpose().conjugate(); +// refS = refUp + refLo; +// refS.diagonal() *= 0.5; +// mS = mUp + mLo; +// 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, refS); +// VERIFY_IS_APPROX(x=mS*b, refX=refS*b); +// VERIFY_IS_APPROX(x=mUp.template marked<UpperTriangular|SelfAdjoint>()*b, refX=refS*b); +// VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b); +// VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b); +// } + +} + +void test_sparse_product() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( sparse_product(SparseMatrix<double>(8, 8)) ); + CALL_SUBTEST( sparse_product(SparseMatrix<std::complex<double> >(16, 16)) ); + CALL_SUBTEST( sparse_product(SparseMatrix<double>(33, 33)) ); + + CALL_SUBTEST( sparse_product(DynamicSparseMatrix<double>(8, 8)) ); + } +} diff --git a/test/sparse_vector.cpp b/test/sparse_vector.cpp index 8207e522a..934719f2c 100644 --- a/test/sparse_vector.cpp +++ b/test/sparse_vector.cpp @@ -84,6 +84,7 @@ template<typename Scalar> void sparse_vector(int rows, int cols) VERIFY_IS_APPROX(v1-=v2, refV1-=refV2); VERIFY_IS_APPROX(v1.dot(v2), refV1.dot(refV2)); + VERIFY_IS_APPROX(v1.dot(refV2), refV1.dot(refV2)); } |