SparseLU和SparseQR在複雜矩陣中獲得differen結果

我試圖使用稀疏QR和LU來解決我的FEA程序中的複雜情況，並且似乎QR和BiCGSTAB方法無法獲得正確的結果。同時BiCGSTAB與IncompleteLUT無關。本徵版本3.3.1使用MinGW-x86_64的GCC 6.2SparseLU和SparseQR在複雜矩陣中獲得differen結果

最小碼

#include <vector> 
#include <complex> 
#include <iostream> 
#include <Eigen/Eigen> 

using namespace std::complex_literals; 

struct mat_cell 
{ 
    int row; 
    int col; 
    double value; 
}; 
// matrix data. 
mat_cell mat01[]={ 
    { 0, 0, 40432.2974517006}, { 0, 6, -20216.1487258503}, { 0, 12, -20216.1487258503}, 
    { 1, 1, 180.518062136147}, { 1, 7, -90.2590310680736}, { 1, 11, -9025.90310680736}, 
    { 1, 13, -90.2590310680736}, { 1, 17, 9025.90310680736}, 
    { 2, 2, 180.518062136147}, { 2, 8, -90.2590310680736}, { 2, 10, 9025.90310680736}, 
    { 2, 14, -90.2590310680736}, { 2, 16, -9025.90310680736}, 
    { 3, 3, 456735.213970955}, { 3, 9, -228367.606985477}, { 3, 15, -228367.606985477}, 
    { 4, 4, 2421773.15749991}, { 4, 8, -9025.90310680736}, { 4, 10, 594294.042611519}, 
    { 4, 14, 9025.90310680736}, { 4, 16, 594294.042611519}, 
    { 5, 5, 2421773.15749991}, { 5, 7, 9025.90310680736}, { 5, 11, 594294.042611519}, 
    { 5, 13, -9025.90310680736}, { 5, 17, 594294.042611519}, 
    { 6, 0, -20216.1487258503}, { 6, 6, 40432.2974517006}, { 6, 24, -20216.1487258503}, 
    { 7, 1, -90.2590310680736}, { 7, 5, 9025.90310680736}, { 7, 7, 180.518062136147}, 
    { 7, 25, -90.2590310680736}, { 7, 29, -9025.90310680736}, 
    { 8, 2, -90.2590310680736}, { 8, 4, -9025.90310680736}, { 8, 8, 180.518062136147}, 
    { 8, 26, -90.2590310680736}, { 8, 28, 9025.90310680736}, 
    { 9, 3, -228367.606985477}, { 9, 9, 456735.213970955}, { 9, 27, -228367.606985477}, 
    {10, 2, 9025.90310680736}, {10, 4, 594294.042611519}, {10, 10, 2421773.15749991}, 
    {10, 26, -9025.90310680736}, {10, 28, 594294.042611519}, 
    {11, 1, -9025.90310680736}, {11, 5, 594294.042611519}, {11, 11, 2421773.15749991}, 
    {11, 25, 9025.90310680736}, {11, 29, 594294.042611519}, 
    {12, 0, -20216.1487258503}, {12, 12, 20216.1487258503}, 
    {13, 1, -90.2590310680736}, {13, 5, -9025.90310680736}, {13, 13, 90.2590310680736}, 
    {13, 17, -9025.90310680736}, 
    {14, 2, -90.2590310680736}, {14, 4, 9025.90310680736}, {14, 14, 90.2590310680736}, 
    {14, 16, 9025.90310680736}, 
    {15, 3, -228367.606985477}, {15, 15, 228367.606985477}, 
    {16, 2, -9025.90310680736}, {16, 4, 594294.042611519}, {16, 14, 9025.90310680736}, 
    {16, 16, 1210886.57874995}, 
    {17, 1, 9025.90310680736}, {17, 5, 594294.042611519}, {17, 13, -9025.90310680736}, 
    {17, 17, 1210886.57874995}, 
    {18, 18, 40432.2974517006}, {18, 24, -20216.1487258503}, 
    {19, 19, 180.518062136147}, {19, 25, -90.2590310680736}, {19, 29, 9025.90310680736}, 
    {20, 20, 180.518062136147}, {20, 26, -90.2590310680736}, {20, 28, -9025.90310680736}, 
    {21, 21, 456735.213970955}, {21, 27, -228367.606985477}, 
    {22, 22, 2421773.15749991}, {22, 26, 9025.90310680736}, {22, 28, 594294.042611519}, 
    {23, 23, 2421773.15749991}, {23, 25, -9025.90310680736}, {23, 29, 594294.042611519}, 
    {24, 6, -20216.1487258503}, {24, 18, -20216.1487258503}, {24, 24, 40432.2974517006}, 
    {25, 7, -90.2590310680736}, {25, 11, 9025.90310680736}, {25, 19, -90.2590310680736}, 
    {25, 23, -9025.90310680736}, {25, 25, 180.518062136147}, 
    {26, 8, -90.2590310680736}, {26, 10, -9025.90310680736}, {26, 20, -90.2590310680736}, 
    {26, 22, 9025.90310680736}, {26, 26, 180.518062136147}, 
    {27, 9, -228367.606985477}, {27, 21, -228367.606985477}, {27, 27, 456735.213970955}, 
    {28, 8, 9025.90310680736}, {28, 10, 594294.042611519}, {28, 20, -9025.90310680736}, 
    {28, 22, 594294.042611519}, {28, 28, 2421773.15749991}, 
    {29, 7, -9025.90310680736}, {29, 11, 594294.042611519}, {29, 19, 9025.90310680736}, 
    {29, 23, 594294.042611519}, {29, 29, 2421773.15749991}}; 

