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-2011 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 #ifndef EIGEN_TESTSPARSE_H
     11 #define EIGEN_TESTSPARSE_H
     12 
     13 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
     14 
     15 #include "main.h"
     16 
     17 #if EIGEN_GNUC_AT_LEAST(4,0) && !defined __ICC && !defined(__clang__)
     18 
     19 #ifdef min
     20 #undef min
     21 #endif
     22 
     23 #ifdef max
     24 #undef max
     25 #endif
     26 
     27 #include <tr1/unordered_map>
     28 #define EIGEN_UNORDERED_MAP_SUPPORT
     29 namespace std {
     30   using std::tr1::unordered_map;
     31 }
     32 #endif
     33 
     34 #ifdef EIGEN_GOOGLEHASH_SUPPORT
     35   #include <google/sparse_hash_map>
     36 #endif
     37 
     38 #include <Eigen/Cholesky>
     39 #include <Eigen/LU>
     40 #include <Eigen/Sparse>
     41 
     42 enum {
     43   ForceNonZeroDiag = 1,
     44   MakeLowerTriangular = 2,
     45   MakeUpperTriangular = 4,
     46   ForceRealDiag = 8
     47 };
     48 
     49 /* Initializes both a sparse and dense matrix with same random values,
     50  * and a ratio of \a density non zero entries.
     51  * \param flags is a union of ForceNonZeroDiag, MakeLowerTriangular and MakeUpperTriangular
     52  *        allowing to control the shape of the matrix.
     53  * \param zeroCoords and nonzeroCoords allows to get the coordinate lists of the non zero,
     54  *        and zero coefficients respectively.
     55  */
     56 template<typename Scalar,int Opt1,int Opt2,typename StorageIndex> void
     57 initSparse(double density,
     58            Matrix<Scalar,Dynamic,Dynamic,Opt1>& refMat,
     59            SparseMatrix<Scalar,Opt2,StorageIndex>& sparseMat,
     60            int flags = 0,
     61            std::vector<Matrix<StorageIndex,2,1> >* zeroCoords = 0,
     62            std::vector<Matrix<StorageIndex,2,1> >* nonzeroCoords = 0)
     63 {
     64   enum { IsRowMajor = SparseMatrix<Scalar,Opt2,StorageIndex>::IsRowMajor };
     65   sparseMat.setZero();
     66   //sparseMat.reserve(int(refMat.rows()*refMat.cols()*density));
     67   sparseMat.reserve(VectorXi::Constant(IsRowMajor ? refMat.rows() : refMat.cols(), int((1.5*density)*(IsRowMajor?refMat.cols():refMat.rows()))));
     68 
     69   for(Index j=0; j<sparseMat.outerSize(); j++)
     70   {
     71     //sparseMat.startVec(j);
     72     for(Index i=0; i<sparseMat.innerSize(); i++)
     73     {
     74       Index ai(i), aj(j);
     75       if(IsRowMajor)
     76         std::swap(ai,aj);
     77       Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
     78       if ((flags&ForceNonZeroDiag) && (i==j))
     79       {
     80         // FIXME: the following is too conservative
     81         v = internal::random<Scalar>()*Scalar(3.);
     82         v = v*v;
     83         if(numext::real(v)>0) v += Scalar(5);
     84         else                  v -= Scalar(5);
     85       }
     86       if ((flags & MakeLowerTriangular) && aj>ai)
     87         v = Scalar(0);
     88       else if ((flags & MakeUpperTriangular) && aj<ai)
     89         v = Scalar(0);
     90 
     91       if ((flags&ForceRealDiag) && (i==j))
     92         v = numext::real(v);
     93 
     94       if (v!=Scalar(0))
     95       {
     96         //sparseMat.insertBackByOuterInner(j,i) = v;
     97         sparseMat.insertByOuterInner(j,i) = v;
     98         if (nonzeroCoords)
     99           nonzeroCoords->push_back(Matrix<StorageIndex,2,1> (ai,aj));
    100       }
    101       else if (zeroCoords)
    102       {
    103         zeroCoords->push_back(Matrix<StorageIndex,2,1> (ai,aj));
    104       }
    105       refMat(ai,aj) = v;
    106     }
    107   }
    108   //sparseMat.finalize();
    109 }
    110 
    111 template<typename Scalar,int Opt1,int Opt2,typename Index> void
    112 initSparse(double density,
    113            Matrix<Scalar,Dynamic,Dynamic, Opt1>& refMat,
    114            DynamicSparseMatrix<Scalar, Opt2, Index>& sparseMat,
    115            int flags = 0,
    116            std::vector<Matrix<Index,2,1> >* zeroCoords = 0,
    117            std::vector<Matrix<Index,2,1> >* nonzeroCoords = 0)
    118 {
    119   enum { IsRowMajor = DynamicSparseMatrix<Scalar,Opt2,Index>::IsRowMajor };
    120   sparseMat.setZero();
    121   sparseMat.reserve(int(refMat.rows()*refMat.cols()*density));
    122   for(int j=0; j<sparseMat.outerSize(); j++)
    123   {
    124     sparseMat.startVec(j); // not needed for DynamicSparseMatrix
    125     for(int i=0; i<sparseMat.innerSize(); i++)
    126     {
    127       int ai(i), aj(j);
    128       if(IsRowMajor)
    129         std::swap(ai,aj);
    130       Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
    131       if ((flags&ForceNonZeroDiag) && (i==j))
    132       {
    133         v = internal::random<Scalar>()*Scalar(3.);
    134         v = v*v + Scalar(5.);
    135       }
    136       if ((flags & MakeLowerTriangular) && aj>ai)
    137         v = Scalar(0);
    138       else if ((flags & MakeUpperTriangular) && aj<ai)
    139         v = Scalar(0);
    140 
    141       if ((flags&ForceRealDiag) && (i==j))
    142         v = numext::real(v);
    143 
    144       if (v!=Scalar(0))
    145       {
    146         sparseMat.insertBackByOuterInner(j,i) = v;
    147         if (nonzeroCoords)
    148           nonzeroCoords->push_back(Matrix<Index,2,1> (ai,aj));
    149       }
    150       else if (zeroCoords)
    151       {
    152         zeroCoords->push_back(Matrix<Index,2,1> (ai,aj));
    153       }
    154       refMat(ai,aj) = v;
    155     }
    156   }
    157   sparseMat.finalize();
    158 }
    159 
    160 template<typename Scalar,int Options,typename Index> void
    161 initSparse(double density,
    162            Matrix<Scalar,Dynamic,1>& refVec,
    163            SparseVector<Scalar,Options,Index>& sparseVec,
    164            std::vector<int>* zeroCoords = 0,
    165            std::vector<int>* nonzeroCoords = 0)
    166 {
    167   sparseVec.reserve(int(refVec.size()*density));
    168   sparseVec.setZero();
    169   for(int i=0; i<refVec.size(); i++)
    170   {
    171     Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
    172     if (v!=Scalar(0))
    173     {
    174       sparseVec.insertBack(i) = v;
    175       if (nonzeroCoords)
    176         nonzeroCoords->push_back(i);
    177     }
    178     else if (zeroCoords)
    179         zeroCoords->push_back(i);
    180     refVec[i] = v;
    181   }
    182 }
    183 
    184 template<typename Scalar,int Options,typename Index> void
    185 initSparse(double density,
    186            Matrix<Scalar,1,Dynamic>& refVec,
    187            SparseVector<Scalar,Options,Index>& sparseVec,
    188            std::vector<int>* zeroCoords = 0,
    189            std::vector<int>* nonzeroCoords = 0)
    190 {
    191   sparseVec.reserve(int(refVec.size()*density));
    192   sparseVec.setZero();
    193   for(int i=0; i<refVec.size(); i++)
    194   {
    195     Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
    196     if (v!=Scalar(0))
    197     {
    198       sparseVec.insertBack(i) = v;
    199       if (nonzeroCoords)
    200         nonzeroCoords->push_back(i);
    201     }
    202     else if (zeroCoords)
    203         zeroCoords->push_back(i);
    204     refVec[i] = v;
    205   }
    206 }
    207 
    208 
    209 #include <unsupported/Eigen/SparseExtra>
    210 #endif // EIGEN_TESTSPARSE_H
    211