TutorialInplaceLU.cpp 1.6 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061
  1. #include <iostream>
  2. struct init {
  3. init() { std::cout << "[" << "init" << "]" << std::endl; }
  4. };
  5. init init_obj;
  6. // [init]
  7. #include <iostream>
  8. #include <Eigen/Dense>
  9. using namespace std;
  10. using namespace Eigen;
  11. int main()
  12. {
  13. MatrixXd A(2,2);
  14. A << 2, -1, 1, 3;
  15. cout << "Here is the input matrix A before decomposition:\n" << A << endl;
  16. cout << "[init]" << endl;
  17. cout << "[declaration]" << endl;
  18. PartialPivLU<Ref<MatrixXd> > lu(A);
  19. cout << "Here is the input matrix A after decomposition:\n" << A << endl;
  20. cout << "[declaration]" << endl;
  21. cout << "[matrixLU]" << endl;
  22. cout << "Here is the matrix storing the L and U factors:\n" << lu.matrixLU() << endl;
  23. cout << "[matrixLU]" << endl;
  24. cout << "[solve]" << endl;
  25. MatrixXd A0(2,2); A0 << 2, -1, 1, 3;
  26. VectorXd b(2); b << 1, 2;
  27. VectorXd x = lu.solve(b);
  28. cout << "Residual: " << (A0 * x - b).norm() << endl;
  29. cout << "[solve]" << endl;
  30. cout << "[modifyA]" << endl;
  31. A << 3, 4, -2, 1;
  32. x = lu.solve(b);
  33. cout << "Residual: " << (A0 * x - b).norm() << endl;
  34. cout << "[modifyA]" << endl;
  35. cout << "[recompute]" << endl;
  36. A0 = A; // save A
  37. lu.compute(A);
  38. x = lu.solve(b);
  39. cout << "Residual: " << (A0 * x - b).norm() << endl;
  40. cout << "[recompute]" << endl;
  41. cout << "[recompute_bis0]" << endl;
  42. MatrixXd A1(2,2);
  43. A1 << 5,-2,3,4;
  44. lu.compute(A1);
  45. cout << "Here is the input matrix A1 after decomposition:\n" << A1 << endl;
  46. cout << "[recompute_bis0]" << endl;
  47. cout << "[recompute_bis1]" << endl;
  48. x = lu.solve(b);
  49. cout << "Residual: " << (A1 * x - b).norm() << endl;
  50. cout << "[recompute_bis1]" << endl;
  51. }