Home | History | Annotate | Download | only in test 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2008-2015 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      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 T>
     13 typename Eigen::internal::enable_if<(T::Flags&RowMajorBit)==RowMajorBit, typename T::RowXpr>::type
     14 innervec(T& A, Index i)
     15 {
     16   return A.row(i);
     17 }
     18 
     19 template<typename T>
     20 typename Eigen::internal::enable_if<(T::Flags&RowMajorBit)==0, typename T::ColXpr>::type
     21 innervec(T& A, Index i)
     22 {
     23   return A.col(i);
     24 }
     25 
     26 template<typename SparseMatrixType> void sparse_block(const SparseMatrixType& ref)
     27 {
     28   const Index rows = ref.rows();
     29   const Index cols = ref.cols();
     30   const Index inner = ref.innerSize();
     31   const Index outer = ref.outerSize();
     32 
     33   typedef typename SparseMatrixType::Scalar Scalar;
     34   typedef typename SparseMatrixType::StorageIndex StorageIndex;
     35 
     36   double density = (std::max)(8./(rows*cols), 0.01);
     37   typedef Matrix<Scalar,Dynamic,Dynamic,SparseMatrixType::IsRowMajor?RowMajor:ColMajor> DenseMatrix;
     38   typedef Matrix<Scalar,Dynamic,1> DenseVector;
     39   typedef Matrix<Scalar,1,Dynamic> RowDenseVector;
     40   typedef SparseVector<Scalar> SparseVectorType;
     41 
     42   Scalar s1 = internal::random<Scalar>();
     43   {
     44     SparseMatrixType m(rows, cols);
     45     DenseMatrix refMat = DenseMatrix::Zero(rows, cols);
     46     initSparse<Scalar>(density, refMat, m);
     47 
     48     VERIFY_IS_APPROX(m, refMat);
     49 
     50     // test InnerIterators and Block expressions
     51     for (int t=0; t<10; ++t)
     52     {
     53       Index j = internal::random<Index>(0,cols-2);
     54       Index i = internal::random<Index>(0,rows-2);
     55       Index w = internal::random<Index>(1,cols-j);
     56       Index h = internal::random<Index>(1,rows-i);
     57 
     58       VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w));
     59       for(Index c=0; c<w; c++)
     60       {
     61         VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c));
     62         for(Index r=0; r<h; r++)
     63         {
     64           VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r));
     65           VERIFY_IS_APPROX(m.block(i,j,h,w).coeff(r,c), refMat.block(i,j,h,w).coeff(r,c));
     66         }
     67       }
     68       for(Index r=0; r<h; r++)
     69       {
     70         VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r));
     71         for(Index c=0; c<w; c++)
     72         {
     73           VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c));
     74           VERIFY_IS_APPROX(m.block(i,j,h,w).coeff(r,c), refMat.block(i,j,h,w).coeff(r,c));
     75         }
     76       }
     77 
     78       VERIFY_IS_APPROX(m.middleCols(j,w), refMat.middleCols(j,w));
     79       VERIFY_IS_APPROX(m.middleRows(i,h), refMat.middleRows(i,h));
     80       for(Index r=0; r<h; r++)
     81       {
     82         VERIFY_IS_APPROX(m.middleCols(j,w).row(r), refMat.middleCols(j,w).row(r));
     83         VERIFY_IS_APPROX(m.middleRows(i,h).row(r), refMat.middleRows(i,h).row(r));
     84         for(Index c=0; c<w; c++)
     85         {
     86           VERIFY_IS_APPROX(m.col(c).coeff(r), refMat.col(c).coeff(r));
     87           VERIFY_IS_APPROX(m.row(r).coeff(c), refMat.row(r).coeff(c));
     88 
     89           VERIFY_IS_APPROX(m.middleCols(j,w).coeff(r,c), refMat.middleCols(j,w).coeff(r,c));
     90           VERIFY_IS_APPROX(m.middleRows(i,h).coeff(r,c), refMat.middleRows(i,h).coeff(r,c));
     91           if(m.middleCols(j,w).coeff(r,c) != Scalar(0))
     92           {
     93             VERIFY_IS_APPROX(m.middleCols(j,w).coeffRef(r,c), refMat.middleCols(j,w).coeff(r,c));
     94           }
     95           if(m.middleRows(i,h).coeff(r,c) != Scalar(0))
     96           {
     97             VERIFY_IS_APPROX(m.middleRows(i,h).coeff(r,c), refMat.middleRows(i,h).coeff(r,c));
     98           }
     99         }
    100       }
    101       for(Index c=0; c<w; c++)
    102       {
    103         VERIFY_IS_APPROX(m.middleCols(j,w).col(c), refMat.middleCols(j,w).col(c));
    104         VERIFY_IS_APPROX(m.middleRows(i,h).col(c), refMat.middleRows(i,h).col(c));
    105       }
    106     }
    107 
    108     for(Index c=0; c<cols; c++)
    109     {
    110       VERIFY_IS_APPROX(m.col(c) + m.col(c), (m + m).col(c));
    111       VERIFY_IS_APPROX(m.col(c) + m.col(c), refMat.col(c) + refMat.col(c));
    112     }
    113 
    114     for(Index r=0; r<rows; r++)
    115     {
    116       VERIFY_IS_APPROX(m.row(r) + m.row(r), (m + m).row(r));
    117       VERIFY_IS_APPROX(m.row(r) + m.row(r), refMat.row(r) + refMat.row(r));
    118     }
    119   }
    120 
    121   // test innerVector()
    122   {
    123     DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
    124     SparseMatrixType m2(rows, cols);
    125     initSparse<Scalar>(density, refMat2, m2);
    126     Index j0 = internal::random<Index>(0,outer-1);
    127     Index j1 = internal::random<Index>(0,outer-1);
    128     Index r0 = internal::random<Index>(0,rows-1);
    129     Index c0 = internal::random<Index>(0,cols-1);
    130 
    131     VERIFY_IS_APPROX(m2.innerVector(j0), innervec(refMat2,j0));
    132     VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), innervec(refMat2,j0)+innervec(refMat2,j1));
    133 
    134     m2.innerVector(j0) *= Scalar(2);
    135     innervec(refMat2,j0) *= Scalar(2);
    136     VERIFY_IS_APPROX(m2, refMat2);
    137 
    138     m2.row(r0) *= Scalar(3);
    139     refMat2.row(r0) *= Scalar(3);
    140     VERIFY_IS_APPROX(m2, refMat2);
    141 
    142     m2.col(c0) *= Scalar(4);
    143     refMat2.col(c0) *= Scalar(4);
    144     VERIFY_IS_APPROX(m2, refMat2);
    145 
    146     m2.row(r0) /= Scalar(3);
    147     refMat2.row(r0) /= Scalar(3);
    148     VERIFY_IS_APPROX(m2, refMat2);
    149 
    150     m2.col(c0) /= Scalar(4);
    151     refMat2.col(c0) /= Scalar(4);
    152     VERIFY_IS_APPROX(m2, refMat2);
    153 
    154     SparseVectorType v1;
    155     VERIFY_IS_APPROX(v1 = m2.col(c0) * 4, refMat2.col(c0)*4);
    156     VERIFY_IS_APPROX(v1 = m2.row(r0) * 4, refMat2.row(r0).transpose()*4);
    157 
    158     SparseMatrixType m3(rows,cols);
    159     m3.reserve(VectorXi::Constant(outer,int(inner/2)));
    160     for(Index j=0; j<outer; ++j)
    161       for(Index k=0; k<(std::min)(j,inner); ++k)
    162         m3.insertByOuterInner(j,k) = internal::convert_index<StorageIndex>(k+1);
    163     for(Index j=0; j<(std::min)(outer, inner); ++j)
    164     {
    165       VERIFY(j==numext::real(m3.innerVector(j).nonZeros()));
    166       if(j>0)
    167         VERIFY(j==numext::real(m3.innerVector(j).lastCoeff()));
    168     }
    169     m3.makeCompressed();
    170     for(Index j=0; j<(std::min)(outer, inner); ++j)
    171     {
    172       VERIFY(j==numext::real(m3.innerVector(j).nonZeros()));
    173       if(j>0)
    174         VERIFY(j==numext::real(m3.innerVector(j).lastCoeff()));
    175     }
    176 
    177     VERIFY(m3.innerVector(j0).nonZeros() == m3.transpose().innerVector(j0).nonZeros());
    178 
    179 //     m2.innerVector(j0) = 2*m2.innerVector(j1);
    180 //     refMat2.col(j0) = 2*refMat2.col(j1);
    181 //     VERIFY_IS_APPROX(m2, refMat2);
    182   }
    183 
    184   // test innerVectors()
    185   {
    186     DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
    187     SparseMatrixType m2(rows, cols);
    188     initSparse<Scalar>(density, refMat2, m2);
    189     if(internal::random<float>(0,1)>0.5f) m2.makeCompressed();
    190     Index j0 = internal::random<Index>(0,outer-2);
    191     Index j1 = internal::random<Index>(0,outer-2);
    192     Index n0 = internal::random<Index>(1,outer-(std::max)(j0,j1));
    193     if(SparseMatrixType::IsRowMajor)
    194       VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(j0,0,n0,cols));
    195     else
    196       VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0));
    197     if(SparseMatrixType::IsRowMajor)
    198       VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0),
    199                        refMat2.middleRows(j0,n0)+refMat2.middleRows(j1,n0));
    200     else
    201       VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0),
    202                       refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0));
    203 
    204     VERIFY_IS_APPROX(m2, refMat2);
    205 
    206     VERIFY(m2.innerVectors(j0,n0).nonZeros() == m2.transpose().innerVectors(j0,n0).nonZeros());
    207 
    208     m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0);
    209     if(SparseMatrixType::IsRowMajor)
    210       refMat2.middleRows(j0,n0) = (refMat2.middleRows(j0,n0) + refMat2.middleRows(j1,n0)).eval();
    211     else
    212       refMat2.middleCols(j0,n0) = (refMat2.middleCols(j0,n0) + refMat2.middleCols(j1,n0)).eval();
    213 
    214     VERIFY_IS_APPROX(m2, refMat2);
    215   }
    216 
    217   // test generic blocks
    218   {
    219     DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
    220     SparseMatrixType m2(rows, cols);
    221     initSparse<Scalar>(density, refMat2, m2);
    222     Index j0 = internal::random<Index>(0,outer-2);
    223     Index j1 = internal::random<Index>(0,outer-2);
    224     Index n0 = internal::random<Index>(1,outer-(std::max)(j0,j1));
    225     if(SparseMatrixType::IsRowMajor)
    226       VERIFY_IS_APPROX(m2.block(j0,0,n0,cols), refMat2.block(j0,0,n0,cols));
    227     else
    228       VERIFY_IS_APPROX(m2.block(0,j0,rows,n0), refMat2.block(0,j0,rows,n0));
    229 
    230     if(SparseMatrixType::IsRowMajor)
    231       VERIFY_IS_APPROX(m2.block(j0,0,n0,cols)+m2.block(j1,0,n0,cols),
    232                       refMat2.block(j0,0,n0,cols)+refMat2.block(j1,0,n0,cols));
    233     else
    234       VERIFY_IS_APPROX(m2.block(0,j0,rows,n0)+m2.block(0,j1,rows,n0),
    235                       refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0));
    236 
    237     Index i = internal::random<Index>(0,m2.outerSize()-1);
    238     if(SparseMatrixType::IsRowMajor) {
    239       m2.innerVector(i) = m2.innerVector(i) * s1;
    240       refMat2.row(i) = refMat2.row(i) * s1;
    241       VERIFY_IS_APPROX(m2,refMat2);
    242     } else {
    243       m2.innerVector(i) = m2.innerVector(i) * s1;
    244       refMat2.col(i) = refMat2.col(i) * s1;
    245       VERIFY_IS_APPROX(m2,refMat2);
    246     }
    247 
    248     Index r0 = internal::random<Index>(0,rows-2);
    249     Index c0 = internal::random<Index>(0,cols-2);
    250     Index r1 = internal::random<Index>(1,rows-r0);
    251     Index c1 = internal::random<Index>(1,cols-c0);
    252 
    253     VERIFY_IS_APPROX(DenseVector(m2.col(c0)), refMat2.col(c0));
    254     VERIFY_IS_APPROX(m2.col(c0), refMat2.col(c0));
    255 
    256     VERIFY_IS_APPROX(RowDenseVector(m2.row(r0)), refMat2.row(r0));
    257     VERIFY_IS_APPROX(m2.row(r0), refMat2.row(r0));
    258 
    259     VERIFY_IS_APPROX(m2.block(r0,c0,r1,c1), refMat2.block(r0,c0,r1,c1));
    260     VERIFY_IS_APPROX((2*m2).block(r0,c0,r1,c1), (2*refMat2).block(r0,c0,r1,c1));
    261 
    262     if(m2.nonZeros()>0)
    263     {
    264       VERIFY_IS_APPROX(m2, refMat2);
    265       SparseMatrixType m3(rows, cols);
    266       DenseMatrix refMat3(rows, cols); refMat3.setZero();
    267       Index n = internal::random<Index>(1,10);
    268       for(Index k=0; k<n; ++k)
    269       {
    270         Index o1 = internal::random<Index>(0,outer-1);
    271         Index o2 = internal::random<Index>(0,outer-1);
    272         if(SparseMatrixType::IsRowMajor)
    273         {
    274           m3.innerVector(o1) = m2.row(o2);
    275           refMat3.row(o1) = refMat2.row(o2);
    276         }
    277         else
    278         {
    279           m3.innerVector(o1) = m2.col(o2);
    280           refMat3.col(o1) = refMat2.col(o2);
    281         }
    282         if(internal::random<bool>())
    283           m3.makeCompressed();
    284       }
    285       if(m3.nonZeros()>0)
    286       VERIFY_IS_APPROX(m3, refMat3);
    287     }
    288   }
    289 }
    290 
    291 void test_sparse_block()
    292 {
    293   for(int i = 0; i < g_repeat; i++) {
    294     int r = Eigen::internal::random<int>(1,200), c = Eigen::internal::random<int>(1,200);
    295     if(Eigen::internal::random<int>(0,4) == 0) {
    296       r = c; // check square matrices in 25% of tries
    297     }
    298     EIGEN_UNUSED_VARIABLE(r+c);
    299     CALL_SUBTEST_1(( sparse_block(SparseMatrix<double>(1, 1)) ));
    300     CALL_SUBTEST_1(( sparse_block(SparseMatrix<double>(8, 8)) ));
    301     CALL_SUBTEST_1(( sparse_block(SparseMatrix<double>(r, c)) ));
    302     CALL_SUBTEST_2(( sparse_block(SparseMatrix<std::complex<double>, ColMajor>(r, c)) ));
    303     CALL_SUBTEST_2(( sparse_block(SparseMatrix<std::complex<double>, RowMajor>(r, c)) ));
    304 
    305     CALL_SUBTEST_3(( sparse_block(SparseMatrix<double,ColMajor,long int>(r, c)) ));
    306     CALL_SUBTEST_3(( sparse_block(SparseMatrix<double,RowMajor,long int>(r, c)) ));
    307 
    308     r = Eigen::internal::random<int>(1,100);
    309     c = Eigen::internal::random<int>(1,100);
    310     if(Eigen::internal::random<int>(0,4) == 0) {
    311       r = c; // check square matrices in 25% of tries
    312     }
    313 
    314     CALL_SUBTEST_4(( sparse_block(SparseMatrix<double,ColMajor,short int>(short(r), short(c))) ));
    315     CALL_SUBTEST_4(( sparse_block(SparseMatrix<double,RowMajor,short int>(short(r), short(c))) ));
    316   }
    317 }
    318