You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

112 lines
4.6 KiB

  1. // This file is part of Eigen, a lightweight C++ template library
  2. // for linear algebra.
  3. //
  4. // Copyright (C) 2008-2010 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. #include "sparse.h"
  10. template<typename Scalar> void
  11. initSPD(double density,
  12. Matrix<Scalar,Dynamic,Dynamic>& refMat,
  13. SparseMatrix<Scalar>& sparseMat)
  14. {
  15. Matrix<Scalar,Dynamic,Dynamic> aux(refMat.rows(),refMat.cols());
  16. initSparse(density,refMat,sparseMat);
  17. refMat = refMat * refMat.adjoint();
  18. for (int k=0; k<2; ++k)
  19. {
  20. initSparse(density,aux,sparseMat,ForceNonZeroDiag);
  21. refMat += aux * aux.adjoint();
  22. }
  23. sparseMat.setZero();
  24. for (int j=0 ; j<sparseMat.cols(); ++j)
  25. for (int i=j ; i<sparseMat.rows(); ++i)
  26. if (refMat(i,j)!=Scalar(0))
  27. sparseMat.insert(i,j) = refMat(i,j);
  28. sparseMat.finalize();
  29. }
  30. template<typename Scalar> void sparse_solvers(int rows, int cols)
  31. {
  32. double density = (std::max)(8./(rows*cols), 0.01);
  33. typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
  34. typedef Matrix<Scalar,Dynamic,1> DenseVector;
  35. // Scalar eps = 1e-6;
  36. DenseVector vec1 = DenseVector::Random(rows);
  37. std::vector<Vector2i> zeroCoords;
  38. std::vector<Vector2i> nonzeroCoords;
  39. // test triangular solver
  40. {
  41. DenseVector vec2 = vec1, vec3 = vec1;
  42. SparseMatrix<Scalar> m2(rows, cols);
  43. DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
  44. // lower - dense
  45. initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
  46. VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
  47. m2.template triangularView<Lower>().solve(vec3));
  48. // upper - dense
  49. initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
  50. VERIFY_IS_APPROX(refMat2.template triangularView<Upper>().solve(vec2),
  51. m2.template triangularView<Upper>().solve(vec3));
  52. VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
  53. m2.conjugate().template triangularView<Upper>().solve(vec3));
  54. {
  55. SparseMatrix<Scalar> cm2(m2);
  56. //Index rows, Index cols, Index nnz, Index* outerIndexPtr, Index* innerIndexPtr, Scalar* valuePtr
  57. MappedSparseMatrix<Scalar> mm2(rows, cols, cm2.nonZeros(), cm2.outerIndexPtr(), cm2.innerIndexPtr(), cm2.valuePtr());
  58. VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
  59. mm2.conjugate().template triangularView<Upper>().solve(vec3));
  60. }
  61. // lower - transpose
  62. initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
  63. VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Upper>().solve(vec2),
  64. m2.transpose().template triangularView<Upper>().solve(vec3));
  65. // upper - transpose
  66. initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
  67. VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Lower>().solve(vec2),
  68. m2.transpose().template triangularView<Lower>().solve(vec3));
  69. SparseMatrix<Scalar> matB(rows, rows);
  70. DenseMatrix refMatB = DenseMatrix::Zero(rows, rows);
  71. // lower - sparse
  72. initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular);
  73. initSparse<Scalar>(density, refMatB, matB);
  74. refMat2.template triangularView<Lower>().solveInPlace(refMatB);
  75. m2.template triangularView<Lower>().solveInPlace(matB);
  76. VERIFY_IS_APPROX(matB.toDense(), refMatB);
  77. // upper - sparse
  78. initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular);
  79. initSparse<Scalar>(density, refMatB, matB);
  80. refMat2.template triangularView<Upper>().solveInPlace(refMatB);
  81. m2.template triangularView<Upper>().solveInPlace(matB);
  82. VERIFY_IS_APPROX(matB, refMatB);
  83. // test deprecated API
  84. initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
  85. VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
  86. m2.template triangularView<Lower>().solve(vec3));
  87. }
  88. }
  89. void test_sparse_solvers()
  90. {
  91. for(int i = 0; i < g_repeat; i++) {
  92. CALL_SUBTEST_1(sparse_solvers<double>(8, 8) );
  93. int s = internal::random<int>(1,300);
  94. CALL_SUBTEST_2(sparse_solvers<std::complex<double> >(s,s) );
  95. CALL_SUBTEST_1(sparse_solvers<double>(s,s) );
  96. }
  97. }