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