I'm trying to use Eigen::LLT to solve H*x=b, this is code:

Code: Select all

typedef Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Matrixf;
typedef Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, 1>> Vectorf;

Matrixf diff = A_adj - A_adj.transpose();
std::cout << "diff norm: " << diff.norm() << std::endl;
Eigen::LLT<Matrixf, Eigen::Upper> lltf_u(A_adj);
Eigen::LLT<Matrixf, Eigen::Lower> lltf_l(A_adj);

float* result_u_ptr = new float[num_rows];
VectorfRef result_u(result_u_ptr, num_rows);
result_u = lltf_u.solve(b_adj);

float* result_l_ptr = new float[num_rows];
VectorfRef result_l(result_l_ptr, num_rows);
result_l = lltf_l.solve(b_adj);


I'm very sure that diff_norm is zero, so A_adj is symmetric.
The lltf_l with _UpLo==Lower gives the result I expected, but lltf_u with _UpLo==Upper gives different result that is obvious wrong.
Note that A_adj and b_adj's data type is float, and A_adj may has large condition number.
what makes me confused is that:
1. why _UpLo==Lower and _UpLo==Upper have different results;
2. why _UpLo==Lower gives resonable result, but _UpLo==Upper not;

sorry, i don't know how to upload txt attachment, if you need the data of A_adj and b_adj, i can give you through email.