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/sparse_product.cpp | |
parent | e0be0206222d3891c1511fc2ec4f4d41c3ccdf84 (diff) |
- add diagonal * sparse product as an expression
- split sparse_basic unit test
- various fixes in sparse module
Diffstat (limited to 'test/sparse_product.cpp')
-rw-r--r-- | test/sparse_product.cpp | 132 |
1 files changed, 132 insertions, 0 deletions
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)) ); + } +} |