c - 为什么 BLAS dsyrk 比简单的 C 实现更精确？

Question

我试图找出为什么 BLAS dsyrk 对称矩阵乘积 A'*A 比执行相同操作的 C 例程更精确。

我是这样测试的:我使用以下 Python 代码通过 mpmath 包非常精确地计算乘积:

#!/usr/bin/env python
# file: make_precise.py

import mpmath

DPS = 100

def write_matrix(A, label):
    fname = './dataset_%s.txt' % label
    with open(fname, 'w') as fid:
        for i in range(A.rows):
            for j in range(A.cols):
                pre = '' if j == 0 else ' '
                fid.write(pre + mpmath.nstr(A[i, j], DPS))
            fid.write('\n')
    print("Created %s" % fname)


def main():
    mpmath.mp.dps = DPS

    A = mpmath.randmatrix(10000, 10)
    write_matrix(A, 'A')

    AA = A.T * A
    write_matrix(AA, 'AA')


if __name__ == '__main__':
    main()

mpmath 包使用大约 100 位小数的精度计算乘积 A'*A。现在在 C 中，我将使用 dsyrk 的 BLAS 计算的乘积的精度与使用标准 C 代码计算的乘积的精度进行比较。 “朴素”的 C 代码基于 dsyrk.c 的第 332 - 350 行.我用来比较实现的代码是:

// file: minimal.c
#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include <math.h>

#define N_OBS 10000
#define N_VAR 10

// generate datasets with:
// python make_precise.py

// compile with:
// gcc minimal.c -o minimal -lcblas -lm

void dsyrk_aa(const double *A, double *AA)
{
    cblas_dsyrk(CblasRowMajor, CblasUpper, CblasTrans, N_VAR, N_OBS, 1.0,
            A, N_VAR, 0., AA, N_VAR);
}


void naive_aa(double *A, double *AA)
{
    int i, j, k;
    double temp;

    for (j=0; j<N_VAR; j++) {
        for (k=j; k<N_VAR; k++) {
            temp = 0.0;
            for (i=0; i<N_OBS; i++) {
                temp += A[i*N_VAR+k] * A[i*N_VAR+j];
            }
            AA[j*N_VAR+k] = temp;
        }
    }
}


double *read_data(const char *name, int rows, int cols)
{
    int i, j;
    double value, *M = NULL;

    char filename[1024];
    sprintf(filename, "./dataset_%s.txt", name);
    FILE *fid = fopen(filename, "r");

    M = malloc(sizeof(double)*rows*cols);

    for (i=0; i<rows; i++) {
        for (j=0; j<cols; j++) {
            fscanf(fid, "%lf", &value);
            M[i*cols+j] = value;
        }
    }

    fclose(fid);
    return M;
}


double dist_from_true(double *A, double *B, int rows, int cols)
{
    int i, j;
    double dist = 0.0;
    for (i=0; i<rows; i++)
        for (j=i; j<cols; j++)
            dist += fabs(A[i*cols+j] - B[i*cols+j]);
    return dist;
}


int main()
{
    double d1, d2;
    double *A = NULL,
           *AA = NULL;

    double *AA1 = calloc(N_VAR*N_VAR, sizeof(double));
    double *AA2 = calloc(N_VAR*N_VAR, sizeof(double));

    A = read_data("A", N_OBS, N_VAR);
    AA = read_data("AA", N_VAR, N_VAR);

    dsyrk_aa(A, AA1);
    naive_aa(A, AA2);

    d1 = dist_from_true(AA, AA1, N_VAR, N_VAR);
    d2 = dist_from_true(AA, AA2, N_VAR, N_VAR);

    free(A);
    free(AA);

    free(AA1);
    free(AA2);

    printf("Dsyrk:  \t%.16f\n", d1);
    printf("Naive:  \t%.16f\n", d2);

    return EXIT_SUCCESS;
}

请注意，我只在两个例程中计算上三角。当然，通过将数据集读入 C 语言，我们会损失一些精度，因为所有内容都将存储为 double。但是，我们将与使用 mpmath 计算的真实产品进行比较，因此我们应该能够比较两种产品的精度。我得到的结果是:

Dsyrk:      0.0000000000923137
Naive:      0.0000000003306013

因此，对于 BLAS，绝对误差比 C 实现小大约 3 倍。这对于多个数据集和多个工作站(均运行 Linux)是可重现的。我知道差异似乎可以忽略不计，但我正在处理更大的数据集，其中错误会随着时间的推移而累积。

我的问题是:这种差异从何而来，我可以做些什么来使 C 实现与 BLAS 实现一样精确？

提前感谢您的宝贵时间!

更新

我重新编译了 ATLAS 以查看在编译 dsyrk 例程时使用了哪些编译器标志。我已将范围缩小到这一行:

