1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra. Eigen itself is part of the KDE project.
3 //
4 // Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #include "sparse.h"
11 
12 template<typename SetterType,typename DenseType, typename Scalar, int Options>
test_random_setter(SparseMatrix<Scalar,Options> & sm,const DenseType & ref,const std::vector<Vector2i> & nonzeroCoords)13 bool test_random_setter(SparseMatrix<Scalar,Options>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords)
14 {
15   typedef SparseMatrix<Scalar,Options> SparseType;
16   {
17     sm.setZero();
18     SetterType w(sm);
19     std::vector<Vector2i> remaining = nonzeroCoords;
20     while(!remaining.empty())
21     {
22       int i = ei_random<int>(0,remaining.size()-1);
23       w(remaining[i].x(),remaining[i].y()) = ref.coeff(remaining[i].x(),remaining[i].y());
24       remaining[i] = remaining.back();
25       remaining.pop_back();
26     }
27   }
28   return sm.isApprox(ref);
29 }
30 
31 template<typename SetterType,typename DenseType, typename T>
test_random_setter(DynamicSparseMatrix<T> & sm,const DenseType & ref,const std::vector<Vector2i> & nonzeroCoords)32 bool test_random_setter(DynamicSparseMatrix<T>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords)
33 {
34   sm.setZero();
35   std::vector<Vector2i> remaining = nonzeroCoords;
36   while(!remaining.empty())
37   {
38     int i = ei_random<int>(0,remaining.size()-1);
39     sm.coeffRef(remaining[i].x(),remaining[i].y()) = ref.coeff(remaining[i].x(),remaining[i].y());
40     remaining[i] = remaining.back();
41     remaining.pop_back();
42   }
43   return sm.isApprox(ref);
44 }
45 
sparse_basic(const SparseMatrixType & ref)46 template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref)
47 {
48   const int rows = ref.rows();
49   const int cols = ref.cols();
50   typedef typename SparseMatrixType::Scalar Scalar;
51   enum { Flags = SparseMatrixType::Flags };
52 
53   double density = std::max(8./(rows*cols), 0.01);
54   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
55   typedef Matrix<Scalar,Dynamic,1> DenseVector;
56   Scalar eps = 1e-6;
57 
58   SparseMatrixType m(rows, cols);
59   DenseMatrix refMat = DenseMatrix::Zero(rows, cols);
60   DenseVector vec1 = DenseVector::Random(rows);
61   Scalar s1 = ei_random<Scalar>();
62 
63   std::vector<Vector2i> zeroCoords;
64   std::vector<Vector2i> nonzeroCoords;
65   initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords);
66 
67   if (zeroCoords.size()==0 || nonzeroCoords.size()==0)
68     return;
69 
70   // test coeff and coeffRef
71   for (int i=0; i<(int)zeroCoords.size(); ++i)
72   {
73     VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps );
74     if(ei_is_same_type<SparseMatrixType,SparseMatrix<Scalar,Flags> >::ret)
75       VERIFY_RAISES_ASSERT( m.coeffRef(zeroCoords[0].x(),zeroCoords[0].y()) = 5 );
76   }
77   VERIFY_IS_APPROX(m, refMat);
78 
79   m.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
80   refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
81 
82   VERIFY_IS_APPROX(m, refMat);
83   /*
84   // test InnerIterators and Block expressions
85   for (int t=0; t<10; ++t)
86   {
87     int j = ei_random<int>(0,cols-1);
88     int i = ei_random<int>(0,rows-1);
89     int w = ei_random<int>(1,cols-j-1);
90     int h = ei_random<int>(1,rows-i-1);
91 
92 //     VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w));
93     for(int c=0; c<w; c++)
94     {
95       VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c));
96       for(int r=0; r<h; r++)
97       {
98 //         VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r));
99       }
100     }
101 //     for(int r=0; r<h; r++)
102 //     {
103 //       VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r));
104 //       for(int c=0; c<w; c++)
105 //       {
106 //         VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c));
107 //       }
108 //     }
109   }
110 
111   for(int c=0; c<cols; c++)
112   {
113     VERIFY_IS_APPROX(m.col(c) + m.col(c), (m + m).col(c));
114     VERIFY_IS_APPROX(m.col(c) + m.col(c), refMat.col(c) + refMat.col(c));
115   }
116 
117   for(int r=0; r<rows; r++)
118   {
119     VERIFY_IS_APPROX(m.row(r) + m.row(r), (m + m).row(r));
120     VERIFY_IS_APPROX(m.row(r) + m.row(r), refMat.row(r) + refMat.row(r));
121   }
122   */
123 
124   // test SparseSetters
125   // coherent setter
126   // TODO extend the MatrixSetter
127 //   {
128 //     m.setZero();
129 //     VERIFY_IS_NOT_APPROX(m, refMat);
130 //     SparseSetter<SparseMatrixType, FullyCoherentAccessPattern> w(m);
131 //     for (int i=0; i<nonzeroCoords.size(); ++i)
132 //     {
133 //       w->coeffRef(nonzeroCoords[i].x(),nonzeroCoords[i].y()) = refMat.coeff(nonzeroCoords[i].x(),nonzeroCoords[i].y());
134 //     }
135 //   }
136 //   VERIFY_IS_APPROX(m, refMat);
137 
138   // random setter
139 //   {
140 //     m.setZero();
141 //     VERIFY_IS_NOT_APPROX(m, refMat);
142 //     SparseSetter<SparseMatrixType, RandomAccessPattern> w(m);
143 //     std::vector<Vector2i> remaining = nonzeroCoords;
144 //     while(!remaining.empty())
145 //     {
146 //       int i = ei_random<int>(0,remaining.size()-1);
147 //       w->coeffRef(remaining[i].x(),remaining[i].y()) = refMat.coeff(remaining[i].x(),remaining[i].y());
148 //       remaining[i] = remaining.back();
149 //       remaining.pop_back();
150 //     }
151 //   }
152 //   VERIFY_IS_APPROX(m, refMat);
153 
154     VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, StdMapTraits> >(m,refMat,nonzeroCoords) ));
155     #ifdef EIGEN_UNORDERED_MAP_SUPPORT
156     VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, StdUnorderedMapTraits> >(m,refMat,nonzeroCoords) ));
157     #endif
158     #ifdef _DENSE_HASH_MAP_H_
159     VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleDenseHashMapTraits> >(m,refMat,nonzeroCoords) ));
160     #endif
161     #ifdef _SPARSE_HASH_MAP_H_
162     VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleSparseHashMapTraits> >(m,refMat,nonzeroCoords) ));
163     #endif
164 
165     // test fillrand
166     {
167       DenseMatrix m1(rows,cols);
168       m1.setZero();
169       SparseMatrixType m2(rows,cols);
170       m2.startFill();
171       for (int j=0; j<cols; ++j)
172       {
173         for (int k=0; k<rows/2; ++k)
174         {
175           int i = ei_random<int>(0,rows-1);
176           if (m1.coeff(i,j)==Scalar(0))
177             m2.fillrand(i,j) = m1(i,j) = ei_random<Scalar>();
178         }
179       }
180       m2.endFill();
181       VERIFY_IS_APPROX(m2,m1);
182     }
183 
184   // test RandomSetter
185   /*{
186     SparseMatrixType m1(rows,cols), m2(rows,cols);
187     DenseMatrix refM1 = DenseMatrix::Zero(rows, rows);
188     initSparse<Scalar>(density, refM1, m1);
189     {
190       Eigen::RandomSetter<SparseMatrixType > setter(m2);
191       for (int j=0; j<m1.outerSize(); ++j)
192         for (typename SparseMatrixType::InnerIterator i(m1,j); i; ++i)
193           setter(i.index(), j) = i.value();
194     }
195     VERIFY_IS_APPROX(m1, m2);
196   }*/
197 //   std::cerr << m.transpose() << "\n\n"  << refMat.transpose() << "\n\n";
198 //   VERIFY_IS_APPROX(m, refMat);
199 
200   // test basic computations
201   {
202     DenseMatrix refM1 = DenseMatrix::Zero(rows, rows);
203     DenseMatrix refM2 = DenseMatrix::Zero(rows, rows);
204     DenseMatrix refM3 = DenseMatrix::Zero(rows, rows);
205     DenseMatrix refM4 = DenseMatrix::Zero(rows, rows);
206     SparseMatrixType m1(rows, rows);
207     SparseMatrixType m2(rows, rows);
208     SparseMatrixType m3(rows, rows);
209     SparseMatrixType m4(rows, rows);
210     initSparse<Scalar>(density, refM1, m1);
211     initSparse<Scalar>(density, refM2, m2);
212     initSparse<Scalar>(density, refM3, m3);
213     initSparse<Scalar>(density, refM4, m4);
214 
215     VERIFY_IS_APPROX(m1+m2, refM1+refM2);
216     VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3);
217     VERIFY_IS_APPROX(m3.cwise()*(m1+m2), refM3.cwise()*(refM1+refM2));
218     VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2);
219 
220     VERIFY_IS_APPROX(m1*=s1, refM1*=s1);
221     VERIFY_IS_APPROX(m1/=s1, refM1/=s1);
222 
223     VERIFY_IS_APPROX(m1+=m2, refM1+=refM2);
224     VERIFY_IS_APPROX(m1-=m2, refM1-=refM2);
225 
226     VERIFY_IS_APPROX(m1.col(0).eigen2_dot(refM2.row(0)), refM1.col(0).eigen2_dot(refM2.row(0)));
227 
228     refM4.setRandom();
229     // sparse cwise* dense
230     VERIFY_IS_APPROX(m3.cwise()*refM4, refM3.cwise()*refM4);
231 //     VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4);
232   }
233 
234   // test innerVector()
235   {
236     DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
237     SparseMatrixType m2(rows, rows);
238     initSparse<Scalar>(density, refMat2, m2);
239     int j0 = ei_random(0,rows-1);
240     int j1 = ei_random(0,rows-1);
241     VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0));
242     VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1));
243     //m2.innerVector(j0) = 2*m2.innerVector(j1);
244     //refMat2.col(j0) = 2*refMat2.col(j1);
245     //VERIFY_IS_APPROX(m2, refMat2);
246   }
247 
248   // test innerVectors()
249   {
250     DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
251     SparseMatrixType m2(rows, rows);
252     initSparse<Scalar>(density, refMat2, m2);
253     int j0 = ei_random(0,rows-2);
254     int j1 = ei_random(0,rows-2);
255     int n0 = ei_random<int>(1,rows-std::max(j0,j1));
256     VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0));
257     VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0),
258                      refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0));
259     //m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0);
260     //refMat2.block(0,j0,rows,n0) = refMat2.block(0,j0,rows,n0) + refMat2.block(0,j1,rows,n0);
261   }
262 
263   // test transpose
264   {
265     DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
266     SparseMatrixType m2(rows, rows);
267     initSparse<Scalar>(density, refMat2, m2);
268     VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval());
269     VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose());
270   }
271 
272   // test prune
273   {
274     SparseMatrixType m2(rows, rows);
275     DenseMatrix refM2(rows, rows);
276     refM2.setZero();
277     int countFalseNonZero = 0;
278     int countTrueNonZero = 0;
279     m2.startFill();
280     for (int j=0; j<m2.outerSize(); ++j)
281       for (int i=0; i<m2.innerSize(); ++i)
282       {
283         float x = ei_random<float>(0,1);
284         if (x<0.1)
285         {
286           // do nothing
287         }
288         else if (x<0.5)
289         {
290           countFalseNonZero++;
291           m2.fill(i,j) = Scalar(0);
292         }
293         else
294         {
295           countTrueNonZero++;
296           m2.fill(i,j) = refM2(i,j) = Scalar(1);
297         }
298       }
299     m2.endFill();
300     VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros());
301     VERIFY_IS_APPROX(m2, refM2);
302     m2.prune(1);
303     VERIFY(countTrueNonZero==m2.nonZeros());
304     VERIFY_IS_APPROX(m2, refM2);
305   }
306 }
307 
test_eigen2_sparse_basic()308 void test_eigen2_sparse_basic()
309 {
310   for(int i = 0; i < g_repeat; i++) {
311     CALL_SUBTEST_1( sparse_basic(SparseMatrix<double>(8, 8)) );
312     CALL_SUBTEST_2( sparse_basic(SparseMatrix<std::complex<double> >(16, 16)) );
313     CALL_SUBTEST_1( sparse_basic(SparseMatrix<double>(33, 33)) );
314 
315     CALL_SUBTEST_3( sparse_basic(DynamicSparseMatrix<double>(8, 8)) );
316   }
317 }
318