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 |
|
---|
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 | }
|
---|