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 #include "sparse.h"
13
sparse_basic(const SparseMatrixType & ref)14 template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref)
15 {
16 typedef typename SparseMatrixType::Index Index;
17 typedef Matrix<Index,2,1> Vector2;
18
19 const Index rows = ref.rows();
20 const Index cols = ref.cols();
21 typedef typename SparseMatrixType::Scalar Scalar;
22 enum { Flags = SparseMatrixType::Flags };
23
24 double density = (std::max)(8./(rows*cols), 0.01);
25 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
26 typedef Matrix<Scalar,Dynamic,1> DenseVector;
27 Scalar eps = 1e-6;
28
29 Scalar s1 = internal::random<Scalar>();
30 {
31 SparseMatrixType m(rows, cols);
32 DenseMatrix refMat = DenseMatrix::Zero(rows, cols);
33 DenseVector vec1 = DenseVector::Random(rows);
34
35 std::vector<Vector2> zeroCoords;
36 std::vector<Vector2> nonzeroCoords;
37 initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords);
38
39 if (zeroCoords.size()==0 || nonzeroCoords.size()==0)
40 return;
41
42 // test coeff and coeffRef
43 for (int i=0; i<(int)zeroCoords.size(); ++i)
44 {
45 VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps );
46 if(internal::is_same<SparseMatrixType,SparseMatrix<Scalar,Flags> >::value)
47 VERIFY_RAISES_ASSERT( m.coeffRef(zeroCoords[0].x(),zeroCoords[0].y()) = 5 );
48 }
49 VERIFY_IS_APPROX(m, refMat);
50
51 m.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
52 refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
53
54 VERIFY_IS_APPROX(m, refMat);
55 /*
56 // test InnerIterators and Block expressions
57 for (int t=0; t<10; ++t)
58 {
59 int j = internal::random<int>(0,cols-1);
60 int i = internal::random<int>(0,rows-1);
61 int w = internal::random<int>(1,cols-j-1);
62 int h = internal::random<int>(1,rows-i-1);
63
64 // VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w));
65 for(int c=0; c<w; c++)
66 {
67 VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c));
68 for(int r=0; r<h; r++)
69 {
70 // VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r));
71 }
72 }
73 // for(int r=0; r<h; r++)
74 // {
75 // VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r));
76 // for(int c=0; c<w; c++)
77 // {
78 // VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c));
79 // }
80 // }
81 }
82
83 for(int c=0; c<cols; c++)
84 {
85 VERIFY_IS_APPROX(m.col(c) + m.col(c), (m + m).col(c));
86 VERIFY_IS_APPROX(m.col(c) + m.col(c), refMat.col(c) + refMat.col(c));
87 }
88
89 for(int r=0; r<rows; r++)
90 {
91 VERIFY_IS_APPROX(m.row(r) + m.row(r), (m + m).row(r));
92 VERIFY_IS_APPROX(m.row(r) + m.row(r), refMat.row(r) + refMat.row(r));
93 }
94 */
95
96 // test assertion
97 VERIFY_RAISES_ASSERT( m.coeffRef(-1,1) = 0 );
98 VERIFY_RAISES_ASSERT( m.coeffRef(0,m.cols()) = 0 );
99 }
100
101 // test insert (inner 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 (Index j=0; j<cols; ++j)
109 {
110 for (Index k=0; k<rows/2; ++k)
111 {
112 Index i = internal::random<Index>(0,rows-1);
113 if (m1.coeff(i,j)==Scalar(0))
114 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>();
115 }
116 }
117 m2.finalize();
118 VERIFY_IS_APPROX(m2,m1);
119 }
120
121 // test insert (fully random)
122 {
123 DenseMatrix m1(rows,cols);
124 m1.setZero();
125 SparseMatrixType m2(rows,cols);
126 if(internal::random<int>()%2)
127 m2.reserve(VectorXi::Constant(m2.outerSize(), 2));
128 for (int k=0; k<rows*cols; ++k)
129 {
130 Index i = internal::random<Index>(0,rows-1);
131 Index j = internal::random<Index>(0,cols-1);
132 if ((m1.coeff(i,j)==Scalar(0)) && (internal::random<int>()%2))
133 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>();
134 else
135 {
136 Scalar v = internal::random<Scalar>();
137 m2.coeffRef(i,j) += v;
138 m1(i,j) += v;
139 }
140 }
141 VERIFY_IS_APPROX(m2,m1);
142 }
143
144 // test insert (un-compressed)
145 for(int mode=0;mode<4;++mode)
146 {
147 DenseMatrix m1(rows,cols);
148 m1.setZero();
149 SparseMatrixType m2(rows,cols);
150 VectorXi r(VectorXi::Constant(m2.outerSize(), ((mode%2)==0) ? m2.innerSize() : std::max<int>(1,m2.innerSize()/8)));
151 m2.reserve(r);
152 for (int k=0; k<rows*cols; ++k)
153 {
154 Index i = internal::random<Index>(0,rows-1);
155 Index j = internal::random<Index>(0,cols-1);
156 if (m1.coeff(i,j)==Scalar(0))
157 m2.insert(i,j) = m1(i,j) = internal::random<Scalar>();
158 if(mode==3)
159 m2.reserve(r);
160 }
161 if(internal::random<int>()%2)
162 m2.makeCompressed();
163 VERIFY_IS_APPROX(m2,m1);
164 }
165
166 // test innerVector()
167 {
168 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
169 SparseMatrixType m2(rows, rows);
170 initSparse<Scalar>(density, refMat2, m2);
171 Index j0 = internal::random<Index>(0,rows-1);
172 Index j1 = internal::random<Index>(0,rows-1);
173 if(SparseMatrixType::IsRowMajor)
174 VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.row(j0));
175 else
176 VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0));
177
178 if(SparseMatrixType::IsRowMajor)
179 VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.row(j0)+refMat2.row(j1));
180 else
181 VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1));
182
183 SparseMatrixType m3(rows,rows);
184 m3.reserve(VectorXi::Constant(rows,rows/2));
185 for(Index j=0; j<rows; ++j)
186 for(Index k=0; k<j; ++k)
187 m3.insertByOuterInner(j,k) = k+1;
188 for(Index j=0; j<rows; ++j)
189 {
190 VERIFY(j==numext::real(m3.innerVector(j).nonZeros()));
191 if(j>0)
192 VERIFY(j==numext::real(m3.innerVector(j).lastCoeff()));
193 }
194 m3.makeCompressed();
195 for(Index j=0; j<rows; ++j)
196 {
197 VERIFY(j==numext::real(m3.innerVector(j).nonZeros()));
198 if(j>0)
199 VERIFY(j==numext::real(m3.innerVector(j).lastCoeff()));
200 }
201
202 //m2.innerVector(j0) = 2*m2.innerVector(j1);
203 //refMat2.col(j0) = 2*refMat2.col(j1);
204 //VERIFY_IS_APPROX(m2, refMat2);
205 }
206
207 // test innerVectors()
208 {
209 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
210 SparseMatrixType m2(rows, rows);
211 initSparse<Scalar>(density, refMat2, m2);
212 if(internal::random<float>(0,1)>0.5) m2.makeCompressed();
213
214 Index j0 = internal::random<Index>(0,rows-2);
215 Index j1 = internal::random<Index>(0,rows-2);
216 Index n0 = internal::random<Index>(1,rows-(std::max)(j0,j1));
217 if(SparseMatrixType::IsRowMajor)
218 VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(j0,0,n0,cols));
219 else
220 VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0));
221 if(SparseMatrixType::IsRowMajor)
222 VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0),
223 refMat2.middleRows(j0,n0)+refMat2.middleRows(j1,n0));
224 else
225 VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0),
226 refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0));
227
228 VERIFY_IS_APPROX(m2, refMat2);
229
230 m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0);
231 if(SparseMatrixType::IsRowMajor)
232 refMat2.middleRows(j0,n0) = (refMat2.middleRows(j0,n0) + refMat2.middleRows(j1,n0)).eval();
233 else
234 refMat2.middleCols(j0,n0) = (refMat2.middleCols(j0,n0) + refMat2.middleCols(j1,n0)).eval();
235
236 VERIFY_IS_APPROX(m2, refMat2);
237
238 }
239
240 // test basic computations
241 {
242 DenseMatrix refM1 = DenseMatrix::Zero(rows, rows);
243 DenseMatrix refM2 = DenseMatrix::Zero(rows, rows);
244 DenseMatrix refM3 = DenseMatrix::Zero(rows, rows);
245 DenseMatrix refM4 = DenseMatrix::Zero(rows, rows);
246 SparseMatrixType m1(rows, rows);
247 SparseMatrixType m2(rows, rows);
248 SparseMatrixType m3(rows, rows);
249 SparseMatrixType m4(rows, rows);
250 initSparse<Scalar>(density, refM1, m1);
251 initSparse<Scalar>(density, refM2, m2);
252 initSparse<Scalar>(density, refM3, m3);
253 initSparse<Scalar>(density, refM4, m4);
254
255 VERIFY_IS_APPROX(m1+m2, refM1+refM2);
256 VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3);
257 VERIFY_IS_APPROX(m3.cwiseProduct(m1+m2), refM3.cwiseProduct(refM1+refM2));
258 VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2);
259
260 VERIFY_IS_APPROX(m1*=s1, refM1*=s1);
261 VERIFY_IS_APPROX(m1/=s1, refM1/=s1);
262
263 VERIFY_IS_APPROX(m1+=m2, refM1+=refM2);
264 VERIFY_IS_APPROX(m1-=m2, refM1-=refM2);
265
266 if(SparseMatrixType::IsRowMajor)
267 VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.row(0).dot(refM2.row(0)));
268 else
269 VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.col(0).dot(refM2.row(0)));
270
271 VERIFY_IS_APPROX(m1.conjugate(), refM1.conjugate());
272 VERIFY_IS_APPROX(m1.real(), refM1.real());
273
274 refM4.setRandom();
275 // sparse cwise* dense
276 VERIFY_IS_APPROX(m3.cwiseProduct(refM4), refM3.cwiseProduct(refM4));
277 // VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4);
278
279 // test aliasing
280 VERIFY_IS_APPROX((m1 = -m1), (refM1 = -refM1));
281 VERIFY_IS_APPROX((m1 = m1.transpose()), (refM1 = refM1.transpose().eval()));
282 VERIFY_IS_APPROX((m1 = -m1.transpose()), (refM1 = -refM1.transpose().eval()));
283 VERIFY_IS_APPROX((m1 += -m1), (refM1 += -refM1));
284 }
285
286 // test transpose
287 {
288 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
289 SparseMatrixType m2(rows, rows);
290 initSparse<Scalar>(density, refMat2, m2);
291 VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval());
292 VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose());
293
294 VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint());
295 }
296
297
298
299 // test generic blocks
300 {
301 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
302 SparseMatrixType m2(rows, rows);
303 initSparse<Scalar>(density, refMat2, m2);
304 Index j0 = internal::random<Index>(0,rows-2);
305 Index j1 = internal::random<Index>(0,rows-2);
306 Index n0 = internal::random<Index>(1,rows-(std::max)(j0,j1));
307 if(SparseMatrixType::IsRowMajor)
308 VERIFY_IS_APPROX(m2.block(j0,0,n0,cols), refMat2.block(j0,0,n0,cols));
309 else
310 VERIFY_IS_APPROX(m2.block(0,j0,rows,n0), refMat2.block(0,j0,rows,n0));
311
312 if(SparseMatrixType::IsRowMajor)
313 VERIFY_IS_APPROX(m2.block(j0,0,n0,cols)+m2.block(j1,0,n0,cols),
314 refMat2.block(j0,0,n0,cols)+refMat2.block(j1,0,n0,cols));
315 else
316 VERIFY_IS_APPROX(m2.block(0,j0,rows,n0)+m2.block(0,j1,rows,n0),
317 refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0));
318
319 Index i = internal::random<Index>(0,m2.outerSize()-1);
320 if(SparseMatrixType::IsRowMajor) {
321 m2.innerVector(i) = m2.innerVector(i) * s1;
322 refMat2.row(i) = refMat2.row(i) * s1;
323 VERIFY_IS_APPROX(m2,refMat2);
324 } else {
325 m2.innerVector(i) = m2.innerVector(i) * s1;
326 refMat2.col(i) = refMat2.col(i) * s1;
327 VERIFY_IS_APPROX(m2,refMat2);
328 }
329 }
330
331 // test prune
332 {
333 SparseMatrixType m2(rows, rows);
334 DenseMatrix refM2(rows, rows);
335 refM2.setZero();
336 int countFalseNonZero = 0;
337 int countTrueNonZero = 0;
338 for (Index j=0; j<m2.outerSize(); ++j)
339 {
340 m2.startVec(j);
341 for (Index i=0; i<m2.innerSize(); ++i)
342 {
343 float x = internal::random<float>(0,1);
344 if (x<0.1)
345 {
346 // do nothing
347 }
348 else if (x<0.5)
349 {
350 countFalseNonZero++;
351 m2.insertBackByOuterInner(j,i) = Scalar(0);
352 }
353 else
354 {
355 countTrueNonZero++;
356 m2.insertBackByOuterInner(j,i) = Scalar(1);
357 if(SparseMatrixType::IsRowMajor)
358 refM2(j,i) = Scalar(1);
359 else
360 refM2(i,j) = Scalar(1);
361 }
362 }
363 }
364 m2.finalize();
365 VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros());
366 VERIFY_IS_APPROX(m2, refM2);
367 m2.prune(Scalar(1));
368 VERIFY(countTrueNonZero==m2.nonZeros());
369 VERIFY_IS_APPROX(m2, refM2);
370 }
371
372 // test setFromTriplets
373 {
374 typedef Triplet<Scalar,Index> TripletType;
375 std::vector<TripletType> triplets;
376 int ntriplets = rows*cols;
377 triplets.reserve(ntriplets);
378 DenseMatrix refMat(rows,cols);
379 refMat.setZero();
380 for(int i=0;i<ntriplets;++i)
381 {
382 Index r = internal::random<Index>(0,rows-1);
383 Index c = internal::random<Index>(0,cols-1);
384 Scalar v = internal::random<Scalar>();
385 triplets.push_back(TripletType(r,c,v));
386 refMat(r,c) += v;
387 }
388 SparseMatrixType m(rows,cols);
389 m.setFromTriplets(triplets.begin(), triplets.end());
390 VERIFY_IS_APPROX(m, refMat);
391 }
392
393 // test triangularView
394 {
395 DenseMatrix refMat2(rows, rows), refMat3(rows, rows);
396 SparseMatrixType m2(rows, rows), m3(rows, rows);
397 initSparse<Scalar>(density, refMat2, m2);
398 refMat3 = refMat2.template triangularView<Lower>();
399 m3 = m2.template triangularView<Lower>();
400 VERIFY_IS_APPROX(m3, refMat3);
401
402 refMat3 = refMat2.template triangularView<Upper>();
403 m3 = m2.template triangularView<Upper>();
404 VERIFY_IS_APPROX(m3, refMat3);
405
406 refMat3 = refMat2.template triangularView<UnitUpper>();
407 m3 = m2.template triangularView<UnitUpper>();
408 VERIFY_IS_APPROX(m3, refMat3);
409
410 refMat3 = refMat2.template triangularView<UnitLower>();
411 m3 = m2.template triangularView<UnitLower>();
412 VERIFY_IS_APPROX(m3, refMat3);
413
414 refMat3 = refMat2.template triangularView<StrictlyUpper>();
415 m3 = m2.template triangularView<StrictlyUpper>();
416 VERIFY_IS_APPROX(m3, refMat3);
417
418 refMat3 = refMat2.template triangularView<StrictlyLower>();
419 m3 = m2.template triangularView<StrictlyLower>();
420 VERIFY_IS_APPROX(m3, refMat3);
421 }
422
423 // test selfadjointView
424 if(!SparseMatrixType::IsRowMajor)
425 {
426 DenseMatrix refMat2(rows, rows), refMat3(rows, rows);
427 SparseMatrixType m2(rows, rows), m3(rows, rows);
428 initSparse<Scalar>(density, refMat2, m2);
429 refMat3 = refMat2.template selfadjointView<Lower>();
430 m3 = m2.template selfadjointView<Lower>();
431 VERIFY_IS_APPROX(m3, refMat3);
432 }
433
434 // test sparseView
435 {
436 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
437 SparseMatrixType m2(rows, rows);
438 initSparse<Scalar>(density, refMat2, m2);
439 VERIFY_IS_APPROX(m2.eval(), refMat2.sparseView().eval());
440 }
441
442 // test diagonal
443 {
444 DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
445 SparseMatrixType m2(rows, rows);
446 initSparse<Scalar>(density, refMat2, m2);
447 VERIFY_IS_APPROX(m2.diagonal(), refMat2.diagonal().eval());
448 }
449
450 // test conservative resize
451 {
452 std::vector< std::pair<Index,Index> > inc;
453 inc.push_back(std::pair<Index,Index>(-3,-2));
454 inc.push_back(std::pair<Index,Index>(0,0));
455 inc.push_back(std::pair<Index,Index>(3,2));
456 inc.push_back(std::pair<Index,Index>(3,0));
457 inc.push_back(std::pair<Index,Index>(0,3));
458
459 for(size_t i = 0; i< inc.size(); i++) {
460 Index incRows = inc[i].first;
461 Index incCols = inc[i].second;
462 SparseMatrixType m1(rows, cols);
463 DenseMatrix refMat1 = DenseMatrix::Zero(rows, cols);
464 initSparse<Scalar>(density, refMat1, m1);
465
466 m1.conservativeResize(rows+incRows, cols+incCols);
467 refMat1.conservativeResize(rows+incRows, cols+incCols);
468 if (incRows > 0) refMat1.bottomRows(incRows).setZero();
469 if (incCols > 0) refMat1.rightCols(incCols).setZero();
470
471 VERIFY_IS_APPROX(m1, refMat1);
472
473 // Insert new values
474 if (incRows > 0)
475 m1.insert(m1.rows()-1, 0) = refMat1(refMat1.rows()-1, 0) = 1;
476 if (incCols > 0)
477 m1.insert(0, m1.cols()-1) = refMat1(0, refMat1.cols()-1) = 1;
478
479 VERIFY_IS_APPROX(m1, refMat1);
480
481
482 }
483 }
484
485 // test Identity matrix
486 {
487 DenseMatrix refMat1 = DenseMatrix::Identity(rows, rows);
488 SparseMatrixType m1(rows, rows);
489 m1.setIdentity();
490 VERIFY_IS_APPROX(m1, refMat1);
491 }
492 }
493
test_sparse_basic()494 void test_sparse_basic()
495 {
496 for(int i = 0; i < g_repeat; i++) {
497 int s = Eigen::internal::random<int>(1,50);
498 EIGEN_UNUSED_VARIABLE(s);
499 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(8, 8)) ));
500 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, ColMajor>(s, s)) ));
501 CALL_SUBTEST_2(( sparse_basic(SparseMatrix<std::complex<double>, RowMajor>(s, s)) ));
502 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double>(s, s)) ));
503 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,ColMajor,long int>(s, s)) ));
504 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,RowMajor,long int>(s, s)) ));
505
506 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,ColMajor,short int>(short(s), short(s))) ));
507 CALL_SUBTEST_1(( sparse_basic(SparseMatrix<double,RowMajor,short int>(short(s), short(s))) ));
508 }
509 }
510