Macros | Functions | Variables
test/sparse_product.cpp File Reference
#include "sparse.h"

Macros

#define EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN    { on_temporary_creation(); }
 
#define VERIFY_EVALUATION_COUNT(XPR, N)
 

Functions

void on_temporary_creation ()
 
template<typename SparseMatrixType >
void sparse_product ()
 
template<typename SparseMatrixType , typename DenseMatrixType >
void sparse_product_regression_test ()
 
template<typename Scalar >
void bug_942 ()
 
template<typename Real >
void test_mixing_types ()
 
template<int OrderA, int OrderB, int OrderC>
void test_mixed_storage_imp ()
 
void test_mixed_storage ()
 
 EIGEN_DECLARE_TEST (sparse_product)
 

Variables

static long int nb_temporaries
 

Macro Definition Documentation

◆ EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN

#define EIGEN_SPARSE_CREATE_TEMPORARY_PLUGIN    { on_temporary_creation(); }

◆ VERIFY_EVALUATION_COUNT

#define VERIFY_EVALUATION_COUNT (   XPR,
  N 
)
Value:
{ \
nb_temporaries = 0; \
CALL_SUBTEST(XPR); \
if (nb_temporaries != N) std::cerr << "nb_temporaries == " << nb_temporaries << "\n"; \
VERIFY((#XPR) && nb_temporaries == N); \
}
N
@ N
Definition: constructor.cpp:22
nb_temporaries
static long int nb_temporaries
Definition: test/sparse_product.cpp:16

Function Documentation

◆ bug_942()

template<typename Scalar >
void bug_942 ( )
377  {
378  typedef Matrix<Scalar, Dynamic, 1> Vector;
379  typedef SparseMatrix<Scalar, ColMajor> ColSpMat;
380  typedef SparseMatrix<Scalar, RowMajor> RowSpMat;
381  ColSpMat cmA(1, 1);
382  cmA.insert(0, 0) = 1;
383 
384  RowSpMat rmA(1, 1);
385  rmA.insert(0, 0) = 1;
386 
387  Vector d(1);
388  d[0] = 2;
389 
390  double res = 2;
391 
392  VERIFY_IS_APPROX((cmA * d.asDiagonal()).eval().coeff(0, 0), res);
393  VERIFY_IS_APPROX((d.asDiagonal() * rmA).eval().coeff(0, 0), res);
394  VERIFY_IS_APPROX((rmA * d.asDiagonal()).eval().coeff(0, 0), res);
395  VERIFY_IS_APPROX((d.asDiagonal() * cmA).eval().coeff(0, 0), res);
396 }
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
Eigen::Vector
Matrix< Type, Size, 1 > Vector
\cpp11 Size×1 vector of type Type.
Definition: Eigen/Eigen/src/Core/Matrix.h:515
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:13

References res, and VERIFY_IS_APPROX.

◆ EIGEN_DECLARE_TEST()

EIGEN_DECLARE_TEST ( sparse_product  )
538  {
539  for (int i = 0; i < g_repeat; i++) {
540  CALL_SUBTEST_1((sparse_product<SparseMatrix<double, ColMajor> >()));
541  CALL_SUBTEST_1((sparse_product<SparseMatrix<double, RowMajor> >()));
542  CALL_SUBTEST_1((bug_942<double>()));
543  CALL_SUBTEST_2((sparse_product<SparseMatrix<std::complex<double>, ColMajor> >()));
544  CALL_SUBTEST_2((sparse_product<SparseMatrix<std::complex<double>, RowMajor> >()));
545  CALL_SUBTEST_3((sparse_product<SparseMatrix<float, ColMajor, long int> >()));
546  CALL_SUBTEST_4((
547  sparse_product_regression_test<SparseMatrix<double, RowMajor>, Matrix<double, Dynamic, Dynamic, RowMajor> >()));
548 
549  CALL_SUBTEST_5((test_mixing_types<float>()));
550  CALL_SUBTEST_5((test_mixed_storage()));
551  }
552 }
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::ColMajor
@ ColMajor
Definition: Constants.h:318
Eigen::RowMajor
@ RowMajor
Definition: Constants.h:320
Eigen::g_repeat
static int g_repeat
Definition: main.h:191
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
sparse_product_regression_test
void sparse_product_regression_test()
Definition: test/sparse_product.cpp:357
test_mixed_storage
void test_mixed_storage()
Definition: test/sparse_product.cpp:527
sparse_product
void sparse_product()
Definition: test/sparse_product.cpp:37

References CALL_SUBTEST_1, CALL_SUBTEST_2, CALL_SUBTEST_3, CALL_SUBTEST_4, CALL_SUBTEST_5, Eigen::ColMajor, Eigen::g_repeat, i, Eigen::RowMajor, sparse_product(), sparse_product_regression_test(), and test_mixed_storage().

◆ on_temporary_creation()

void on_temporary_creation ( )
inline
18  {
19  // here's a great place to set a breakpoint when debugging failures in this test!
20  nb_temporaries++;
21 }

References nb_temporaries.

◆ sparse_product()

template<typename SparseMatrixType >
void sparse_product ( )
37  {
38  typedef typename SparseMatrixType::StorageIndex StorageIndex;
39  Index n = 100;
40  const Index rows = internal::random<Index>(1, n);
41  const Index cols = internal::random<Index>(1, n);
42  const Index depth = internal::random<Index>(1, n);
43  typedef typename SparseMatrixType::Scalar Scalar;
44  enum { Flags = SparseMatrixType::Flags };
45 
46  double density = (std::max)(8. / (rows * cols), 0.2);
47  typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix;
48  typedef Matrix<Scalar, Dynamic, 1> DenseVector;
49  typedef Matrix<Scalar, 1, Dynamic> RowDenseVector;
50  typedef SparseVector<Scalar, 0, StorageIndex> ColSpVector;
51  typedef SparseVector<Scalar, RowMajor, StorageIndex> RowSpVector;
52 
53  Scalar s1 = internal::random<Scalar>();
54  Scalar s2 = internal::random<Scalar>();
55 
56  // test matrix-matrix product
57  {
58  DenseMatrix refMat2 = DenseMatrix::Zero(rows, depth);
59  DenseMatrix refMat2t = DenseMatrix::Zero(depth, rows);
60  DenseMatrix refMat3 = DenseMatrix::Zero(depth, cols);
61  DenseMatrix refMat3t = DenseMatrix::Zero(cols, depth);
62  DenseMatrix refMat4 = DenseMatrix::Zero(rows, cols);
63  DenseMatrix refMat4t = DenseMatrix::Zero(cols, rows);
64  DenseMatrix refMat5 = DenseMatrix::Random(depth, cols);
65  DenseMatrix refMat6 = DenseMatrix::Random(rows, rows);
66  DenseMatrix dm4 = DenseMatrix::Zero(rows, rows);
67  // DenseVector dv1 = DenseVector::Random(rows);
68  SparseMatrixType m2(rows, depth);
69  SparseMatrixType m2t(depth, rows);
70  SparseMatrixType m3(depth, cols);
71  SparseMatrixType m3t(cols, depth);
72  SparseMatrixType m4(rows, cols);
73  SparseMatrixType m4t(cols, rows);
74  SparseMatrixType m6(rows, rows);
75  initSparse(density, refMat2, m2);
76  initSparse(density, refMat2t, m2t);
77  initSparse(density, refMat3, m3);
78  initSparse(density, refMat3t, m3t);
79  initSparse(density, refMat4, m4);
80  initSparse(density, refMat4t, m4t);
81  initSparse(density, refMat6, m6);
82 
83  // int c = internal::random<int>(0,depth-1);
84 
85  // sparse * sparse
86  VERIFY_IS_APPROX(m4 = m2 * m3, refMat4 = refMat2 * refMat3);
87  VERIFY_IS_APPROX(m4 = m2t.transpose() * m3, refMat4 = refMat2t.transpose() * refMat3);
88  VERIFY_IS_APPROX(m4 = m2t.transpose() * m3t.transpose(), refMat4 = refMat2t.transpose() * refMat3t.transpose());
89  VERIFY_IS_APPROX(m4 = m2 * m3t.transpose(), refMat4 = refMat2 * refMat3t.transpose());
90 
91  VERIFY_IS_APPROX(m4 = m2 * m3 / s1, refMat4 = refMat2 * refMat3 / s1);
92  VERIFY_IS_APPROX(m4 = m2 * m3 * s1, refMat4 = refMat2 * refMat3 * s1);
93  VERIFY_IS_APPROX(m4 = s2 * m2 * m3 * s1, refMat4 = s2 * refMat2 * refMat3 * s1);
94  VERIFY_IS_APPROX(m4 = (m2 + m2) * m3, refMat4 = (refMat2 + refMat2) * refMat3);
95  VERIFY_IS_APPROX(m4 = m2 * m3.leftCols(cols / 2), refMat4 = refMat2 * refMat3.leftCols(cols / 2));
96  VERIFY_IS_APPROX(m4 = m2 * (m3 + m3).leftCols(cols / 2),
97  refMat4 = refMat2 * (refMat3 + refMat3).leftCols(cols / 2));
98 
99  VERIFY_IS_APPROX(m4 = (m2 * m3).pruned(0), refMat4 = refMat2 * refMat3);
100  VERIFY_IS_APPROX(m4 = (m2t.transpose() * m3).pruned(0), refMat4 = refMat2t.transpose() * refMat3);
101  VERIFY_IS_APPROX(m4 = (m2t.transpose() * m3t.transpose()).pruned(0),
102  refMat4 = refMat2t.transpose() * refMat3t.transpose());
103  VERIFY_IS_APPROX(m4 = (m2 * m3t.transpose()).pruned(0), refMat4 = refMat2 * refMat3t.transpose());
104 
105 #ifndef EIGEN_SPARSE_PRODUCT_IGNORE_TEMPORARY_COUNT
106  // make sure the right product implementation is called:
107  if ((!SparseMatrixType::IsRowMajor) && m2.rows() <= m3.cols()) {
108  VERIFY_EVALUATION_COUNT(m4 = m2 * m3, 2); // 2 for transposing and get a sorted result.
109  VERIFY_EVALUATION_COUNT(m4 = (m2 * m3).pruned(0), 1);
110  VERIFY_EVALUATION_COUNT(m4 = (m2 * m3).eval().pruned(0), 4);
111  }
112 #endif
113 
114  // and that pruning is effective:
115  {
116  DenseMatrix Ad(2, 2);
117  Ad << -1, 1, 1, 1;
118  SparseMatrixType As(Ad.sparseView()), B(2, 2);
119  VERIFY_IS_EQUAL((As * As.transpose()).eval().nonZeros(), 4);
120  VERIFY_IS_EQUAL((Ad * Ad.transpose()).eval().sparseView().eval().nonZeros(), 2);
121  VERIFY_IS_EQUAL((As * As.transpose()).pruned(1e-6).eval().nonZeros(), 2);
122  }
123 
124  // dense ?= sparse * sparse
125  VERIFY_IS_APPROX(dm4 = m2 * m3, refMat4 = refMat2 * refMat3);
126  VERIFY_IS_APPROX(dm4 += m2 * m3, refMat4 += refMat2 * refMat3);
127  VERIFY_IS_APPROX(dm4 -= m2 * m3, refMat4 -= refMat2 * refMat3);
128  VERIFY_IS_APPROX(dm4 = m2t.transpose() * m3, refMat4 = refMat2t.transpose() * refMat3);
129  VERIFY_IS_APPROX(dm4 += m2t.transpose() * m3, refMat4 += refMat2t.transpose() * refMat3);
130  VERIFY_IS_APPROX(dm4 -= m2t.transpose() * m3, refMat4 -= refMat2t.transpose() * refMat3);
131  VERIFY_IS_APPROX(dm4 = m2t.transpose() * m3t.transpose(), refMat4 = refMat2t.transpose() * refMat3t.transpose());
132  VERIFY_IS_APPROX(dm4 += m2t.transpose() * m3t.transpose(), refMat4 += refMat2t.transpose() * refMat3t.transpose());
133  VERIFY_IS_APPROX(dm4 -= m2t.transpose() * m3t.transpose(), refMat4 -= refMat2t.transpose() * refMat3t.transpose());
134  VERIFY_IS_APPROX(dm4 = m2 * m3t.transpose(), refMat4 = refMat2 * refMat3t.transpose());
135  VERIFY_IS_APPROX(dm4 += m2 * m3t.transpose(), refMat4 += refMat2 * refMat3t.transpose());
136  VERIFY_IS_APPROX(dm4 -= m2 * m3t.transpose(), refMat4 -= refMat2 * refMat3t.transpose());
137  VERIFY_IS_APPROX(dm4 = m2 * m3 * s1, refMat4 = refMat2 * refMat3 * s1);
138 
139  // test aliasing
140  m4 = m2;
141  refMat4 = refMat2;
142  VERIFY_IS_APPROX(m4 = m4 * m3, refMat4 = refMat4 * refMat3);
143 
144  // sparse * dense matrix
145  VERIFY_IS_APPROX(dm4 = m2 * refMat3, refMat4 = refMat2 * refMat3);
146  VERIFY_IS_APPROX(dm4 = m2 * refMat3t.transpose(), refMat4 = refMat2 * refMat3t.transpose());
147  VERIFY_IS_APPROX(dm4 = m2t.transpose() * refMat3, refMat4 = refMat2t.transpose() * refMat3);
148  VERIFY_IS_APPROX(dm4 = m2t.transpose() * refMat3t.transpose(),
149  refMat4 = refMat2t.transpose() * refMat3t.transpose());
150 
151  VERIFY_IS_APPROX(dm4 = m2 * refMat3, refMat4 = refMat2 * refMat3);
152  VERIFY_IS_APPROX(dm4 = dm4 + m2 * refMat3, refMat4 = refMat4 + refMat2 * refMat3);
153  VERIFY_IS_APPROX(dm4 += m2 * refMat3, refMat4 += refMat2 * refMat3);
154  VERIFY_IS_APPROX(dm4 -= m2 * refMat3, refMat4 -= refMat2 * refMat3);
155  VERIFY_IS_APPROX(dm4.noalias() += m2 * refMat3, refMat4 += refMat2 * refMat3);
156  VERIFY_IS_APPROX(dm4.noalias() -= m2 * refMat3, refMat4 -= refMat2 * refMat3);
157  VERIFY_IS_APPROX(dm4 = m2 * (refMat3 + refMat3), refMat4 = refMat2 * (refMat3 + refMat3));
158  VERIFY_IS_APPROX(dm4 = m2t.transpose() * (refMat3 + refMat5) * 0.5,
159  refMat4 = refMat2t.transpose() * (refMat3 + refMat5) * 0.5);
160 
161  // sparse * dense vector
162  VERIFY_IS_APPROX(dm4.col(0) = m2 * refMat3.col(0), refMat4.col(0) = refMat2 * refMat3.col(0));
163  VERIFY_IS_APPROX(dm4.col(0) = m2 * refMat3t.transpose().col(0),
164  refMat4.col(0) = refMat2 * refMat3t.transpose().col(0));
165  VERIFY_IS_APPROX(dm4.col(0) = m2t.transpose() * refMat3.col(0),
166  refMat4.col(0) = refMat2t.transpose() * refMat3.col(0));
167  VERIFY_IS_APPROX(dm4.col(0) = m2t.transpose() * refMat3t.transpose().col(0),
168  refMat4.col(0) = refMat2t.transpose() * refMat3t.transpose().col(0));
169 
170  // dense * sparse
171  VERIFY_IS_APPROX(dm4 = refMat2 * m3, refMat4 = refMat2 * refMat3);
172  VERIFY_IS_APPROX(dm4 = dm4 + refMat2 * m3, refMat4 = refMat4 + refMat2 * refMat3);
173  VERIFY_IS_APPROX(dm4 += refMat2 * m3, refMat4 += refMat2 * refMat3);
174  VERIFY_IS_APPROX(dm4 -= refMat2 * m3, refMat4 -= refMat2 * refMat3);
175  VERIFY_IS_APPROX(dm4.noalias() += refMat2 * m3, refMat4 += refMat2 * refMat3);
176  VERIFY_IS_APPROX(dm4.noalias() -= refMat2 * m3, refMat4 -= refMat2 * refMat3);
177  VERIFY_IS_APPROX(dm4 = refMat2 * m3t.transpose(), refMat4 = refMat2 * refMat3t.transpose());
178  VERIFY_IS_APPROX(dm4 = refMat2t.transpose() * m3, refMat4 = refMat2t.transpose() * refMat3);
179  VERIFY_IS_APPROX(dm4 = refMat2t.transpose() * m3t.transpose(),
180  refMat4 = refMat2t.transpose() * refMat3t.transpose());
181 
182  // sparse * dense and dense * sparse outer product
183  {
184  Index c = internal::random<Index>(0, depth - 1);
185  Index r = internal::random<Index>(0, rows - 1);
186  Index c1 = internal::random<Index>(0, cols - 1);
187  Index r1 = internal::random<Index>(0, depth - 1);
188  DenseMatrix dm5 = DenseMatrix::Random(depth, cols);
189 
190  VERIFY_IS_APPROX(m4 = m2.col(c) * dm5.col(c1).transpose(), refMat4 = refMat2.col(c) * dm5.col(c1).transpose());
191  VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array() != 0).count());
192  VERIFY_IS_APPROX(m4 = m2.middleCols(c, 1) * dm5.col(c1).transpose(),
193  refMat4 = refMat2.col(c) * dm5.col(c1).transpose());
194  VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array() != 0).count());
195  VERIFY_IS_APPROX(dm4 = m2.col(c) * dm5.col(c1).transpose(), refMat4 = refMat2.col(c) * dm5.col(c1).transpose());
196 
197  VERIFY_IS_APPROX(m4 = dm5.col(c1) * m2.col(c).transpose(), refMat4 = dm5.col(c1) * refMat2.col(c).transpose());
198  VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array() != 0).count());
199  VERIFY_IS_APPROX(m4 = dm5.col(c1) * m2.middleCols(c, 1).transpose(),
200  refMat4 = dm5.col(c1) * refMat2.col(c).transpose());
201  VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array() != 0).count());
202  VERIFY_IS_APPROX(dm4 = dm5.col(c1) * m2.col(c).transpose(), refMat4 = dm5.col(c1) * refMat2.col(c).transpose());
203 
204  VERIFY_IS_APPROX(m4 = dm5.row(r1).transpose() * m2.col(c).transpose(),
205  refMat4 = dm5.row(r1).transpose() * refMat2.col(c).transpose());
206  VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array() != 0).count());
207  VERIFY_IS_APPROX(dm4 = dm5.row(r1).transpose() * m2.col(c).transpose(),
208  refMat4 = dm5.row(r1).transpose() * refMat2.col(c).transpose());
209 
210  VERIFY_IS_APPROX(m4 = m2.row(r).transpose() * dm5.col(c1).transpose(),
211  refMat4 = refMat2.row(r).transpose() * dm5.col(c1).transpose());
212  VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array() != 0).count());
213  VERIFY_IS_APPROX(m4 = m2.middleRows(r, 1).transpose() * dm5.col(c1).transpose(),
214  refMat4 = refMat2.row(r).transpose() * dm5.col(c1).transpose());
215  VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array() != 0).count());
216  VERIFY_IS_APPROX(dm4 = m2.row(r).transpose() * dm5.col(c1).transpose(),
217  refMat4 = refMat2.row(r).transpose() * dm5.col(c1).transpose());
218 
219  VERIFY_IS_APPROX(m4 = dm5.col(c1) * m2.row(r), refMat4 = dm5.col(c1) * refMat2.row(r));
220  VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array() != 0).count());
221  VERIFY_IS_APPROX(m4 = dm5.col(c1) * m2.middleRows(r, 1), refMat4 = dm5.col(c1) * refMat2.row(r));
222  VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array() != 0).count());
223  VERIFY_IS_APPROX(dm4 = dm5.col(c1) * m2.row(r), refMat4 = dm5.col(c1) * refMat2.row(r));
224 
225  VERIFY_IS_APPROX(m4 = dm5.row(r1).transpose() * m2.row(r), refMat4 = dm5.row(r1).transpose() * refMat2.row(r));
226  VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array() != 0).count());
227  VERIFY_IS_APPROX(dm4 = dm5.row(r1).transpose() * m2.row(r), refMat4 = dm5.row(r1).transpose() * refMat2.row(r));
228  }
229 
230  VERIFY_IS_APPROX(m6 = m6 * m6, refMat6 = refMat6 * refMat6);
231 
232  // sparse matrix * sparse vector
233  ColSpVector cv0(cols), cv1;
234  DenseVector dcv0(cols), dcv1;
235  initSparse(2 * density, dcv0, cv0);
236 
237  RowSpVector rv0(depth), rv1;
238  RowDenseVector drv0(depth), drv1(rv1);
239  initSparse(2 * density, drv0, rv0);
240 
241  VERIFY_IS_APPROX(cv1 = m3 * cv0, dcv1 = refMat3 * dcv0);
242  VERIFY_IS_APPROX(rv1 = rv0 * m3, drv1 = drv0 * refMat3);
243  VERIFY_IS_APPROX(cv1 = m3t.adjoint() * cv0, dcv1 = refMat3t.adjoint() * dcv0);
244  VERIFY_IS_APPROX(cv1 = rv0 * m3, dcv1 = drv0 * refMat3);
245  VERIFY_IS_APPROX(rv1 = m3 * cv0, drv1 = refMat3 * dcv0);
246  }
247 
248  // test matrix - diagonal product
249  {
250  DenseMatrix refM2 = DenseMatrix::Zero(rows, cols);
251  DenseMatrix refM3 = DenseMatrix::Zero(rows, cols);
252  DenseMatrix d3 = DenseMatrix::Zero(rows, cols);
253  DiagonalMatrix<Scalar, Dynamic> d1(DenseVector::Random(cols));
254  DiagonalMatrix<Scalar, Dynamic> d2(DenseVector::Random(rows));
255  SparseMatrixType m2(rows, cols);
256  SparseMatrixType m3(rows, cols);
257  initSparse<Scalar>(density, refM2, m2);
258  initSparse<Scalar>(density, refM3, m3);
259  VERIFY_IS_APPROX(m3 = m2 * d1, refM3 = refM2 * d1);
260  VERIFY_IS_APPROX(m3 = m2.transpose() * d2, refM3 = refM2.transpose() * d2);
261  VERIFY_IS_APPROX(m3 = d2 * m2, refM3 = d2 * refM2);
262  VERIFY_IS_APPROX(m3 = d1 * m2.transpose(), refM3 = d1 * refM2.transpose());
263 
264  // also check with a SparseWrapper:
265  DenseVector v1 = DenseVector::Random(cols);
266  DenseVector v2 = DenseVector::Random(rows);
267  DenseVector v3 = DenseVector::Random(rows);
268  VERIFY_IS_APPROX(m3 = m2 * v1.asDiagonal(), refM3 = refM2 * v1.asDiagonal());
269  VERIFY_IS_APPROX(m3 = m2.transpose() * v2.asDiagonal(), refM3 = refM2.transpose() * v2.asDiagonal());
270  VERIFY_IS_APPROX(m3 = v2.asDiagonal() * m2, refM3 = v2.asDiagonal() * refM2);
271  VERIFY_IS_APPROX(m3 = v1.asDiagonal() * m2.transpose(), refM3 = v1.asDiagonal() * refM2.transpose());
272 
273  VERIFY_IS_APPROX(m3 = v2.asDiagonal() * m2 * v1.asDiagonal(), refM3 = v2.asDiagonal() * refM2 * v1.asDiagonal());
274 
275  VERIFY_IS_APPROX(v2 = m2 * v1.asDiagonal() * v1, refM2 * v1.asDiagonal() * v1);
276  VERIFY_IS_APPROX(v3 = v2.asDiagonal() * m2 * v1, v2.asDiagonal() * refM2 * v1);
277 
278  // evaluate to a dense matrix to check the .row() and .col() iterator functions
279  VERIFY_IS_APPROX(d3 = m2 * d1, refM3 = refM2 * d1);
280  VERIFY_IS_APPROX(d3 = m2.transpose() * d2, refM3 = refM2.transpose() * d2);
281  VERIFY_IS_APPROX(d3 = d2 * m2, refM3 = d2 * refM2);
282  VERIFY_IS_APPROX(d3 = d1 * m2.transpose(), refM3 = d1 * refM2.transpose());
283  }
284 
285  // test self-adjoint and triangular-view products
286  {
287  DenseMatrix b = DenseMatrix::Random(rows, rows);
288  DenseMatrix x = DenseMatrix::Random(rows, rows);
289  DenseMatrix refX = DenseMatrix::Random(rows, rows);
290  DenseMatrix refUp = DenseMatrix::Zero(rows, rows);
291  DenseMatrix refLo = DenseMatrix::Zero(rows, rows);
292  DenseMatrix refS = DenseMatrix::Zero(rows, rows);
293  DenseMatrix refA = DenseMatrix::Zero(rows, rows);
294  SparseMatrixType mUp(rows, rows);
295  SparseMatrixType mLo(rows, rows);
296  SparseMatrixType mS(rows, rows);
297  SparseMatrixType mA(rows, rows);
298  initSparse<Scalar>(density, refA, mA);
299  do {
300  initSparse<Scalar>(density, refUp, mUp, ForceRealDiag | /*ForceNonZeroDiag|*/ MakeUpperTriangular);
301  } while (refUp.isZero());
302  refLo = refUp.adjoint();
303  mLo = mUp.adjoint();
304  refS = refUp + refLo;
305  refS.diagonal() *= 0.5;
306  mS = mUp + mLo;
307  // TODO be able to address the diagonal....
308  for (int k = 0; k < mS.outerSize(); ++k)
309  for (typename SparseMatrixType::InnerIterator it(mS, k); it; ++it)
310  if (it.index() == k) it.valueRef() *= Scalar(0.5);
311 
312  VERIFY_IS_APPROX(refS.adjoint(), refS);
313  VERIFY_IS_APPROX(mS.adjoint(), mS);
314  VERIFY_IS_APPROX(mS, refS);
315  VERIFY_IS_APPROX(x = mS * b, refX = refS * b);
316 
317  // sparse selfadjointView with dense matrices
318  VERIFY_IS_APPROX(x = mUp.template selfadjointView<Upper>() * b, refX = refS * b);
319  VERIFY_IS_APPROX(x = mLo.template selfadjointView<Lower>() * b, refX = refS * b);
320  VERIFY_IS_APPROX(x = mS.template selfadjointView<Upper | Lower>() * b, refX = refS * b);
321 
322  VERIFY_IS_APPROX(x = b * mUp.template selfadjointView<Upper>(), refX = b * refS);
323  VERIFY_IS_APPROX(x = b * mLo.template selfadjointView<Lower>(), refX = b * refS);
324  VERIFY_IS_APPROX(x = b * mS.template selfadjointView<Upper | Lower>(), refX = b * refS);
325 
326  VERIFY_IS_APPROX(x.noalias() += mUp.template selfadjointView<Upper>() * b, refX += refS * b);
327  VERIFY_IS_APPROX(x.noalias() -= mLo.template selfadjointView<Lower>() * b, refX -= refS * b);
328  VERIFY_IS_APPROX(x.noalias() += mS.template selfadjointView<Upper | Lower>() * b, refX += refS * b);
329 
330  // sparse selfadjointView with sparse matrices
331  SparseMatrixType mSres(rows, rows);
332  VERIFY_IS_APPROX(mSres = mLo.template selfadjointView<Lower>() * mS,
333  refX = refLo.template selfadjointView<Lower>() * refS);
334  VERIFY_IS_APPROX(mSres = mS * mLo.template selfadjointView<Lower>(),
335  refX = refS * refLo.template selfadjointView<Lower>());
336 
337  // sparse triangularView with dense matrices
338  VERIFY_IS_APPROX(x = mA.template triangularView<Upper>() * b, refX = refA.template triangularView<Upper>() * b);
339  VERIFY_IS_APPROX(x = mA.template triangularView<Lower>() * b, refX = refA.template triangularView<Lower>() * b);
340  VERIFY_IS_APPROX(x = b * mA.template triangularView<Upper>(), refX = b * refA.template triangularView<Upper>());
341  VERIFY_IS_APPROX(x = b * mA.template triangularView<Lower>(), refX = b * refA.template triangularView<Lower>());
342 
343  // sparse triangularView with sparse matrices
344  VERIFY_IS_APPROX(mSres = mA.template triangularView<Lower>() * mS,
345  refX = refA.template triangularView<Lower>() * refS);
346  VERIFY_IS_APPROX(mSres = mS * mA.template triangularView<Lower>(),
347  refX = refS * refA.template triangularView<Lower>());
348  VERIFY_IS_APPROX(mSres = mA.template triangularView<Upper>() * mS,
349  refX = refA.template triangularView<Upper>() * refS);
350  VERIFY_IS_APPROX(mSres = mS * mA.template triangularView<Upper>(),
351  refX = refS * refA.template triangularView<Upper>());
352  }
353 }
DenseMatrix
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Definition: BenchSparseUtil.h:23
DenseVector
Matrix< Scalar, Dynamic, 1 > DenseVector
Definition: BenchSparseUtil.h:24
n
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
m2
MatrixType m2(n_dims)
v2
Map< RowVectorXf > v2(M2.data(), M2.size())
v1
M1<< 1, 2, 3, 4, 5, 6, 7, 8, 9;Map< RowVectorXf > v1(M1.data(), M1.size())
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
b
Scalar * b
Definition: benchVecAdd.cpp:17
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:45
B
Matrix< SCALARB, Dynamic, Dynamic, opt_B > B
Definition: bench_gemm.cpp:48
Eigen::DiagonalMatrix
Represents a diagonal matrix with its storage.
Definition: DiagonalMatrix.h:172
Eigen::SparseVector
a sparse vector class
Definition: SparseVector.h:62
max
#define max(a, b)
Definition: datatypes.h:23
k
char char char int int * k
Definition: level2_impl.h:374
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:367
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:83
Flags
Extend namespace for flags.
Definition: fsi_chan_precond_driver.cc:56
UniformPSDSelfTest.density
density
Definition: UniformPSDSelfTest.py:19
UniformPSDSelfTest.r
r
Definition: UniformPSDSelfTest.py:20
calibrate.c
int c
Definition: calibrate.py:100
oomph::PseudoSolidHelper::Zero
double Zero
Definition: pseudosolid_node_update_elements.cc:35
plotDoE.x
list x
Definition: plotDoE.py:28
ForceRealDiag
@ ForceRealDiag
Definition: sparse.h:32
MakeUpperTriangular
@ MakeUpperTriangular
Definition: sparse.h:32
initSparse
void initSparse(double density, Matrix< Scalar, Dynamic, Dynamic, Opt1 > &refMat, SparseMatrix< Scalar, Opt2, StorageIndex > &sparseMat, int flags=0, std::vector< Matrix< StorageIndex, 2, 1 > > *zeroCoords=0, std::vector< Matrix< StorageIndex, 2, 1 > > *nonzeroCoords=0)
Definition: sparse.h:42
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:47
VERIFY_EVALUATION_COUNT
#define VERIFY_EVALUATION_COUNT(XPR, N)
Definition: test/sparse_product.cpp:28

