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 //
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 Index> void
57 initSparse(double density,
58            Matrix<Scalar,Dynamic,Dynamic,Opt1>& refMat,
59            SparseMatrix<Scalar,Opt2,Index>& sparseMat,
60            int flags = 0,
61            std::vector<Matrix<Index,2,1> >* zeroCoords = 0,
62            std::vector<Matrix<Index,2,1> >* nonzeroCoords = 0)
63 {
64   enum { IsRowMajor = SparseMatrix<Scalar,Opt2,Index>::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       int 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         v = internal::random<Scalar>()*Scalar(3.);
81         v = v*v + Scalar(5.);
82       }
83       if ((flags & MakeLowerTriangular) && aj>ai)
84         v = Scalar(0);
85       else if ((flags & MakeUpperTriangular) && aj<ai)
86         v = Scalar(0);
87 
88       if ((flags&ForceRealDiag) && (i==j))
89         v = numext::real(v);
90 
91       if (v!=Scalar(0))
92       {
93         //sparseMat.insertBackByOuterInner(j,i) = v;
94         sparseMat.insertByOuterInner(j,i) = v;
95         if (nonzeroCoords)
96           nonzeroCoords->push_back(Matrix<Index,2,1> (ai,aj));
97       }
98       else if (zeroCoords)
99       {
100         zeroCoords->push_back(Matrix<Index,2,1> (ai,aj));
101       }
102       refMat(ai,aj) = v;
103     }
104   }
105   //sparseMat.finalize();
106 }
107 
108 template<typename Scalar,int Opt1,int Opt2,typename Index> void
109 initSparse(double density,
110            Matrix<Scalar,Dynamic,Dynamic, Opt1>& refMat,
111            DynamicSparseMatrix<Scalar, Opt2, Index>& sparseMat,
112            int flags = 0,
113            std::vector<Matrix<Index,2,1> >* zeroCoords = 0,
114            std::vector<Matrix<Index,2,1> >* nonzeroCoords = 0)
115 {
116   enum { IsRowMajor = DynamicSparseMatrix<Scalar,Opt2,Index>::IsRowMajor };
117   sparseMat.setZero();
118   sparseMat.reserve(int(refMat.rows()*refMat.cols()*density));
119   for(int j=0; j<sparseMat.outerSize(); j++)
120   {
121     sparseMat.startVec(j); // not needed for DynamicSparseMatrix
122     for(int i=0; i<sparseMat.innerSize(); i++)
123     {
124       int ai(i), aj(j);
125       if(IsRowMajor)
126         std::swap(ai,aj);
127       Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
128       if ((flags&ForceNonZeroDiag) && (i==j))
129       {
130         v = internal::random<Scalar>()*Scalar(3.);
131         v = v*v + Scalar(5.);
132       }
133       if ((flags & MakeLowerTriangular) && aj>ai)
134         v = Scalar(0);
135       else if ((flags & MakeUpperTriangular) && aj<ai)
136         v = Scalar(0);
137 
138       if ((flags&ForceRealDiag) && (i==j))
139         v = numext::real(v);
140 
141       if (v!=Scalar(0))
142       {
143         sparseMat.insertBackByOuterInner(j,i) = v;
144         if (nonzeroCoords)
145           nonzeroCoords->push_back(Matrix<Index,2,1> (ai,aj));
146       }
147       else if (zeroCoords)
148       {
149         zeroCoords->push_back(Matrix<Index,2,1> (ai,aj));
150       }
151       refMat(ai,aj) = v;
152     }
153   }
154   sparseMat.finalize();
155 }
156 
157 template<typename Scalar,int Options,typename Index> void
158 initSparse(double density,
159            Matrix<Scalar,Dynamic,1>& refVec,
160            SparseVector<Scalar,Options,Index>& sparseVec,
161            std::vector<int>* zeroCoords = 0,
162            std::vector<int>* nonzeroCoords = 0)
163 {
164   sparseVec.reserve(int(refVec.size()*density));
165   sparseVec.setZero();
166   for(Index i=0; i<refVec.size(); i++)
167   {
168     Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
169     if (v!=Scalar(0))
170     {
171       sparseVec.insertBack(i) = v;
172       if (nonzeroCoords)
173         nonzeroCoords->push_back(i);
174     }
175     else if (zeroCoords)
176         zeroCoords->push_back(i);
177     refVec[i] = v;
178   }
179 }
180 
181 template<typename Scalar,int Options,typename Index> void
182 initSparse(double density,
183            Matrix<Scalar,1,Dynamic>& refVec,
184            SparseVector<Scalar,Options,Index>& sparseVec,
185            std::vector<int>* zeroCoords = 0,
186            std::vector<int>* nonzeroCoords = 0)
187 {
188   sparseVec.reserve(int(refVec.size()*density));
189   sparseVec.setZero();
190   for(int i=0; i<refVec.size(); i++)
191   {
192     Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
193     if (v!=Scalar(0))
194     {
195       sparseVec.insertBack(i) = v;
196       if (nonzeroCoords)
197         nonzeroCoords->push_back(i);
198     }
199     else if (zeroCoords)
200         zeroCoords->push_back(i);
201     refVec[i] = v;
202   }
203 }
204 
205 
206 #include <unsupported/Eigen/SparseExtra>
207 #endif // EIGEN_TESTSPARSE_H
208