is_symmetric.cpp 2.5 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879
  1. // This file is part of libigl, a simple c++ geometry processing library.
  2. //
  3. // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
  4. //
  5. // This Source Code Form is subject to the terms of the Mozilla Public License
  6. // v. 2.0. If a copy of the MPL was not distributed with this file, You can
  7. // obtain one at http://mozilla.org/MPL/2.0/.
  8. #include "is_symmetric.h"
  9. #include "find.h"
  10. template <typename T>
  11. IGL_INLINE bool igl::is_symmetric(const Eigen::SparseMatrix<T>& A)
  12. {
  13. if(A.rows() != A.cols())
  14. {
  15. return false;
  16. }
  17. assert(A.size() != 0);
  18. Eigen::SparseMatrix<T> AT = A.transpose();
  19. Eigen::SparseMatrix<T> AmAT = A-AT;
  20. //// Eigen screws up something with LLT if you try to do
  21. //SparseMatrix<T> AmAT = A-A.transpose();
  22. //// Eigen crashes at runtime if you try to do
  23. // return (A-A.transpose()).nonZeros() == 0;
  24. return AmAT.nonZeros() == 0;
  25. }
  26. template <typename DerivedA>
  27. IGL_INLINE bool igl::is_symmetric(
  28. const Eigen::PlainObjectBase<DerivedA>& A)
  29. {
  30. if(A.rows() != A.cols())
  31. {
  32. return false;
  33. }
  34. assert(A.size() != 0);
  35. const typename Eigen::PlainObjectBase<DerivedA>& AT = A.transpose();
  36. const typename Eigen::PlainObjectBase<DerivedA>& AmAT = A-AT;
  37. //// Eigen screws up something with LLT if you try to do
  38. //SparseMatrix<T> AmAT = A-A.transpose();
  39. //// Eigen crashes at runtime if you try to do
  40. // return (A-A.transpose()).nonZeros() == 0;
  41. return AmAT.nonZeros() == 0;
  42. }
  43. template <typename AType, typename epsilonT>
  44. IGL_INLINE bool igl::is_symmetric(
  45. const Eigen::SparseMatrix<AType>& A,
  46. const epsilonT epsilon)
  47. {
  48. using namespace Eigen;
  49. using namespace std;
  50. if(A.rows() != A.cols())
  51. {
  52. return false;
  53. }
  54. assert(A.size() != 0);
  55. SparseMatrix<AType> AT = A.transpose();
  56. SparseMatrix<AType> AmAT = A-AT;
  57. VectorXi AmATI,AmATJ;
  58. Matrix<AType,Dynamic,1> AmATV;
  59. find(AmAT,AmATI,AmATJ,AmATV);
  60. if(AmATI.size() == 0)
  61. {
  62. return true;
  63. }
  64. return AmATV.maxCoeff() < epsilon && AmATV.minCoeff() > -epsilon;
  65. }
  66. #ifdef IGL_STATIC_LIBRARY
  67. // Explicit template specialization
  68. // generated by autoexplicit.sh
  69. template bool igl::is_symmetric<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&);
  70. // generated by autoexplicit.sh
  71. template bool igl::is_symmetric<double>(Eigen::SparseMatrix<double, 0, int> const&);
  72. template bool igl::is_symmetric<double, double>(Eigen::SparseMatrix<double, 0, int> const&, double);
  73. template bool igl::is_symmetric<double, int>(Eigen::SparseMatrix<double, 0, int> const&, int);
  74. #endif