References b, calibrate::c, cols, UniformPSDSelfTest::density, e(), eval(), ForceRealDiag, initSparse(), k, m2(), MakeUpperTriangular, max, n, UniformPSDSelfTest::r, rows, v1(), v2(), VERIFY_EVALUATION_COUNT, VERIFY_IS_APPROX, VERIFY_IS_EQUAL, plotDoE::x, and oomph::PseudoSolidHelper::Zero.

Referenced by EIGEN_DECLARE_TEST().

◆ sparse_product_regression_test()

template<typename SparseMatrixType , typename DenseMatrixType >
void sparse_product_regression_test ( )
357  {
358  // This code does not compile with afflicted versions of the bug
359  SparseMatrixType sm1(3, 2);
360  DenseMatrixType m2(2, 2);
361  sm1.setZero();
362  m2.setZero();
363 
364  DenseMatrixType m3 = sm1 * m2;
365 
366  // This code produces a segfault with afflicted versions of another SparseTimeDenseProduct
367  // bug
368 
369  SparseMatrixType sm2(20000, 2);
370  sm2.setZero();
371  DenseMatrixType m4(sm2 * m2);
372 
373  VERIFY_IS_APPROX(m4(0, 0), 0.0);
374 }

References m2(), and VERIFY_IS_APPROX.

