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