int main(int argc, char *argv[]) 
{ 
    int nn{30}; 

    Eigen::MatrixXcd A_dens = Eigen::MatrixXcd::Zero(nn, nn); 
    Eigen::VectorXcd rhs = Eigen::VectorXcd::Zero(nn); 

    Eigen::SparseMatrix<std::complex<double>> A_sp(nn, nn); 
    std::vector<Eigen::Triplet<std::complex<double>>> triList; 

    double yita{0.02};// small imag. 
    for(auto const cell: mat01){ 
     A_dens(cell.row, cell.col) = cell.value*(1.+yita*1.0i); 
     triList.push_back({cell.row, cell.col, cell.value*(1.+yita*1.0i)}); 
    } 
    A_sp.setFromTriplets(triList.begin(), triList.end()); 
    triList.clear(); 
    A_sp.makeCompressed(); 

    int ix[]={12, 13, 14}; 
    double scale{1.e60};// Large than 1e38. 

    for(auto const j: ix){ 
     A_dens(j, j) *= scale; 
     A_sp.coeffRef(j, j) *= scale; 
    } 

    rhs(ix[1]) = 0.618*A_sp.coeff(ix[1], ix[1]); 

    // solve by dense LU method. 
    Eigen::VectorXcd x_lu = A_dens.lu().solve(rhs); 

    // define sparse solver. 
    Eigen::SparseLU<Eigen::SparseMatrix<std::complex<double>>, Eigen::COLAMDOrdering<int>> solver_lu; 
    Eigen::SparseQR<Eigen::SparseMatrix<std::complex<double>>, Eigen::COLAMDOrdering<int>> solver_qr; 
    Eigen::BiCGSTAB<Eigen::SparseMatrix<std::complex<double>>> solver_bi; 
    Eigen::BiCGSTAB<Eigen::SparseMatrix<std::complex<double>>, Eigen::IncompleteLUT<std::complex<double>, int>> solver_bi_2; 

    solver_lu.compute(A_sp); 
    if(solver_lu.info()!=Eigen::ComputationInfo::Success)std::cout << "SparseLU decomposition failed!\n"; 
    Eigen::VectorXcd x_sp_lu = solver_lu.solve(rhs); 
    if(solver_lu.info()!=Eigen::ComputationInfo::Success)std::cout << "SparseLU solve failed!\n"; 

    solver_qr.compute(A_sp); 
    if(solver_qr.info()!=Eigen::ComputationInfo::Success)std::cout << "SparseQR decomposition failed!\n"; 
    Eigen::VectorXcd x_sp_qr = solver_qr.solve(rhs); 
    if(solver_qr.info()!=Eigen::ComputationInfo::Success)std::cout << "SparseQR solve failed!\n"; 

    solver_bi.compute(A_sp); 
    if(solver_bi.info()!=Eigen::ComputationInfo::Success)std::cout << "SparseBi decomposition failed!\n"; 
    Eigen::VectorXcd x_sp_bi = solver_bi.solve(rhs); 
    if(solver_bi.info()!=Eigen::ComputationInfo::Success)std::cout << "SparseBi solve failed!\n"; 

    solver_bi_2.compute(A_sp); 
    if(solver_bi_2.info()!=Eigen::ComputationInfo::Success)std::cout << "SparseBi2 decomposition failed!\n"; 
    Eigen::VectorXcd x_sp_bi_2 = solver_bi_2.solve(rhs); 
    if(solver_bi_2.info()!=Eigen::ComputationInfo::Success)std::cout << "SparseBi2 solve failed!\n"; 

    std::cout << "No | Dense LU | SparseLU | SparseQR | BiCGSTAB |BiCGSTAB+ILUT|\n"; 
    std::cout << "---|---|---|---|---|---|\n"; 
    for(int i=0; i<nn; i++){ 
     std::cout << i << "|"; 
     std::cout << x_lu(i) << "|"; 
     std::cout << x_sp_lu(i) << "|"; 
     std::cout << x_sp_qr(i) << "|"; 
     std::cout << x_sp_bi(i) << "|"; 
     std::cout << x_sp_bi_2(i) << "|\n"; 
    } 
}

方法 X（1）X（5）

DenseLU（0.435087，-1.73121e- 017）（0.0008897,7.91857e-020）

SparseLU（0.435087,3.61979e-017）（0.0008897，-1.2936e-019）

SparseQR（0,0）（0,0）

的BiCGSTAB（0.187474，-8.66607e-019）（0.00139743，-2.34841e-021）

的BiCGSTAB + ILUT（0.435068,1.58791e-017）（0.000889823，-1.00545e-019）

More detailed result compare picture

來源

2016-12-15 pasuka

你能用較小的矩陣重現這個問題嗎？您正在使用的矩陣的條件編號是什麼？矩陣是可逆的（如果不是，你不會得到一個獨特的解決方案）。 –

你看到的結果有多不同？你是否真切地期望他們有多接近？ –

您好@MichaelAnderson，您可以在結果向量x（1）和x（5）中看到LU和BiCGSTAB + ILUT方法都可以，但QR和BiCGSTAB似乎不正確。 – pasuka

你的矩陣是奇異的高度。對於雙精度數字，由於奇怪的3列縮放，所以矩陣的排名僅爲3。因此，解決您的問題有無限的空間。如果你看一下相對誤差：(A*x-b).norm()/b.norm()，您可以：

No | Dense LU | SparseLU | SparseQR | BiCGSTAB |BiCGSTAB+ILUT 
--- | ---------- | ---------- | ---------- | ---------- |------------- 
res | 3.6633e-74 | 1.4915e-74 | 3.1977e-18 | 1.9095e-59 | 2.67692e-63

這意味着所有的結果都是「正確」相對於雙精度浮點數的精度。

來源

2016-12-15 08:16:46 ggael

SparseLU和SparseQR在複雜矩陣中獲得differen結果

回答

相關問題