sparse.h 6.1 KB

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