我正在尝试转换一个算法,它是用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获取值的十六进制表示形式来比较它们。使用列主版本生成的向量如下所示:
0x1.8402515a17beap-3 -0x1.8e67415bce3aep-1
而对一个排成一列的少校来说则是这样的:
0x1.8402515a17bep-3 -0x1.8e67415bce3bp-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
的双精度计算来说,这种差异似乎足够小对于-0x1.8e67415bce3aep-1
和-0x1.8e67415bce3bp-1
,差异也小于1.0e-15。 -0.778131525475250773737911913486
-0.778131525475250995782516838517
0.000000000000000222044604925031
为了获得更好的一致性,可能需要四倍(或更高)的精度。
关于c - 从专栏专业转移到行专业,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/31007672/