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