• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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@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 
sparse_product(const SparseMatrixType & ref)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 
test_eigen2_sparse_product()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