Перенос из столбца Major в Row Major

Question

Я пытаюсь преобразовать алгоритм, написанный в Fortran, и использует основной порядок столбцов для C, используя строковое упорядочение строк. Алгоритм использует вызовы gemv blas.

Я изменил звонки для строки основного макета, как в интерфейсе cblas:

• флага переключения транспонированной
• подкачки М и 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 

Как это объяснимо и то, что может быть сделано, чтобы сделать алгоритмы работы равна ?

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. Что касается -0x1.8e67415bce3aep-1 и -0x1.8e67415bce3bp-1, разница также ниже 1.0e-15.

-0.778131525475250773737911913486 
    -0.778131525475250995782516838517 
    0.000000000000000222044604925031 

Чтобы получить еще лучшее соглашение, вероятно, необходима четверная (или более высокая) точность.