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
sparse_basic(const SparseMatrixType & ref)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 m1 = m4; refM1 = refM4;
232 }
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>
big_sparse_triplet(Index rows,Index cols,double density)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 Scalar v = internal::random<Scalar>();
634 triplets.push_back(TripletType(r,c,v));
635 sum += v;
636 }
637 SparseMatrixType m(rows,cols);
638 m.setFromTriplets(triplets.begin(), triplets.end());
639 VERIFY(m.nonZeros() <= ntriplets);
640 VERIFY_IS_APPROX(sum, m.sum());
641 }
642
643
test_sparse_basic()644 void test_sparse_basic()
645 {
646 for(int i = 0; i < g_repeat; i++) {
647 int r = Eigen::internal::random<int>(1,200), c = Eigen::internal::random<int>(1,200);
648 if(Eigen::internal::random<int>(0,4) == 0) {
649 r = c; // check square matrices in 25% of tries
650 }
651 EIGEN_UNUSED_VARIABLE(r+c);
652 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(1, 1)) ));
653 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(8, 8)) ));
654 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, ColMajor>(r, c)) ));
655 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, RowMajor>(r, c)) ));
656 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(r, c)) ));
657 CALL_SUBTEST_5(( sparse_basic(SparseMatrix<double,ColMajor,long int>(r, c)) ));
658 CALL_SUBTEST_5(( sparse_basic(SparseMatrix<double,RowMajor,long int>(r, c)) ));
659
660 r = Eigen::internal::random<int>(1,100);
661 c = Eigen::internal::random<int>(1,100);
662 if(Eigen::internal::random<int>(0,4) == 0) {
663 r = c; // check square matrices in 25% of tries
664 }
665
666 CALL_SUBTEST_6(( sparse_basic(SparseMatrix<double,ColMajor,short int>(short(r), short(c))) ));
667 CALL_SUBTEST_6(( sparse_basic(SparseMatrix<double,RowMajor,short int>(short(r), short(c))) ));
668 }
669
670 // Regression test for bug 900: (manually insert higher values here, if you have enough RAM):
671 CALL_SUBTEST_3((big_sparse_triplet<SparseMatrix<float, RowMajor, int> >(10000, 10000, 0.125)));
672 CALL_SUBTEST_4((big_sparse_triplet<SparseMatrix<double, ColMajor, long int> >(10000, 10000, 0.125)));
673
674 // Regression test for bug 1105
675 #ifdef EIGEN_TEST_PART_7
676 {
677 int n = Eigen::internal::random<int>(200,600);
678 SparseMatrix<std::complex<double>,0, long> mat(n, n);
679 std::complex<double> val;
680
681 for(int i=0; i<n; ++i)
682 {
683 mat.coeffRef(i, i%(n/10)) = val;
684 VERIFY(mat.data().allocatedSize()<20*n);
685 }
686 }
687 #endif
688 }
689