min_quad_with_fixed.h 5.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159
  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. #ifndef IGL_MIN_QUAD_WITH_FIXED_H
  9. #define IGL_MIN_QUAD_WITH_FIXED_H
  10. #include "igl_inline.h"
  11. #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
  12. #include <Eigen/Core>
  13. #include <Eigen/Dense>
  14. #include <Eigen/Sparse>
  15. // Bug in unsupported/Eigen/SparseExtra needs iostream first
  16. #include <iostream>
  17. #include <unsupported/Eigen/SparseExtra>
  18. namespace igl
  19. {
  20. template <typename T>
  21. struct min_quad_with_fixed_data;
  22. // Known Bugs: rows of Aeq **must** be linearly independent. Should be using
  23. // QR decomposition otherwise:
  24. // http://www.okstate.edu/sas/v8/sashtml/ormp/chap5/sect32.htm
  25. //
  26. // MIN_QUAD_WITH_FIXED Minimize quadratic energy Z'*A*Z + Z'*B + C with
  27. // constraints that Z(known) = Y, optionally also subject to the constraints
  28. // Aeq*Z = Beq
  29. //
  30. // Templates:
  31. // T should be a eigen matrix primitive type like int or double
  32. // Inputs:
  33. // A n by n matrix of quadratic coefficients
  34. // known list of indices to known rows in Z
  35. // Y list of fixed values corresponding to known rows in Z
  36. // Aeq m by n list of linear equality constraint coefficients
  37. // pd flag specifying whether A(unknown,unknown) is positive definite
  38. // Outputs:
  39. // data factorization struct with all necessary information to solve
  40. // using min_quad_with_fixed_solve
  41. // Returns true on success, false on error
  42. //
  43. // Benchmark: For a harmonic solve on a mesh with 325K facets, matlab 2.2
  44. // secs, igl/min_quad_with_fixed.h 7.1 secs
  45. //
  46. template <typename T, typename Derivedknown>
  47. IGL_INLINE bool min_quad_with_fixed_precompute(
  48. const Eigen::SparseMatrix<T>& A,
  49. const Eigen::PlainObjectBase<Derivedknown> & known,
  50. const Eigen::SparseMatrix<T>& Aeq,
  51. const bool pd,
  52. min_quad_with_fixed_data<T> & data
  53. );
  54. // Solves a system previously factored using min_quad_with_fixed_precompute
  55. //
  56. // Template:
  57. // T type of sparse matrix (e.g. double)
  58. // DerivedY type of Y (e.g. derived from VectorXd or MatrixXd)
  59. // DerivedZ type of Z (e.g. derived from VectorXd or MatrixXd)
  60. // Inputs:
  61. // data factorization struct with all necessary precomputation to solve
  62. // B n by 1 column of linear coefficients
  63. // Y b by 1 list of constant fixed values
  64. // Beq m by 1 list of linear equality constraint constant values
  65. // Outputs:
  66. // Z n by cols solution
  67. // sol #unknowns+#lagrange by cols solution to linear system
  68. // Returns true on success, false on error
  69. template <
  70. typename T,
  71. typename DerivedB,
  72. typename DerivedY,
  73. typename DerivedBeq,
  74. typename DerivedZ,
  75. typename Derivedsol>
  76. IGL_INLINE bool min_quad_with_fixed_solve(
  77. const min_quad_with_fixed_data<T> & data,
  78. const Eigen::PlainObjectBase<DerivedB> & B,
  79. const Eigen::PlainObjectBase<DerivedY> & Y,
  80. const Eigen::PlainObjectBase<DerivedBeq> & Beq,
  81. Eigen::PlainObjectBase<DerivedZ> & Z,
  82. Eigen::PlainObjectBase<Derivedsol> & sol);
  83. // Wrapper without sol
  84. template <
  85. typename T,
  86. typename DerivedB,
  87. typename DerivedY,
  88. typename DerivedBeq,
  89. typename DerivedZ>
  90. IGL_INLINE bool min_quad_with_fixed_solve(
  91. const min_quad_with_fixed_data<T> & data,
  92. const Eigen::PlainObjectBase<DerivedB> & B,
  93. const Eigen::PlainObjectBase<DerivedY> & Y,
  94. const Eigen::PlainObjectBase<DerivedBeq> & Beq,
  95. Eigen::PlainObjectBase<DerivedZ> & Z);
  96. }
  97. template <typename T>
  98. struct igl::min_quad_with_fixed_data
  99. {
  100. // Size of original system: number of unknowns + number of knowns
  101. int n;
  102. // Whether A(unknown,unknown) is positive definite
  103. bool Auu_pd;
  104. // Whether A(unknown,unknown) is symmetric
  105. bool Auu_sym;
  106. // Indices of known variables
  107. Eigen::VectorXi known;
  108. // Indices of unknown variables
  109. Eigen::VectorXi unknown;
  110. // Indices of lagrange variables
  111. Eigen::VectorXi lagrange;
  112. // Indices of unknown variable followed by Indices of lagrange variables
  113. Eigen::VectorXi unknown_lagrange;
  114. // Matrix multiplied against Y when constructing right hand side
  115. Eigen::SparseMatrix<T> preY;
  116. enum SolverType
  117. {
  118. LLT = 0,
  119. LDLT = 1,
  120. LU = 2,
  121. QR_LLT = 3,
  122. NUM_SOLVER_TYPES = 4
  123. } solver_type;
  124. // Solvers
  125. Eigen::SimplicialLLT <Eigen::SparseMatrix<T > > llt;
  126. Eigen::SimplicialLDLT<Eigen::SparseMatrix<T > > ldlt;
  127. Eigen::SparseLU<Eigen::SparseMatrix<T, Eigen::ColMajor>, Eigen::COLAMDOrdering<int> > lu;
  128. // QR factorization
  129. // Are rows of Aeq linearly independent?
  130. bool Aeq_li;
  131. // Columns of Aeq corresponding to unknowns
  132. int neq;
  133. Eigen::SparseQR<Eigen::SparseMatrix<T>, Eigen::COLAMDOrdering<int> > AeqTQR;
  134. Eigen::SparseMatrix<T> Aeqk;
  135. Eigen::SparseMatrix<T> Aequ;
  136. Eigen::SparseMatrix<T> Auu;
  137. Eigen::SparseMatrix<T> AeqTQ1;
  138. Eigen::SparseMatrix<T> AeqTQ1T;
  139. Eigen::SparseMatrix<T> AeqTQ2;
  140. Eigen::SparseMatrix<T> AeqTQ2T;
  141. Eigen::SparseMatrix<T> AeqTR1;
  142. Eigen::SparseMatrix<T> AeqTR1T;
  143. Eigen::SparseMatrix<T> AeqTE;
  144. Eigen::SparseMatrix<T> AeqTET;
  145. // Debug
  146. Eigen::SparseMatrix<T> NA;
  147. Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> NB;
  148. };
  149. #ifdef IGL_HEADER_ONLY
  150. # include "min_quad_with_fixed.cpp"
  151. #endif
  152. #endif