使用犰狳进行基准测试时，SuperLu 和 LaPack 的比较失败答案

【问题标题】：Comparison of SuperLu and LaPack fails when benchmarking with armadillo使用犰狳进行基准测试时，SuperLu 和 LaPack 的比较失败
【发布时间】：2017-05-26 19:58:12
【问题描述】：

在犰狳中使用spsolve() 时，我想比较 SuperLu 的稀疏求解器和使用 LaPack 的密集求解器的速度。于是我写了这个程序：

#include "stdafx.h"
#include <iostream>
#include <Windows.h>
#include <armadillo\armadillo>

#define SIZE 2500
#define ROUNDS 2500

int main()
{
    //Time measurement stuff
    LARGE_INTEGER frequency;
    LARGE_INTEGER t1, t2, t3, t4;
    QueryPerformanceFrequency(&frequency);
    //Other stuff
    arma::cx_colvec b = arma::randu<arma::cx_colvec>(SIZE);
    arma::cx_colvec b1 = b, b2 = b;
    arma::sp_cx_mat A = arma::sp_cx_mat(SIZE, SIZE);
    A.diag(-2).fill(-1);
    A.diag(-1).fill(16);
    A.diag(0).fill(-30);
    A.diag(1).fill(16);
    A.diag(2).fill(-1);


    arma::cx_colvec c = arma::zeros<arma::cx_colvec>(SIZE), d = arma::zeros<arma::cx_colvec>(SIZE); 
    QueryPerformanceCounter(&t1);
    for (size_t i = 0; i < ROUNDS; i++)
    {
        if(arma::spsolve(c, A, b1, "superlu") == false)
        {
            std::cout << "Error in round 1\n";
            break;
        }
        b1 = c;
    }
    QueryPerformanceCounter(&t2);
    QueryPerformanceCounter(&t3);
    for (size_t i = 0; i < ROUNDS; i++)
    {
        if(arma::spsolve(d, A, b2, "lapack") == false)
        {
            std::cout << "Error in round 2\n";
            break;
        }
        b2 = d;
    }
    QueryPerformanceCounter(&t4);
    std::cout << "Superlu took " << (t2.QuadPart - t1.QuadPart)*1000.0 / frequency.QuadPart << '\n';
    std::cout << "Lapack took " << (t4.QuadPart - t3.QuadPart)*1000.0 / frequency.QuadPart << '\n';
    std::cout << "Both results are equal: " << arma::approx_equal(b1, b2, "abstol", 1e-5) << '\n';
    return 0;
}

现在对于 SIZE 和 ROUND 的小值，函数 approx_equal 返回 true，但对于较大的值，根据 approx_equal，结果 b1 和 b2 不再相等。为什么？问题可能是我的超级库无法正常工作吗？

【问题讨论】：

标签： c++ lapack armadillo

【解决方案1】：

我不会责怪 SuperLU 库。这里的“问题”似乎是矩阵A 的最小特征值随着SIZE 的值越来越大而变得越来越小。现在，for 循环重复地将inv(A) 应用于给定向量。由于您开始的向量是随机的，因此它将具有A 的特征向量的一些非零“混合”，对应于最小特征值。如果反演多次重复，该分量会显着放大，因此向量b1/b2 的各个分量会变大。

例如，对于SIZE=2000 和ROUNDS=2，我得到解的最大分量（绝对值）在10^9 附近。您的测试似乎规定了10^-5 的绝对容差。但是，对于如此大的数字，这意味着 14 个有效数字必须完全匹配，这几乎是双精度的限制。在我看来，鉴于这里比较的数字的性质，测试会更有意义，例如，approx_equal(b1, b2, "reldiff", 1E-8) 的相对误差。

另外，应该检查解决方案是否真的有意义——对于大量的ROUNDS，它迟早会溢出。例如已经使用SIZE=2000 和ROUNDS=80，我在b1/b2 向量中得到无穷大...

【讨论】：