Referenced by EIGEN_DECLARE_TEST().

◆ test_mixed_storage()

void test_mixed_storage ( )
527  {
528  test_mixed_storage_imp<RowMajor, RowMajor, RowMajor>();
529  test_mixed_storage_imp<RowMajor, RowMajor, ColMajor>();
530  test_mixed_storage_imp<RowMajor, ColMajor, RowMajor>();
531  test_mixed_storage_imp<RowMajor, ColMajor, ColMajor>();
532  test_mixed_storage_imp<ColMajor, RowMajor, RowMajor>();
533  test_mixed_storage_imp<ColMajor, RowMajor, ColMajor>();
534  test_mixed_storage_imp<ColMajor, ColMajor, RowMajor>();
535  test_mixed_storage_imp<ColMajor, ColMajor, ColMajor>();
536 }

Referenced by EIGEN_DECLARE_TEST().

◆ test_mixed_storage_imp()

template<int OrderA, int OrderB, int OrderC>
void test_mixed_storage_imp ( )
488  {
489  typedef float Real;
490  typedef Matrix<Real, Dynamic, Dynamic> DenseMat;
491 
492  // Case: Large inputs but small result
493  {
494  SparseMatrix<Real, OrderA> A(8, 512);
495  SparseMatrix<Real, OrderB> B(512, 8);
496  DenseMat refA(8, 512);
497  DenseMat refB(512, 8);
498 
499  initSparse<Real>(0.1, refA, A);
500  initSparse<Real>(0.1, refB, B);
501 
502  SparseMatrix<Real, OrderC, std::int8_t> result;
503  SparseMatrix<Real, OrderC> result_large;
504  DenseMat refResult;
505 
506  VERIFY_IS_APPROX(result = (A * B), refResult = refA * refB);
507  }
508 
509  // Case: Small input but large result
510  {
511  SparseMatrix<Real, OrderA, std::int8_t> A(127, 8);
512  SparseMatrix<Real, OrderB, std::int8_t> B(8, 127);
513  DenseMat refA(127, 8);
514  DenseMat refB(8, 127);
515 
516  initSparse<Real>(0.01, refA, A);
517  initSparse<Real>(0.01, refB, B);
518 
519  SparseMatrix<Real, OrderC> result;
520  SparseMatrix<Real, OrderC> result_large;
521  DenseMat refResult;
522 
523  VERIFY_IS_APPROX(result = (A * B), refResult = refA * refB);
524  }
525 }
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:47
Real
boost::multiprecision::number< boost::multiprecision::cpp_dec_float< 100 >, boost::multiprecision::et_on > Real
Definition: boostmultiprec.cpp:77
oomph::Matrix
Definition: matrices.h:74

