1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. Eigen itself is part of the KDE project. 3 // 4 // Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro (at) gmail.com> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #include "sparse.h" 11 12 template<typename SparseMatrixType> void sparse_product(const SparseMatrixType& ref) 13 { 14 const int rows = ref.rows(); 15 const int cols = ref.cols(); 16 typedef typename SparseMatrixType::Scalar Scalar; 17 enum { Flags = SparseMatrixType::Flags }; 18 19 double density = std::max(8./(rows*cols), 0.01); 20 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 21 typedef Matrix<Scalar,Dynamic,1> DenseVector; 22 23 // test matrix-matrix product 24 { 25 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); 26 DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows); 27 DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows); 28 DenseMatrix dm4 = DenseMatrix::Zero(rows, rows); 29 SparseMatrixType m2(rows, rows); 30 SparseMatrixType m3(rows, rows); 31 SparseMatrixType m4(rows, rows); 32 initSparse<Scalar>(density, refMat2, m2); 33 initSparse<Scalar>(density, refMat3, m3); 34 initSparse<Scalar>(density, refMat4, m4); 35 VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3); 36 VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3); 37 VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); 38 VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); 39 40 // sparse * dense 41 VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3); 42 VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose()); 43 VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3, refMat4=refMat2.transpose()*refMat3); 44 VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); 45 46 // dense * sparse 47 VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3); 48 VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); 49 VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3); 50 VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); 51 52 VERIFY_IS_APPROX(m3=m3*m3, refMat3=refMat3*refMat3); 53 } 54 55 // test matrix - diagonal product 56 if(false) // it compiles, but the precision is terrible. probably doesn't matter in this branch.... 57 { 58 DenseMatrix refM2 = DenseMatrix::Zero(rows, rows); 59 DenseMatrix refM3 = DenseMatrix::Zero(rows, rows); 60 DiagonalMatrix<DenseVector> d1(DenseVector::Random(rows)); 61 SparseMatrixType m2(rows, rows); 62 SparseMatrixType m3(rows, rows); 63 initSparse<Scalar>(density, refM2, m2); 64 initSparse<Scalar>(density, refM3, m3); 65 VERIFY_IS_APPROX(m3=m2*d1, refM3=refM2*d1); 66 VERIFY_IS_APPROX(m3=m2.transpose()*d1, refM3=refM2.transpose()*d1); 67 VERIFY_IS_APPROX(m3=d1*m2, refM3=d1*refM2); 68 VERIFY_IS_APPROX(m3=d1*m2.transpose(), refM3=d1 * refM2.transpose()); 69 } 70 71 // test self adjoint products 72 { 73 DenseMatrix b = DenseMatrix::Random(rows, rows); 74 DenseMatrix x = DenseMatrix::Random(rows, rows); 75 DenseMatrix refX = DenseMatrix::Random(rows, rows); 76 DenseMatrix refUp = DenseMatrix::Zero(rows, rows); 77 DenseMatrix refLo = DenseMatrix::Zero(rows, rows); 78 DenseMatrix refS = DenseMatrix::Zero(rows, rows); 79 SparseMatrixType mUp(rows, rows); 80 SparseMatrixType mLo(rows, rows); 81 SparseMatrixType mS(rows, rows); 82 do { 83 initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular); 84 } while (refUp.isZero()); 85 refLo = refUp.transpose().conjugate(); 86 mLo = mUp.transpose().conjugate(); 87 refS = refUp + refLo; 88 refS.diagonal() *= 0.5; 89 mS = mUp + mLo; 90 for (int k=0; k<mS.outerSize(); ++k) 91 for (typename SparseMatrixType::InnerIterator it(mS,k); it; ++it) 92 if (it.index() == k) 93 it.valueRef() *= 0.5; 94 95 VERIFY_IS_APPROX(refS.adjoint(), refS); 96 VERIFY_IS_APPROX(mS.transpose().conjugate(), mS); 97 VERIFY_IS_APPROX(mS, refS); 98 VERIFY_IS_APPROX(x=mS*b, refX=refS*b); 99 VERIFY_IS_APPROX(x=mUp.template marked<UpperTriangular|SelfAdjoint>()*b, refX=refS*b); 100 VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b); 101 VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b); 102 } 103 104 } 105 106 void test_eigen2_sparse_product() 107 { 108 for(int i = 0; i < g_repeat; i++) { 109 CALL_SUBTEST_1( sparse_product(SparseMatrix<double>(8, 8)) ); 110 CALL_SUBTEST_2( sparse_product(SparseMatrix<std::complex<double> >(16, 16)) ); 111 CALL_SUBTEST_1( sparse_product(SparseMatrix<double>(33, 33)) ); 112 113 CALL_SUBTEST_3( sparse_product(DynamicSparseMatrix<double>(8, 8)) ); 114 } 115 } 116
