我的矩阵乘法不像BLAS那样精确

Question

我用C++编写了一个矩阵乘法例程，并在特征解算器中使用了它。然后我意识到我的实现有精度问题，即使我已经使用了双精度浮点。我的实现结果似乎不够精确。我用BLAS测试了它，发现它的结果与BLAS的不匹配。下面的简单代码说明了这个问题，输出如下。

0: 0
1: 0
2: -1.42109e-14
test: test.cpp:59: int main(): Assertion `C1[i] == C2[i]' failed.

我可以简单地调用BLAS例程进行矩阵乘法，但我想知道为什么会发生这种情况。我不知道这是否是与架构相关的问题。我使用Intel Xeon E5-4620。

顺便说一句，我试着用"long double“定义sum，仍然得到了不同的结果。

谢谢你，爸爸

#include <stdlib.h>
#include <stdio.h>
#include <stddef.h>
#include <cblas.h>
#include <assert.h>

#include <vector>

void init_vec(std::vector<double> &vec)
{
    for (size_t i = 0; i < vec.size(); i++)
        vec[i] = ((double) random()) / ((double) random());
}

void mv(const std::vector<double> &A, const std::vector<double> &B,
        std::vector<double> &C)
{
    size_t num_rows = A.size() / B.size();
    size_t num_cols = B.size();
    for (size_t i = 0; i < num_rows; i++) {
        double sum = 0;
        for (size_t j = 0; j < num_cols; j++)
            sum += A[j * num_rows + i] * B[j];
        C[i] = sum;
    }
}

int main()
{
    std::vector<double> A(100 * 10);
    std::vector<double> B(10);
    std::vector<double> C1(100);
    std::vector<double> C2(100);
    init_vec(A);
    init_vec(B);

    mv(A, B, C1);
    cblas_dgemv(CblasColMajor, CblasNoTrans, 100, 10, 1, A.data(), 100, B.data(), 1, 0, C2.data(), 1);
    for (size_t i = 0; i < C1.size(); i++) {
        printf("%ld: %g\n", i, C1[i] - C2[i]);
        assert(C1[i] == C2[i]);
    }
}

Answer 1

我意识到现在的问题所在。BLAS确实使用了80位浮点进行乘法运算。下面代码的结果与BLAS匹配。

void mv(const std::vector<double> &A, const std::vector<double> &B,
       std::vector<double> &C)
{
    size_t num_rows = A.size() / B.size();
    size_t num_cols = B.size();
    for (size_t i = 0; i < num_rows; i++) {
        double sum = 0;
        for (size_t j = 0; j < num_cols; j++) {
            long double tmp1 = A[j * num_rows + i];
            long double tmp2 = B[j];
            sum += tmp2 * tmp1;
        }
        C[i] = sum;
    }
}

80位浮点似乎只用于乘法运算。求和仍然使用64位。我希望这能在未来帮助到一些人。