sparse_basic.cpp
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
6 // Copyright (C) 2013 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
7 //
8 // This Source Code Form is subject to the terms of the Mozilla
9 // Public License v. 2.0. If a copy of the MPL was not distributed
10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11 
12 static long g_realloc_count = 0;
13 #define EIGEN_SPARSE_COMPRESSED_STORAGE_REALLOCATE_PLUGIN g_realloc_count++;
14 
15 #include "sparse.h"
16 
17 template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref)
18 {
19  typedef typename SparseMatrixType::StorageIndex StorageIndex;
20  typedef Matrix<StorageIndex,2,1> Vector2;
21 
22  const Index rows = ref.rows();
23  const Index cols = ref.cols();
24  //const Index inner = ref.innerSize();
25  //const Index outer = ref.outerSize();
26 
27  typedef typename SparseMatrixType::Scalar Scalar;
28  typedef typename SparseMatrixType::RealScalar RealScalar;
29  enum { Flags = SparseMatrixType::Flags };
30 
31  double density = (std::max)(8./(rows*cols), 0.01);
32  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
33  typedef Matrix<Scalar,Dynamic,1> DenseVector;
34  Scalar eps = 1e-6;
35 
36  Scalar s1 = internal::random<Scalar>();
37  {
38  SparseMatrixType m(rows, cols);
39  DenseMatrix refMat = DenseMatrix::Zero(rows, cols);
40  DenseVector vec1 = DenseVector::Random(rows);
41 
42  std::vector<Vector2> zeroCoords;
43  std::vector<Vector2> nonzeroCoords;
44  initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords);
45 
46  // test coeff and coeffRef
47  for (std::size_t i=0; i<zeroCoords.size(); ++i)
48  {
49  VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps );
50  if(internal::is_same<SparseMatrixType,SparseMatrix<Scalar,Flags> >::value)
51  VERIFY_RAISES_ASSERT( m.coeffRef(zeroCoords[i].x(),zeroCoords[i].y()) = 5 );
52  }
53  VERIFY_IS_APPROX(m, refMat);
54 
55  if(!nonzeroCoords.empty()) {
56  m.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
57  refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
58  }
59 
60  VERIFY_IS_APPROX(m, refMat);
61 
62  // test assertion
63  VERIFY_RAISES_ASSERT( m.coeffRef(-1,1) = 0 );
64  VERIFY_RAISES_ASSERT( m.coeffRef(0,m.cols()) = 0 );
65  }
66 
67  // test insert (inner random)
68  {
69  DenseMatrix m1(rows,cols);
70  m1.setZero();
71  SparseMatrixType m2(rows,cols);
72  bool call_reserve = internal::random<int>()%2;
73  Index nnz = internal::random<int>(1,int(rows)/2);
74  if(call_reserve)
75  {
76  if(internal::random<int>()%2)
77  m2.reserve(VectorXi::Constant(m2.outerSize(), int(nnz)));
78  else
79  m2.reserve(m2.outerSize() * nnz);
80  }
81  g_realloc_count = 0;
82  for (Index j=0; j<cols; ++j)
83  {
84  for (Index k=0; k<nnz; ++k)
85  {
86  Index i = internal::random<Index>(0,rows-1);
87  if (m1.coeff(i,j)==Scalar(0))
88  m2.insert(i,j) = m1(i,j) = internal::random<Scalar>();
89  }
90  }
91 
92  if(call_reserve && !SparseMatrixType::IsRowMajor)
93  {
94  VERIFY(g_realloc_count==0);
95  }
96 
97  m2.finalize();
98  VERIFY_IS_APPROX(m2,m1);
99  }
100 
101  // test insert (fully random)
102  {
103  DenseMatrix m1(rows,cols);
104  m1.setZero();
105  SparseMatrixType m2(rows,cols);
106  if(internal::random<int>()%2)
107  m2.reserve(VectorXi::Constant(m2.outerSize(), 2));
108  for (int k=0; k<rows*cols; ++k)
109  {
110  Index i = internal::random<Index>(0,rows-1);
111  Index j = internal::random<Index>(0,cols-1);
112  if ((m1.coeff(i,j)==Scalar(0)) && (internal::random<int>()%2))
113  m2.insert(i,j) = m1(i,j) = internal::random<Scalar>();
114  else
115  {
116  Scalar v = internal::random<Scalar>();
117  m2.coeffRef(i,j) += v;
118  m1(i,j) += v;
119  }
120  }
121  VERIFY_IS_APPROX(m2,m1);
122  }
123 
124  // test insert (un-compressed)
125  for(int mode=0;mode<4;++mode)
126  {
127  DenseMatrix m1(rows,cols);
128  m1.setZero();
129  SparseMatrixType m2(rows,cols);
130  VectorXi r(VectorXi::Constant(m2.outerSize(), ((mode%2)==0) ? int(m2.innerSize()) : std::max<int>(1,int(m2.innerSize())/8)));
131  m2.reserve(r);
132  for (Index k=0; k<rows*cols; ++k)
133  {
134  Index i = internal::random<Index>(0,rows-1);
135  Index j = internal::random<Index>(0,cols-1);
136  if (m1.coeff(i,j)==Scalar(0))
137  m2.insert(i,j) = m1(i,j) = internal::random<Scalar>();
138  if(mode==3)
139  m2.reserve(r);
140  }
141  if(internal::random<int>()%2)
142  m2.makeCompressed();
143  VERIFY_IS_APPROX(m2,m1);
144  }
145 
146  // test basic computations
147  {
148  DenseMatrix refM1 = DenseMatrix::Zero(rows, cols);
149  DenseMatrix refM2 = DenseMatrix::Zero(rows, cols);
150  DenseMatrix refM3 = DenseMatrix::Zero(rows, cols);
151  DenseMatrix refM4 = DenseMatrix::Zero(rows, cols);
152  SparseMatrixType m1(rows, cols);
153  SparseMatrixType m2(rows, cols);
154  SparseMatrixType m3(rows, cols);
155  SparseMatrixType m4(rows, cols);
156  initSparse<Scalar>(density, refM1, m1);
157  initSparse<Scalar>(density, refM2, m2);
158  initSparse<Scalar>(density, refM3, m3);
159  initSparse<Scalar>(density, refM4, m4);
160 
161  if(internal::random<bool>())
162  m1.makeCompressed();
163 
164  Index m1_nnz = m1.nonZeros();
165 
166  VERIFY_IS_APPROX(m1*s1, refM1*s1);
167  VERIFY_IS_APPROX(m1+m2, refM1+refM2);
168  VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3);
169  VERIFY_IS_APPROX(m3.cwiseProduct(m1+m2), refM3.cwiseProduct(refM1+refM2));
170  VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2);
171  VERIFY_IS_APPROX(m4=m1/s1, refM1/s1);
172  VERIFY_IS_EQUAL(m4.nonZeros(), m1_nnz);
173 
174  if(SparseMatrixType::IsRowMajor)
175  VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.row(0).dot(refM2.row(0)));
176  else
177  VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.col(0)), refM1.col(0).dot(refM2.col(0)));
178 
179  DenseVector rv = DenseVector::Random(m1.cols());
180  DenseVector cv = DenseVector::Random(m1.rows());
181  Index r = internal::random<Index>(0,m1.rows()-2);
182  Index c = internal::random<Index>(0,m1.cols()-1);
183  VERIFY_IS_APPROX(( m1.template block<1,Dynamic>(r,0,1,m1.cols()).dot(rv)) , refM1.row(r).dot(rv));
184  VERIFY_IS_APPROX(m1.row(r).dot(rv), refM1.row(r).dot(rv));
185  VERIFY_IS_APPROX(m1.col(c).dot(cv), refM1.col(c).dot(cv));
186 
187  VERIFY_IS_APPROX(m1.conjugate(), refM1.conjugate());
188  VERIFY_IS_APPROX(m1.real(), refM1.real());
189 
190  refM4.setRandom();
191  // sparse cwise* dense
192  VERIFY_IS_APPROX(m3.cwiseProduct(refM4), refM3.cwiseProduct(refM4));
193  // dense cwise* sparse
194  VERIFY_IS_APPROX(refM4.cwiseProduct(m3), refM4.cwiseProduct(refM3));
195 // VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4);
196 
197  VERIFY_IS_APPROX(refM4 + m3, refM4 + refM3);
198  VERIFY_IS_APPROX(m3 + refM4, refM3 + refM4);
199  VERIFY_IS_APPROX(refM4 - m3, refM4 - refM3);
200  VERIFY_IS_APPROX(m3 - refM4, refM3 - refM4);
201  VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + RealScalar(0.5)*m3).eval(), RealScalar(0.5)*refM4 + RealScalar(0.5)*refM3);
202  VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + m3*RealScalar(0.5)).eval(), RealScalar(0.5)*refM4 + RealScalar(0.5)*refM3);
203  VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + m3.cwiseProduct(m3)).eval(), RealScalar(0.5)*refM4 + refM3.cwiseProduct(refM3));
204 
205  VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + RealScalar(0.5)*m3).eval(), RealScalar(0.5)*refM4 + RealScalar(0.5)*refM3);
206  VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + m3*RealScalar(0.5)).eval(), RealScalar(0.5)*refM4 + RealScalar(0.5)*refM3);
207  VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + (m3+m3)).eval(), RealScalar(0.5)*refM4 + (refM3+refM3));
208  VERIFY_IS_APPROX(((refM3+m3)+RealScalar(0.5)*m3).eval(), RealScalar(0.5)*refM3 + (refM3+refM3));
209  VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + (refM3+m3)).eval(), RealScalar(0.5)*refM4 + (refM3+refM3));
210  VERIFY_IS_APPROX((RealScalar(0.5)*refM4 + (m3+refM3)).eval(), RealScalar(0.5)*refM4 + (refM3+refM3));
211 
212 
213  VERIFY_IS_APPROX(m1.sum(), refM1.sum());
214 
215  m4 = m1; refM4 = m4;
216 
217  VERIFY_IS_APPROX(m1*=s1, refM1*=s1);
218  VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz);
219  VERIFY_IS_APPROX(m1/=s1, refM1/=s1);
220  VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz);
221 
222  VERIFY_IS_APPROX(m1+=m2, refM1+=refM2);
223  VERIFY_IS_APPROX(m1-=m2, refM1-=refM2);
224 
225  if (rows>=2 && cols>=2)
226  {
227  VERIFY_RAISES_ASSERT( m1 += m1.innerVector(0) );
228  VERIFY_RAISES_ASSERT( m1 -= m1.innerVector(0) );
229  VERIFY_RAISES_ASSERT( refM1 -= m1.innerVector(0) );
230  VERIFY_RAISES_ASSERT( refM1 += m1.innerVector(0) );
231  }
232  m1 = m4; refM1 = refM4;
233 
234  // test aliasing
235  VERIFY_IS_APPROX((m1 = -m1), (refM1 = -refM1));
236  VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz);
237  m1 = m4; refM1 = refM4;
238  VERIFY_IS_APPROX((m1 = m1.transpose()), (refM1 = refM1.transpose().eval()));
239  VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz);
240  m1 = m4; refM1 = refM4;
241  VERIFY_IS_APPROX((m1 = -m1.transpose()), (refM1 = -refM1.transpose().eval()));
242  VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz);
243  m1 = m4; refM1 = refM4;
244  VERIFY_IS_APPROX((m1 += -m1), (refM1 += -refM1));
245  VERIFY_IS_EQUAL(m1.nonZeros(), m1_nnz);
246  m1 = m4; refM1 = refM4;
247 
248  if(m1.isCompressed())
249  {
250  VERIFY_IS_APPROX(m1.coeffs().sum(), m1.sum());
251  m1.coeffs() += s1;
252  for(Index j = 0; j<m1.outerSize(); ++j)
253  for(typename SparseMatrixType::InnerIterator it(m1,j); it; ++it)
254  refM1(it.row(), it.col()) += s1;
255  VERIFY_IS_APPROX(m1, refM1);
256  }
257 
258  // and/or
259  {
260  typedef SparseMatrix<bool, SparseMatrixType::Options, typename SparseMatrixType::StorageIndex> SpBool;
261  SpBool mb1 = m1.real().template cast<bool>();
262  SpBool mb2 = m2.real().template cast<bool>();
263  VERIFY_IS_EQUAL(mb1.template cast<int>().sum(), refM1.real().template cast<bool>().count());
264  VERIFY_IS_EQUAL((mb1 && mb2).template cast<int>().sum(), (refM1.real().template cast<bool>() && refM2.real().template cast<bool>()).count());
265  VERIFY_IS_EQUAL((mb1 || mb2).template cast<int>().sum(), (refM1.real().template cast<bool>() || refM2.real().template cast<bool>()).count());
266  SpBool mb3 = mb1 && mb2;
267  if(mb1.coeffs().all() && mb2.coeffs().all())
268  {
269  VERIFY_IS_EQUAL(mb3.nonZeros(), (refM1.real().template cast<bool>() && refM2.real().template cast<bool>()).count());
270  }
271  }
272  }
273 
274  // test reverse iterators
275  {
276  DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
277  SparseMatrixType m2(rows, cols);
278  initSparse<Scalar>(density, refMat2, m2);
279  std::vector<Scalar> ref_value(m2.innerSize());
280  std::vector<Index> ref_index(m2.innerSize());
281  if(internal::random<bool>())
282  m2.makeCompressed();
283  for(Index j = 0; j<m2.outerSize(); ++j)
284  {
285  Index count_forward = 0;
286 
287  for(typename SparseMatrixType::InnerIterator it(m2,j); it; ++it)
288  {
289  ref_value[ref_value.size()-1-count_forward] = it.value();
290  ref_index[ref_index.size()-1-count_forward] = it.index();
291  count_forward++;
292  }
293  Index count_reverse = 0;
294  for(typename SparseMatrixType::ReverseInnerIterator it(m2,j); it; --it)
295  {
296  VERIFY_IS_APPROX( std::abs(ref_value[ref_value.size()-count_forward+count_reverse])+1, std::abs(it.value())+1);
297  VERIFY_IS_EQUAL( ref_index[ref_index.size()-count_forward+count_reverse] , it.index());
298  count_reverse++;
299  }
300  VERIFY_IS_EQUAL(count_forward, count_reverse);
301  }
302  }
303 
304  // test transpose
305  {
306  DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
307  SparseMatrixType m2(rows, cols);
308  initSparse<Scalar>(density, refMat2, m2);
309  VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval());
310  VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose());
311 
312  VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint());
313 
314  // check isApprox handles opposite storage order
315  typename Transpose<SparseMatrixType>::PlainObject m3(m2);
316  VERIFY(m2.isApprox(m3));
317  }
318 
319  // test prune
320  {
321  SparseMatrixType m2(rows, cols);
322  DenseMatrix refM2(rows, cols);
323  refM2.setZero();
324  int countFalseNonZero = 0;
325  int countTrueNonZero = 0;
326  m2.reserve(VectorXi::Constant(m2.outerSize(), int(m2.innerSize())));
327  for (Index j=0; j<m2.cols(); ++j)
328  {
329  for (Index i=0; i<m2.rows(); ++i)
330  {
331  float x = internal::random<float>(0,1);
332  if (x<0.1f)
333  {
334  // do nothing
335  }
336  else if (x<0.5f)
337  {
338  countFalseNonZero++;
339  m2.insert(i,j) = Scalar(0);
340  }
341  else
342  {
343  countTrueNonZero++;
344  m2.insert(i,j) = Scalar(1);
345  refM2(i,j) = Scalar(1);
346  }
347  }
348  }
349  if(internal::random<bool>())
350  m2.makeCompressed();
351  VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros());
352  if(countTrueNonZero>0)
353  VERIFY_IS_APPROX(m2, refM2);
354  m2.prune(Scalar(1));
355  VERIFY(countTrueNonZero==m2.nonZeros());
356  VERIFY_IS_APPROX(m2, refM2);
357  }
358 
359  // test setFromTriplets
360  {
361  typedef Triplet<Scalar,StorageIndex> TripletType;
362  std::vector<TripletType> triplets;
363  Index ntriplets = rows*cols;
364  triplets.reserve(ntriplets);
365  DenseMatrix refMat_sum = DenseMatrix::Zero(rows,cols);
366  DenseMatrix refMat_prod = DenseMatrix::Zero(rows,cols);
367  DenseMatrix refMat_last = DenseMatrix::Zero(rows,cols);
368 
369  for(Index i=0;i<ntriplets;++i)
370  {
371  StorageIndex r = internal::random<StorageIndex>(0,StorageIndex(rows-1));
372  StorageIndex c = internal::random<StorageIndex>(0,StorageIndex(cols-1));
373  Scalar v = internal::random<Scalar>();
374  triplets.push_back(TripletType(r,c,v));
375  refMat_sum(r,c) += v;
376  if(std::abs(refMat_prod(r,c))==0)
377  refMat_prod(r,c) = v;
378  else
379  refMat_prod(r,c) *= v;
380  refMat_last(r,c) = v;
381  }
382  SparseMatrixType m(rows,cols);
383  m.setFromTriplets(triplets.begin(), triplets.end());
384  VERIFY_IS_APPROX(m, refMat_sum);
385 
386  m.setFromTriplets(triplets.begin(), triplets.end(), std::multiplies<Scalar>());
387  VERIFY_IS_APPROX(m, refMat_prod);
388 #if (defined(__cplusplus) && __cplusplus >= 201103L)
389  m.setFromTriplets(triplets.begin(), triplets.end(), [] (Scalar,Scalar b) { return b; });
390  VERIFY_IS_APPROX(m, refMat_last);
391 #endif
392  }
393 
394  // test Map
395  {
396  DenseMatrix refMat2(rows, cols), refMat3(rows, cols);
397  SparseMatrixType m2(rows, cols), m3(rows, cols);
398  initSparse<Scalar>(density, refMat2, m2);
399  initSparse<Scalar>(density, refMat3, m3);
400  {
401  Map<SparseMatrixType> mapMat2(m2.rows(), m2.cols(), m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr());
402  Map<SparseMatrixType> mapMat3(m3.rows(), m3.cols(), m3.nonZeros(), m3.outerIndexPtr(), m3.innerIndexPtr(), m3.valuePtr(), m3.innerNonZeroPtr());
403  VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
404  VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
405  }
406  {
407  MappedSparseMatrix<Scalar,SparseMatrixType::Options,StorageIndex> mapMat2(m2.rows(), m2.cols(), m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr());
408  MappedSparseMatrix<Scalar,SparseMatrixType::Options,StorageIndex> mapMat3(m3.rows(), m3.cols(), m3.nonZeros(), m3.outerIndexPtr(), m3.innerIndexPtr(), m3.valuePtr(), m3.innerNonZeroPtr());
409  VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
410  VERIFY_IS_APPROX(mapMat2+mapMat3, refMat2+refMat3);
411  }
412 
413  Index i = internal::random<Index>(0,rows-1);
414  Index j = internal::random<Index>(0,cols-1);
415  m2.coeffRef(i,j) = 123;
416  if(internal::random<bool>())
417  m2.makeCompressed();
418  Map<SparseMatrixType> mapMat2(rows, cols, m2.nonZeros(), m2.outerIndexPtr(), m2.innerIndexPtr(), m2.valuePtr(), m2.innerNonZeroPtr());
419  VERIFY_IS_EQUAL(m2.coeff(i,j),Scalar(123));
420  VERIFY_IS_EQUAL(mapMat2.coeff(i,j),Scalar(123));
421  mapMat2.coeffRef(i,j) = -123;
422  VERIFY_IS_EQUAL(m2.coeff(i,j),Scalar(-123));
423  }
424 
425  // test triangularView
426  {
427  DenseMatrix refMat2(rows, cols), refMat3(rows, cols);
428  SparseMatrixType m2(rows, cols), m3(rows, cols);
429  initSparse<Scalar>(density, refMat2, m2);
430  refMat3 = refMat2.template triangularView<Lower>();
431  m3 = m2.template triangularView<Lower>();
432  VERIFY_IS_APPROX(m3, refMat3);
433 
434  refMat3 = refMat2.template triangularView<Upper>();
435  m3 = m2.template triangularView<Upper>();
436  VERIFY_IS_APPROX(m3, refMat3);
437 
438  {
439  refMat3 = refMat2.template triangularView<UnitUpper>();
440  m3 = m2.template triangularView<UnitUpper>();
441  VERIFY_IS_APPROX(m3, refMat3);
442 
443  refMat3 = refMat2.template triangularView<UnitLower>();
444  m3 = m2.template triangularView<UnitLower>();
445  VERIFY_IS_APPROX(m3, refMat3);
446  }
447 
448  refMat3 = refMat2.template triangularView<StrictlyUpper>();
449  m3 = m2.template triangularView<StrictlyUpper>();
450  VERIFY_IS_APPROX(m3, refMat3);
451 
452  refMat3 = refMat2.template triangularView<StrictlyLower>();
453  m3 = m2.template triangularView<StrictlyLower>();
454  VERIFY_IS_APPROX(m3, refMat3);
455 
456  // check sparse-triangular to dense
457  refMat3 = m2.template triangularView<StrictlyUpper>();
458  VERIFY_IS_APPROX(refMat3, DenseMatrix(refMat2.template triangularView<StrictlyUpper>()));
459  }
460 
461  // test selfadjointView
462  if(!SparseMatrixType::IsRowMajor)
463  {
464  DenseMatrix refMat2(rows, rows), refMat3(rows, rows);
465  SparseMatrixType m2(rows, rows), m3(rows, rows);
466  initSparse<Scalar>(density, refMat2, m2);
467  refMat3 = refMat2.template selfadjointView<Lower>();
468  m3 = m2.template selfadjointView<Lower>();
469  VERIFY_IS_APPROX(m3, refMat3);
470 
471  refMat3 += refMat2.template selfadjointView<Lower>();
472  m3 += m2.template selfadjointView<Lower>();
473  VERIFY_IS_APPROX(m3, refMat3);
474 
475  refMat3 -= refMat2.template selfadjointView<Lower>();
476  m3 -= m2.template selfadjointView<Lower>();
477  VERIFY_IS_APPROX(m3, refMat3);
478 
479  // selfadjointView only works for square matrices:
480  SparseMatrixType m4(rows, rows+1);
481  VERIFY_RAISES_ASSERT(m4.template selfadjointView<Lower>());
482  VERIFY_RAISES_ASSERT(m4.template selfadjointView<Upper>());
483  }
484 
485  // test sparseView
486  {
487  DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
488  SparseMatrixType m2(rows, rows);
489  initSparse<Scalar>(density, refMat2, m2);
490  VERIFY_IS_APPROX(m2.eval(), refMat2.sparseView().eval());
491 
492  // sparse view on expressions:
493  VERIFY_IS_APPROX((s1*m2).eval(), (s1*refMat2).sparseView().eval());
494  VERIFY_IS_APPROX((m2+m2).eval(), (refMat2+refMat2).sparseView().eval());
495  VERIFY_IS_APPROX((m2*m2).eval(), (refMat2.lazyProduct(refMat2)).sparseView().eval());
496  VERIFY_IS_APPROX((m2*m2).eval(), (refMat2*refMat2).sparseView().eval());
497  }
498 
499  // test diagonal
500  {
501  DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
502  SparseMatrixType m2(rows, cols);
503  initSparse<Scalar>(density, refMat2, m2);
504  VERIFY_IS_APPROX(m2.diagonal(), refMat2.diagonal().eval());
505  DenseVector d = m2.diagonal();
506  VERIFY_IS_APPROX(d, refMat2.diagonal().eval());
507  d = m2.diagonal().array();
508  VERIFY_IS_APPROX(d, refMat2.diagonal().eval());
509  VERIFY_IS_APPROX(const_cast<const SparseMatrixType&>(m2).diagonal(), refMat2.diagonal().eval());
510 
511  initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag);
512  m2.diagonal() += refMat2.diagonal();
513  refMat2.diagonal() += refMat2.diagonal();
514  VERIFY_IS_APPROX(m2, refMat2);
515  }
516 
517  // test diagonal to sparse
518  {
519  DenseVector d = DenseVector::Random(rows);
520  DenseMatrix refMat2 = d.asDiagonal();
521  SparseMatrixType m2(rows, rows);
522  m2 = d.asDiagonal();
523  VERIFY_IS_APPROX(m2, refMat2);
524  SparseMatrixType m3(d.asDiagonal());
525  VERIFY_IS_APPROX(m3, refMat2);
526  refMat2 += d.asDiagonal();
527  m2 += d.asDiagonal();
528  VERIFY_IS_APPROX(m2, refMat2);
529  }
530 
531  // test conservative resize
532  {
533  std::vector< std::pair<StorageIndex,StorageIndex> > inc;
534  if(rows > 3 && cols > 2)
535  inc.push_back(std::pair<StorageIndex,StorageIndex>(-3,-2));
536  inc.push_back(std::pair<StorageIndex,StorageIndex>(0,0));
537  inc.push_back(std::pair<StorageIndex,StorageIndex>(3,2));
538  inc.push_back(std::pair<StorageIndex,StorageIndex>(3,0));
539  inc.push_back(std::pair<StorageIndex,StorageIndex>(0,3));
540 
541  for(size_t i = 0; i< inc.size(); i++) {
542  StorageIndex incRows = inc[i].first;
543  StorageIndex incCols = inc[i].second;
544  SparseMatrixType m1(rows, cols);
545  DenseMatrix refMat1 = DenseMatrix::Zero(rows, cols);
546  initSparse<Scalar>(density, refMat1, m1);
547 
548  m1.conservativeResize(rows+incRows, cols+incCols);
549  refMat1.conservativeResize(rows+incRows, cols+incCols);
550  if (incRows > 0) refMat1.bottomRows(incRows).setZero();
551  if (incCols > 0) refMat1.rightCols(incCols).setZero();
552 
553  VERIFY_IS_APPROX(m1, refMat1);
554 
555  // Insert new values
556  if (incRows > 0)
557  m1.insert(m1.rows()-1, 0) = refMat1(refMat1.rows()-1, 0) = 1;
558  if (incCols > 0)
559  m1.insert(0, m1.cols()-1) = refMat1(0, refMat1.cols()-1) = 1;
560 
561  VERIFY_IS_APPROX(m1, refMat1);
562 
563 
564  }
565  }
566 
567  // test Identity matrix
568  {
569  DenseMatrix refMat1 = DenseMatrix::Identity(rows, rows);
570  SparseMatrixType m1(rows, rows);
571  m1.setIdentity();
572  VERIFY_IS_APPROX(m1, refMat1);
573  for(int k=0; k<rows*rows/4; ++k)
574  {
575  Index i = internal::random<Index>(0,rows-1);
576  Index j = internal::random<Index>(0,rows-1);
577  Scalar v = internal::random<Scalar>();
578  m1.coeffRef(i,j) = v;
579  refMat1.coeffRef(i,j) = v;
580  VERIFY_IS_APPROX(m1, refMat1);
581  if(internal::random<Index>(0,10)<2)
582  m1.makeCompressed();
583  }
584  m1.setIdentity();
585  refMat1.setIdentity();
586  VERIFY_IS_APPROX(m1, refMat1);
587  }
588 
589  // test array/vector of InnerIterator
590  {
591  typedef typename SparseMatrixType::InnerIterator IteratorType;
592 
593  DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
594  SparseMatrixType m2(rows, cols);
595  initSparse<Scalar>(density, refMat2, m2);
596  IteratorType static_array[2];
597  static_array[0] = IteratorType(m2,0);
598  static_array[1] = IteratorType(m2,m2.outerSize()-1);
599  VERIFY( static_array[0] || m2.innerVector(static_array[0].outer()).nonZeros() == 0 );
600  VERIFY( static_array[1] || m2.innerVector(static_array[1].outer()).nonZeros() == 0 );
601  if(static_array[0] && static_array[1])
602  {
603  ++(static_array[1]);
604  static_array[1] = IteratorType(m2,0);
605  VERIFY( static_array[1] );
606  VERIFY( static_array[1].index() == static_array[0].index() );
607  VERIFY( static_array[1].outer() == static_array[0].outer() );
608  VERIFY( static_array[1].value() == static_array[0].value() );
609  }
610 
611  std::vector<IteratorType> iters(2);
612  iters[0] = IteratorType(m2,0);
613  iters[1] = IteratorType(m2,m2.outerSize()-1);
614  }
615 }
616 
617 
618 template<typename SparseMatrixType>
619 void big_sparse_triplet(Index rows, Index cols, double density) {
620  typedef typename SparseMatrixType::StorageIndex StorageIndex;
621  typedef typename SparseMatrixType::Scalar Scalar;
622  typedef Triplet<Scalar,Index> TripletType;
623  std::vector<TripletType> triplets;
624  double nelements = density * rows*cols;
625  VERIFY(nelements>=0 && nelements < NumTraits<StorageIndex>::highest());
626  Index ntriplets = Index(nelements);
627  triplets.reserve(ntriplets);
628  Scalar sum = Scalar(0);
629  for(Index i=0;i<ntriplets;++i)
630  {
631  Index r = internal::random<Index>(0,rows-1);
632  Index c = internal::random<Index>(0,cols-1);
633  // use positive values to prevent numerical cancellation errors in sum
634  Scalar v = numext::abs(internal::random<Scalar>());
635  triplets.push_back(TripletType(r,c,v));
636  sum += v;
637  }
638  SparseMatrixType m(rows,cols);
639  m.setFromTriplets(triplets.begin(), triplets.end());
640  VERIFY(m.nonZeros() <= ntriplets);
641  VERIFY_IS_APPROX(sum, m.sum());
642 }
643 
644 
645 void test_sparse_basic()
646 {
647  for(int i = 0; i < g_repeat; i++) {
648  int r = Eigen::internal::random<int>(1,200), c = Eigen::internal::random<int>(1,200);
649  if(Eigen::internal::random<int>(0,4) == 0) {
650  r = c; // check square matrices in 25% of tries
651  }
652  EIGEN_UNUSED_VARIABLE(r+c);
653  CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(1, 1)) ));
654  CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(8, 8)) ));
655  CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, ColMajor>(r, c)) ));
656  CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, RowMajor>(r, c)) ));
657  CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(r, c)) ));
658  CALL_SUBTEST_5(( sparse_basic(SparseMatrix<double,ColMajor,long int>(r, c)) ));
659  CALL_SUBTEST_5(( sparse_basic(SparseMatrix<double,RowMajor,long int>(r, c)) ));
660 
661  r = Eigen::internal::random<int>(1,100);
662  c = Eigen::internal::random<int>(1,100);
663  if(Eigen::internal::random<int>(0,4) == 0) {
664  r = c; // check square matrices in 25% of tries
665  }
666 
667  CALL_SUBTEST_6(( sparse_basic(SparseMatrix<double,ColMajor,short int>(short(r), short(c))) ));
668  CALL_SUBTEST_6(( sparse_basic(SparseMatrix<double,RowMajor,short int>(short(r), short(c))) ));
669  }
670 
671  // Regression test for bug 900: (manually insert higher values here, if you have enough RAM):
672  CALL_SUBTEST_3((big_sparse_triplet<SparseMatrix<float, RowMajor, int> >(10000, 10000, 0.125)));
673  CALL_SUBTEST_4((big_sparse_triplet<SparseMatrix<double, ColMajor, long int> >(10000, 10000, 0.125)));
674 
675  // Regression test for bug 1105
676 #ifdef EIGEN_TEST_PART_7
677  {
678  int n = Eigen::internal::random<int>(200,600);
679  SparseMatrix<std::complex<double>,0, long> mat(n, n);
680  std::complex<double> val;
681 
682  for(int i=0; i<n; ++i)
683  {
684  mat.coeffRef(i, i%(n/10)) = val;
685  VERIFY(mat.data().allocatedSize()<20*n);
686  }
687  }
688 #endif
689 }
test_docs.d
d
Definition: test_docs.py:28
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
DenseMatrix
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Definition: BenchSparseUtil.h:23
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:33
VERIFY_RAISES_ASSERT
#define VERIFY_RAISES_ASSERT(a)
Definition: main.h:285
max
#define max(a, b)
Definition: datatypes.h:20
Eigen::ColMajor
Definition: Constants.h:320
y
Scalar * y
Definition: level1_cplx_impl.h:102
b
Scalar * b
Definition: benchVecAdd.cpp:17
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
int
return int(ret)+1
Eigen::SparseMatrix
A versatible sparse matrix representation.
Definition: SparseMatrix.h:96
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
Eigen::Transpose
Expression of the transpose of a matrix.
Definition: Transpose.h:52
m2
MatrixType m2(n_dims)
v
ArrayXcf v
Definition: Cwise_arg.cpp:1
n
int n
Definition: BiCGSTAB_simple.cpp:1
diagonal
void diagonal(const MatrixType &m)
Definition: diagonal.cpp:12
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
big_sparse_triplet
void big_sparse_triplet(Index rows, Index cols, double density)
Definition: sparse_basic.cpp:619
DenseVector
Matrix< Scalar, Dynamic, 1 > DenseVector
Definition: BenchSparseUtil.h:24
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
dot
Scalar EIGEN_BLAS_FUNC() dot(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
Definition: level1_real_impl.h:27
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:331
Eigen::SparseMatrixBase::real
EIGEN_DEVICE_FUNC RealReturnType real() const
Definition: SparseMatrixBase.h:87
m1
Matrix3d m1
Definition: IOFormat.cpp:2
test_callbacks.value
value
Definition: test_callbacks.py:116
Eigen::g_repeat
static int g_repeat
Definition: main.h:144
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:211
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
mpfr::sum
const mpreal sum(const mpreal tab[], const unsigned long int n, int &status, mp_rnd_t mode=mpreal::get_default_rnd())
Definition: mpreal.h:2381
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
VERIFY_IS_MUCH_SMALLER_THAN
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition: main.h:335
Eigen::Triplet
A small structure to hold a non zero as a triplet (i,j,value).
Definition: SparseUtil.h:154
vec1
RowVectorXd vec1(3)
gtsam::Vector2
Eigen::Vector2d Vector2
Definition: Vector.h:42
VERIFY
#define VERIFY(a)
Definition: main.h:325
triangularView< Lower >
A triangularView< Lower >().adjoint().solveInPlace(B)
test_eigen.ref
ref
Definition: test_eigen.py:9
Eigen::RowMajor
Definition: Constants.h:322
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
sparse_basic
void sparse_basic(const SparseMatrixType &ref)
Definition: sparse_basic.cpp:17
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:38
g_realloc_count
static long g_realloc_count
Definition: sparse_basic.cpp:12
test_sparse_basic
void test_sparse_basic()
Definition: sparse_basic.cpp:645
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
abs
#define abs(x)
Definition: datatypes.h:17
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
EIGEN_UNUSED_VARIABLE
#define EIGEN_UNUSED_VARIABLE(var)
Definition: Macros.h:618

gtsam
Author(s):
autogenerated on Sat May 8 2021 02:44:17