c - 从列专业转移到行专业

Question

我正在尝试将一个用 Fortran 编写并使用列主排序的算法转换为使用行主排序的 C。该算法使用 gemv blas 调用。

我改变了在 cblas 界面中对行主要布局的调用:

• 切换转置标志
• 交换M和N
• 改变主要维度

但该算法的行为并不相同。我得到不同的结果。我创建了一个显示行为的最小样本。

#include <stdio.h>

void dgemv_( const char * t, const int * m, const int * n, const double * alpha, const double * A, const int *lda, const double * X, const int * incx, 
    const double * beta, double * Y, const int *incy );

int main()
{
    const int M = 2, N = 2;
    const int one = 1;
    const double alpha = -1.0, beta = 1.0;
    const char trans = 'T';
    const char noTrans = 'N';

    double Yc[4] = { 0x1.42c7bd3b6266cp+4, 0x1.6c6ff393729dp+4, 0x1.acee1f3938c0bp-2, 0x1.b0cd5ba440d93p+0 };
    double Yr[4] = { 0x1.42c7bd3b6266cp+4, 0x1.acee1f3938c0bp-2, 0x1.6c6ff393729dp+4, 0x1.b0cd5ba440d93p+0 };

    double A[2] = { 0x1.11acee560242ap-2, 0x1p+0 };

    double Bc[2] = { 0x1.8p+2, 0x1.cp+2 };
    double Br[2] = { 0x1.8p+2, 0x1.cp+2 };

    dgemv_( &noTrans, &M, &N, &alpha, Yc, &M, A, &one, &beta, Bc, &one );

    printf("Result Column Major\n");
    printf("%a %a\n", Bc[0], Bc[1]);

    dgemv_( &trans, &N, &M, &alpha, Yr, &N, A, &one, &beta, Br, &one );

    printf("Result Row Major\n");
    printf("%a %a\n", Br[0], Br[1]);
}

我使用格式字符串 %a 来获取值的十六进制表示来比较它们。使用列主要版本的结果 vector 如下所示:

0x1.8402515a17beap-3 -0x1.8e67415bce3aep-1

虽然一个连续的专业看起来像这样:

0x1.8402515a17bep-3 -0x1.8e67415bce3bp-1

这如何解释以及可以做些什么来使算法工作相同？

Answer 1

如果结果用十进制表示比较

double x = 0x1.8402515a17beap-3, y = 0x1.8402515a17bep-3;
printf( "%40.30f\n", x );
printf( "%40.30f\n", y );
printf( "%40.30f\n", x - y );

他们同意最多 15 个有效数字

    0.189457545816338168709336287066
    0.189457545816337891153580130776
    0.000000000000000277555756156289

因此对于使用 double 进行 double 计算来说，差异似乎足够小。至于-0x1.8e67415bce3aep-1和-0x1.8e67415bce3bp-1，差值也在1.0e-15以下。

   -0.778131525475250773737911913486
   -0.778131525475250995782516838517
    0.000000000000000222044604925031

为了获得更好的一致性，可能需要四倍(或更高)精度。