References VERIFY_IS_APPROX.

◆ test_mixing_types()

template<typename Real >
void test_mixing_types ( )
399  {
400  typedef std::complex<Real> Cplx;
401  typedef SparseMatrix<Real> SpMatReal;
402  typedef SparseMatrix<Cplx> SpMatCplx;
403  typedef SparseMatrix<Cplx, RowMajor> SpRowMatCplx;
404  typedef Matrix<Real, Dynamic, Dynamic> DenseMatReal;
405  typedef Matrix<Cplx, Dynamic, Dynamic> DenseMatCplx;
406 
407  Index n = internal::random<Index>(1, 100);
408  double density = (std::max)(8. / static_cast<double>(n * n), 0.2);
409 
410  SpMatReal sR1(n, n);
411  SpMatCplx sC1(n, n), sC2(n, n), sC3(n, n);
412  SpRowMatCplx sCR(n, n);
413  DenseMatReal dR1(n, n);
414  DenseMatCplx dC1(n, n), dC2(n, n), dC3(n, n);
415 
416  initSparse<Real>(density, dR1, sR1);
417  initSparse<Cplx>(density, dC1, sC1);
418  initSparse<Cplx>(density, dC2, sC2);
419 
420  VERIFY_IS_APPROX(sC2 = (sR1 * sC1), dC3 = dR1.template cast<Cplx>() * dC1);
421  VERIFY_IS_APPROX(sC2 = (sC1 * sR1), dC3 = dC1 * dR1.template cast<Cplx>());
422  VERIFY_IS_APPROX(sC2 = (sR1.transpose() * sC1), dC3 = dR1.template cast<Cplx>().transpose() * dC1);
423  VERIFY_IS_APPROX(sC2 = (sC1.transpose() * sR1), dC3 = dC1.transpose() * dR1.template cast<Cplx>());
424  VERIFY_IS_APPROX(sC2 = (sR1 * sC1.transpose()), dC3 = dR1.template cast<Cplx>() * dC1.transpose());
425  VERIFY_IS_APPROX(sC2 = (sC1 * sR1.transpose()), dC3 = dC1 * dR1.template cast<Cplx>().transpose());
426  VERIFY_IS_APPROX(sC2 = (sR1.transpose() * sC1.transpose()),
427  dC3 = dR1.template cast<Cplx>().transpose() * dC1.transpose());
428  VERIFY_IS_APPROX(sC2 = (sC1.transpose() * sR1.transpose()),
429  dC3 = dC1.transpose() * dR1.template cast<Cplx>().transpose());
430 
431  VERIFY_IS_APPROX(sCR = (sR1 * sC1), dC3 = dR1.template cast<Cplx>() * dC1);
432  VERIFY_IS_APPROX(sCR = (sC1 * sR1), dC3 = dC1 * dR1.template cast<Cplx>());
433  VERIFY_IS_APPROX(sCR = (sR1.transpose() * sC1), dC3 = dR1.template cast<Cplx>().transpose() * dC1);
434  VERIFY_IS_APPROX(sCR = (sC1.transpose() * sR1), dC3 = dC1.transpose() * dR1.template cast<Cplx>());
435  VERIFY_IS_APPROX(sCR = (sR1 * sC1.transpose()), dC3 = dR1.template cast<Cplx>() * dC1.transpose());
436  VERIFY_IS_APPROX(sCR = (sC1 * sR1.transpose()), dC3 = dC1 * dR1.template cast<Cplx>().transpose());
437  VERIFY_IS_APPROX(sCR = (sR1.transpose() * sC1.transpose()),
438  dC3 = dR1.template cast<Cplx>().transpose() * dC1.transpose());
439  VERIFY_IS_APPROX(sCR = (sC1.transpose() * sR1.transpose()),
440  dC3 = dC1.transpose() * dR1.template cast<Cplx>().transpose());
441 
442  VERIFY_IS_APPROX(sC2 = (sR1 * sC1).pruned(), dC3 = dR1.template cast<Cplx>() * dC1);
443  VERIFY_IS_APPROX(sC2 = (sC1 * sR1).pruned(), dC3 = dC1 * dR1.template cast<Cplx>());
444  VERIFY_IS_APPROX(sC2 = (sR1.transpose() * sC1).pruned(), dC3 = dR1.template cast<Cplx>().transpose() * dC1);
445  VERIFY_IS_APPROX(sC2 = (sC1.transpose() * sR1).pruned(), dC3 = dC1.transpose() * dR1.template cast<Cplx>());
446  VERIFY_IS_APPROX(sC2 = (sR1 * sC1.transpose()).pruned(), dC3 = dR1.template cast<Cplx>() * dC1.transpose());
447  VERIFY_IS_APPROX(sC2 = (sC1 * sR1.transpose()).pruned(), dC3 = dC1 * dR1.template cast<Cplx>().transpose());
448  VERIFY_IS_APPROX(sC2 = (sR1.transpose() * sC1.transpose()).pruned(),
449  dC3 = dR1.template cast<Cplx>().transpose() * dC1.transpose());
450  VERIFY_IS_APPROX(sC2 = (sC1.transpose() * sR1.transpose()).pruned(),
451  dC3 = dC1.transpose() * dR1.template cast<Cplx>().transpose());
452 
453  VERIFY_IS_APPROX(sCR = (sR1 * sC1).pruned(), dC3 = dR1.template cast<Cplx>() * dC1);
454  VERIFY_IS_APPROX(sCR = (sC1 * sR1).pruned(), dC3 = dC1 * dR1.template cast<Cplx>());
455  VERIFY_IS_APPROX(sCR = (sR1.transpose() * sC1).pruned(), dC3 = dR1.template cast<Cplx>().transpose() * dC1);
456  VERIFY_IS_APPROX(sCR = (sC1.transpose() * sR1).pruned(), dC3 = dC1.transpose() * dR1.template cast<Cplx>());
457  VERIFY_IS_APPROX(sCR = (sR1 * sC1.transpose()).pruned(), dC3 = dR1.template cast<Cplx>() * dC1.transpose());
458  VERIFY_IS_APPROX(sCR = (sC1 * sR1.transpose()).pruned(), dC3 = dC1 * dR1.template cast<Cplx>().transpose());
459  VERIFY_IS_APPROX(sCR = (sR1.transpose() * sC1.transpose()).pruned(),
460  dC3 = dR1.template cast<Cplx>().transpose() * dC1.transpose());
461  VERIFY_IS_APPROX(sCR = (sC1.transpose() * sR1.transpose()).pruned(),
462  dC3 = dC1.transpose() * dR1.template cast<Cplx>().transpose());
463 
464  VERIFY_IS_APPROX(dC2 = (sR1 * sC1), dC3 = dR1.template cast<Cplx>() * dC1);
465  VERIFY_IS_APPROX(dC2 = (sC1 * sR1), dC3 = dC1 * dR1.template cast<Cplx>());
466  VERIFY_IS_APPROX(dC2 = (sR1.transpose() * sC1), dC3 = dR1.template cast<Cplx>().transpose() * dC1);
467  VERIFY_IS_APPROX(dC2 = (sC1.transpose() * sR1), dC3 = dC1.transpose() * dR1.template cast<Cplx>());
468  VERIFY_IS_APPROX(dC2 = (sR1 * sC1.transpose()), dC3 = dR1.template cast<Cplx>() * dC1.transpose());
469  VERIFY_IS_APPROX(dC2 = (sC1 * sR1.transpose()), dC3 = dC1 * dR1.template cast<Cplx>().transpose());
470  VERIFY_IS_APPROX(dC2 = (sR1.transpose() * sC1.transpose()),
471  dC3 = dR1.template cast<Cplx>().transpose() * dC1.transpose());
472  VERIFY_IS_APPROX(dC2 = (sC1.transpose() * sR1.transpose()),
473  dC3 = dC1.transpose() * dR1.template cast<Cplx>().transpose());
474 
475  VERIFY_IS_APPROX(dC2 = dR1 * sC1, dC3 = dR1.template cast<Cplx>() * sC1);
476  VERIFY_IS_APPROX(dC2 = sR1 * dC1, dC3 = sR1.template cast<Cplx>() * dC1);
477  VERIFY_IS_APPROX(dC2 = dC1 * sR1, dC3 = dC1 * sR1.template cast<Cplx>());
478  VERIFY_IS_APPROX(dC2 = sC1 * dR1, dC3 = sC1 * dR1.template cast<Cplx>());
479 
480  VERIFY_IS_APPROX(dC2 = dR1.row(0) * sC1, dC3 = dR1.template cast<Cplx>().row(0) * sC1);
481  VERIFY_IS_APPROX(dC2 = sR1 * dC1.col(0), dC3 = sR1.template cast<Cplx>() * dC1.col(0));
482  VERIFY_IS_APPROX(dC2 = dC1.row(0) * sR1, dC3 = dC1.row(0) * sR1.template cast<Cplx>());
483  VERIFY_IS_APPROX(dC2 = sC1 * dR1.col(0), dC3 = sC1 * dR1.template cast<Cplx>().col(0));
484 }

References UniformPSDSelfTest::density, max, n, and VERIFY_IS_APPROX.

Variable Documentation

◆ nb_temporaries

long int nb_temporaries
static

Referenced by on_temporary_creation().