sparse_basic.cpp revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba 
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
14template<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
494void 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