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 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