diff options
author | Gael Guennebaud <g.gael@free.fr> | 2009-02-12 15:57:13 +0000 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2009-02-12 15:57:13 +0000 |
commit | 20a8bb96ebeebaa9353ffd09a8ff0ab0d6df4c3e (patch) | |
tree | d74c6cda6c89e464a458f5992ddab549d34af1e8 /test/sparse_product.cpp | |
parent | 59a1ed0932a06fd91d444e563479e42b3df432f3 (diff) |
fix m = m*m with m sparse (gug found by Frederik Heinz)
Diffstat (limited to 'test/sparse_product.cpp')
-rw-r--r-- | test/sparse_product.cpp | 76 |
1 files changed, 37 insertions, 39 deletions
diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index 95b3546c1..d2c3b1c98 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp @@ -30,14 +30,13 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType& 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); @@ -65,9 +64,10 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType& 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(m3=m3*m3, refMat3=refMat3*refMat3); } - */ - + // test matrix - diagonal product { DenseMatrix refM2 = DenseMatrix::Zero(rows, rows); @@ -77,46 +77,44 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType& 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); -// } + { + 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); + } } @@ -126,7 +124,7 @@ void test_sparse_product() 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)) ); } } |