lu_lagrange.h 2.6 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667
  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_LU_LAGRANGE_H
  9. #define IGL_LU_LAGRANGE_H
  10. #include "igl_inline.h"
  11. #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
  12. #include <Eigen/Dense>
  13. #include <Eigen/Sparse>
  14. namespace igl
  15. {
  16. // KNOWN BUGS: This does not seem to be correct for non-empty C
  17. //
  18. // LU_LAGRANGE Compute a LU decomposition for a special type of
  19. // matrix Q that is symmetric but not positive-definite:
  20. // Q = [A'*A C
  21. // C' 0];
  22. // where A'*A, or ATA, is given as a symmetric positive definite matrix and C
  23. // has full column-rank(?)
  24. //
  25. // [J] = lu_lagrange(ATA,C)
  26. //
  27. // Templates:
  28. // T should be a eigen matrix primitive type like int or double
  29. // Inputs:
  30. // ATA n by n square, symmetric, positive-definite system matrix, usually
  31. // the quadratic coefficients corresponding to the original unknowns in a
  32. // system
  33. // C n by m rectangular matrix corresponding the quadratic coefficients of
  34. // the original unknowns times the lagrange multipliers enforcing linear
  35. // equality constraints
  36. // Outputs:
  37. // L lower triangular matrix such that Q = L*U
  38. // U upper triangular matrix such that Q = L*U
  39. // Returns true on success, false on error
  40. //
  41. // Note: C should *not* have any empty columns. Typically C is the slice of
  42. // the linear constraints matrix Aeq concerning the unknown variables of a
  43. // quadratic optimization. Generally constraints may deal with unknowns as
  44. // well as knowns. Each linear constraint corresponds to a column of Aeq. As
  45. // long as each constraint concerns at least one unknown then the
  46. // corresponding column in C will have at least one non zero entry. If a
  47. // constraint concerns *no* unknowns, you should double check that this is a
  48. // valid constraint. How can you constrain known values to each other? This
  49. // is either a contradiction to the knowns' values or redundent. In either
  50. // case, it's not this functions responsiblilty to handle empty constraints
  51. // so you will get an error.
  52. //
  53. template <typename T>
  54. IGL_INLINE bool lu_lagrange(
  55. const Eigen::SparseMatrix<T> & ATA,
  56. const Eigen::SparseMatrix<T> & C,
  57. Eigen::SparseMatrix<T> & L,
  58. Eigen::SparseMatrix<T> & U);
  59. }
  60. #ifndef IGL_STATIC_LIBRARY
  61. # include "lu_lagrange.cpp"
  62. #endif
  63. #endif