/usr/bin/x86_64-pc-linux-gnu-gcc-6.2.1 -o ATL_drefsyrkLN.o -c -DL2SIZE=33554432 -I/tmp/ATLAS/build/include -I/tmp/ATLAS/build/..//include -I/tmp/ATLAS/build/..//include/contrib -DAdd_ -DF77_INTEGER=int -DStringSunStyle -DATL_OS_Linux -DATL_ARCH_Corei4 -DATL_CPUMHZ=3200 -DATL_AVXMAC -DATL_AVX -DATL_SSE3 -DATL_SSE2 -DATL_SSE1 -DATL_USE64BITS -DATL_GAS_x8664 -m64 -DATL_DYLIBS -DPentiumCPS=3200.000 -DATL_FULL_LAPACK -DATL_NCPU=4 -fomit-frame-pointer -mfpmath=sse -O2 -mavx2 -mfma -fPIC -m64 -fPIC /tmp/ATLAS/build/..//src/blas/reference/level3/ATL_drefsyrkLN.c

我认为重要的标志是:

-m64 -fomit-frame-pointer -mfpmath=sse -O2 -mavx2 -mfma -fPIC

然而，当编译上面的最小示例时:

gcc -o minimal -m64 -fomit-frame-pointer -mfpmath=sse -O2 -mavx2 -mfma -fPIC minimal.c -lcblas -lm

精度结果不受影响。非常感谢任何帮助。

Answer 1

将解决方案作为答案发布以供将来引用。 @rubenvb 的提示很有用。我已经追踪到 cblas_dsyrk 例程在 ATLAS 库中采用的确切路径。这是如下:

1. cblas_dsyrk
2. ATL_syrk
3. ATL_dsprk
4. ATL_dsprk_rK
5. ATL_dprk_kmm
6. ATL_gNBmm
7. ATL_dpKBmm
8. ATL_dgpKBmm
9. ATL_dJIK0x0x56TN56x56x0_a1_bX
10. ATL_dJIK0x0x56TN1x1x4_a1_bX(完成所有工作)

现在，ATLAS 为我的 CPU 使用的“ block 大小”是 NB = 56。因此，为了便于测试，我将测试数据集缩减为它的倍数(上面 minimal.c 中的 N_OBS = 9968)。如果不是这种情况，将为其余行调用类似于列表中最后两个的一些例程。

找到最终函数后，我开始测试结果矩阵 AA 中单个元素的值。我注意到对 ATL_dprk_kmm 的调用实际上执行了三次，初始 A 矩阵的不同 block (行)。这些 block 的大小分别为 3472、3472 和 3024 行。对于这些 block 中的每一个，构建一个单独的结果矩阵AA，最终将其添加到用户的AA矩阵中。这立即引出了解决方案，在原始函数中，我们也以 block 为单位计算 AA:

void newnaive_aa(const double *A, double *AA)
{
    int i, j, k, blk_start, blk_end;
    double temp;
    double *tmpAA = calloc(N_VAR*N_VAR, sizeof(double));

    // block 1
    blk_start = 0; blk_end = 3472;
    for (j=0; j<N_VAR; j++) {
        for (k=j; k<N_VAR; k++) {
            temp = 0.0;
            for (i=blk_start; i<blk_end; i++) {
                temp += A[i*N_VAR+k] * A[i*N_VAR+j];
            }
            tmpAA[j*N_VAR+k] = temp;
        }
    }
    // copy matrix over
    for (j=0; j<N_VAR*N_VAR; j++)
        AA[j] += tmpAA[j];

    // block 2
    blk_start = 3472; blk_end = 6944;
    for (j=0; j<N_VAR; j++) {
        for (k=j; k<N_VAR; k++) {
            temp = 0.0;
            for (i=blk_start; i<blk_end; i++) {
                temp += A[i*N_VAR+k] * A[i*N_VAR+j];
            }
            tmpAA[j*N_VAR+k] = temp;
        }
    }
    // copy matrix over
    for (j=0; j<N_VAR*N_VAR; j++)
        AA[j] += tmpAA[j];

    // block 3
    blk_start = 6944; blk_end = 9968;
    for (j=0; j<N_VAR; j++) {
        for (k=j; k<N_VAR; k++) {
            temp = 0.0;
            for (i=blk_start; i<blk_end; i++) {
                temp += A[i*N_VAR+k] * A[i*N_VAR+j];
            }
            tmpAA[j*N_VAR+k] = temp;
        }
    }
    // copy matrix over
    for (j=0; j<N_VAR*N_VAR; j++)
        AA[j] += tmpAA[j];

    free(tmpAA);
}

通过这个新实现，可以获得与 dsyrk 调用完全相同的精度。

另请注意，对于包含较少行数的数据集(例如 N_OBS = 1000)，原始 naive_aa 的精度已经与 dsyrk 的精度相同结果。这样做的原因现在很明显，因为该数据集大小不需要中间结果。