c – 使用Eigen 3线性系统求解器的错误结果
作者:互联网
我将一个MATLAB代码移植到C与Eigen 3,我决定构造我的线性求解器,而不是从矩阵对象调用它,所以我可以重用它.不幸的是,它没有产生预期的结果.经过几次测试后,我将问题追溯到似乎是线性求解器对象,如下面的相关代码所示:
MatrixXcd M(6,6), Y(6,6), dvalor(6,6);
// Initialization of the matrices...
FullPivHouseholderQR<MatrixXcd> mdivide;
auto result = mldivide.compute(Y).solve(M).eval() *
dvalor * mldivide.compute(M).solve(Y).eval();
std::cout << "result:" << std::endl
<< result << std::endl << std::endl;
std::cout << "using fullPivHouseholderQr:" << std::endl
<< Y.fullPivHouseholderQr().solve(M) * dvalor *
M.fullPivHouseholderQr().solve(Y) << std::endl << std::endl;
std::cout << "using jacobi:" << std::endl
<< Y.jacobiSvd(ComputeThinU | ComputeThinV).solve(M) * dvalor *
M.jacobiSvd(ComputeThinU | ComputeThinV).solve(Y) << std::endl << std::endl;
std::cout << "using inverse:" << std::endl
<< Y.inverse() * M * dvalor * M.inverse() * Y << std::endl << std::endl;
这是结果:
result:
(0.564196,0.00606298) (-0.15282,-0.00179293) (-0.16564,-0.000220726) (-0.16564,-0.000220726) (-0.15282,-0.00179293) (-0.179235,0.00541725)
(2.73563e+184,-8.89017e+185) (1,1.12229e-10) (-6.04297e-11,5.62259e-12) (-6.02256e-11,2.60221e-12) (-5.33448e-11,-3.05694e-11) (-6.55427e-11,8.17723e-12)
(2.11057e+184,-8.81209e+185) (-5.53615e-11,4.54731e-12) (1,1.11384e-10) (-5.99088e-11,9.29495e-12) (-5.51744e-11,1.52693e-12) (-6.21597e-11,-3.11176e-11)
(2.11057e+184,-8.81209e+185) (-0.171319,-0.00113681) (0.134093,-2.05541e-05) (0.865907,2.05542e-05) (0.171319,0.00113681) (-0.140389,-0.000216684)
(2.73563e+184,-8.89017e+185) (0.218904,0.00143522) (-0.171337,3.98397e-05) (0.171337,-3.98397e-05) (0.781096,-0.00143522) (0.179383,0.000262654)
(1.4668e+184,-8.78677e+185) (-5.52732e-11,5.99039e-12) (-5.69133e-11,-3.2223e-11) (-5.97328e-11,9.63533e-12) (-5.51767e-11,4.01574e-12) (1,1.1063e-10)
using fullPivHouseholderQr:
(1,1.1063e-10) (-5.51766e-11,4.01578e-12) (-5.97331e-11,9.63535e-12) (-5.6913e-11,-3.22229e-11) (-5.52733e-11,5.99029e-12) (-6.46967e-11,1.0973e-11)
(-6.54343e-11,6.20279e-12) (1,1.12229e-10) (-6.04297e-11,5.62259e-12) (-6.02256e-11,2.60221e-12) (-5.33448e-11,-3.05694e-11) (-6.55427e-11,8.17723e-12)
(-6.48883e-11,1.07408e-11) (-5.53615e-11,4.54731e-12) (1,1.11384e-10) (-5.99088e-11,9.29495e-12) (-5.51744e-11,1.52693e-12) (-6.21597e-11,-3.11176e-11)
(-6.21592e-11,-3.11175e-11) (-5.51743e-11,1.52697e-12) (-5.99093e-11,9.29482e-12) (1,1.11385e-10) (-5.53614e-11,4.54728e-12) (-6.48881e-11,1.07408e-11)
(-6.55427e-11,8.1773e-12) (-5.33444e-11,-3.05693e-11) (-6.02261e-11,2.60223e-12) (-6.04293e-11,5.62257e-12) (1,1.12229e-10) (-6.54344e-11,6.20272e-12)
(-6.46967e-11,1.0973e-11) (-5.52732e-11,5.99039e-12) (-5.69133e-11,-3.2223e-11) (-5.97328e-11,9.63533e-12) (-5.51767e-11,4.01574e-12) (1,1.1063e-10)
using jacobi:
(1,1.1063e-10) (-5.51761e-11,4.01562e-12) (-5.9733e-11,9.63488e-12) (-5.69134e-11,-3.22232e-11) (-5.5273e-11,5.9903e-12) (-6.46969e-11,1.09735e-11)
(-6.54346e-11,6.20257e-12) (1,1.12228e-10) (-6.04303e-11,5.62232e-12) (-6.0225e-11,2.60224e-12) (-5.33448e-11,-3.05694e-11) (-6.55414e-11,8.17752e-12)
(-6.48874e-11,1.07408e-11) (-5.53617e-11,4.54725e-12) (1,1.11384e-10) (-5.99097e-11,9.29465e-12) (-5.51742e-11,1.52688e-12) (-6.21596e-11,-3.11175e-11)
(-6.21608e-11,-3.11175e-11) (-5.5174e-11,1.52688e-12) (-5.99104e-11,9.29433e-12) (1,1.11384e-10) (-5.53612e-11,4.54728e-12) (-6.48873e-11,1.07413e-11)
(-6.55423e-11,8.17718e-12) (-5.33459e-11,-3.05695e-11) (-6.02258e-11,2.60206e-12) (-6.04294e-11,5.6224e-12) (1,1.12229e-10) (-6.5434e-11,6.20302e-12)
(-6.46966e-11,1.09729e-11) (-5.5273e-11,5.99037e-12) (-5.69134e-11,-3.22233e-11) (-5.97333e-11,9.63543e-12) (-5.51767e-11,4.01582e-12) (1,1.1063e-10)
using inverse:
(1,1.1063e-10) (-5.51768e-11,4.01574e-12) (-5.97329e-11,9.63535e-12) (-5.6914e-11,-3.2223e-11) (-5.52733e-11,5.99031e-12) (-6.46967e-11,1.0973e-11)
(-6.54338e-11,6.20274e-12) (1,1.12229e-10) (-6.04293e-11,5.6225e-12) (-6.02264e-11,2.6022e-12) (-5.33451e-11,-3.05693e-11) (-6.55425e-11,8.1773e-12)
(-6.48877e-11,1.07408e-11) (-5.53615e-11,4.54726e-12) (1,1.11384e-10) (-5.99097e-11,9.29488e-12) (-5.51743e-11,1.52697e-12) (-6.21592e-11,-3.11176e-11)
(-6.21587e-11,-3.11176e-11) (-5.51744e-11,1.52696e-12) (-5.99093e-11,9.29489e-12) (1,1.11384e-10) (-5.53616e-11,4.54727e-12) (-6.4888e-11,1.07408e-11)
(-6.55419e-11,8.17731e-12) (-5.33448e-11,-3.05693e-11) (-6.02257e-11,2.6022e-12) (-6.04301e-11,5.6225e-12) (1,1.12229e-10) (-6.54343e-11,6.20273e-12)
(-6.46964e-11,1.0973e-11) (-5.52732e-11,5.99031e-12) (-5.69131e-11,-3.2223e-11) (-5.97336e-11,9.63535e-12) (-5.51767e-11,4.01574e-12) (1,1.1063e-10)
MATLAB中算法的第一次迭代产生一个近似单位矩阵,可以通过反转矩阵并从矩阵对象中调用求解器在C中观察到.
但是使用求解器对象时结果完全错误.我还尝试了来自eigen的文档页面的其他求解器,结果相同:https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html.
我错过了一个步骤,还是从矩阵中调用求解器以不同方式初始化求解器?
Eigen版本是3.3.4,我还没有使用任何编译器优化,只是C 11:
g++ -I. -std=c++11 main.cpp
矩阵的大小是恒定的,它将在一个大循环内执行.因此,通过拉出解算器可以避免的内存分配越多越好.
谢谢.
解决方法:
小心auto,见common pitfalls,基本上在你的情况下结果是Product<>表达式引用死对象.用MatrixXcd替换auto并解决问题.
标签:c,eigen,linear-algebra,eigen3 来源: https://codeday.me/bug/20190828/1751230.html