Matlab и C++ матричный продукт: тот же вход, разные выходы

Question

Я использую Matlab 2010 в моем ubuntu 12.04 на компьютере x86 и g ++ 4.6.3. это, как я делаю производство и входы:

#include <Src/Tools/Math/Matrix_nxn.h> 
    #include <iostream> 
    using namespace std; 
    int main() 
    { 
     Matrix_nxn<double,4> A1,A2,Tb,aa; 
     A1[0][0] = 0.99958087959447828; A1[0][1] = 1.7725781974830023e-18;A1[0][2] = 0.028949354900049871; A1[0][3] = 0; 
     A1[1][0] = -0.028949354900049871; A1[1][1] = 6.1204654815537932e-17;A1[1][2] = 0.99958087959447828; A1[1][3] = 0; 
     A1[2][0] = 0,   A1[2][1] = -1;   A1[2][2] = 6.1230317691118863e-17;A1[2][3] = 0.21129000000000001; 
     A1[3][0] = 0,   A1[3][1] = 0;   A1[3][2] = 0;    A1[3][3] = 1; 

     A2[0][0] = 0.90634806393366396; A2[0][1] = -0.42253187690835708;A2[0][2] = 0;A2[0][3] = 0; 
     A2[1][0] = 0.42253187690835708; A2[1][1] = 0.90634806393366396; A2[1][2] = 0;A2[1][3] = 0; 
     A2[2][0] = 0;   A2[2][1] = 0;   A2[2][2] = 1;A2[2][3] = 0; 
     A2[3][0] = 0;   A2[3][1] = 0;   A2[3][2] = 0;A2[3][3] = 1; 

     Tb[0][0] = 0.99956387949834924; Tb[0][1] = -0.00016363183229951183; Tb[0][2] = -0.029530052943282908; Tb[0][3] = 0; 
     Tb[1][0] = 0;   Tb[1][1] = 0.99998464792303143; Tb[1][2] = -0.0055411116439683869;Tb[1][3] = 0; 
     Tb[2][0] = 0.029530506297888514;Tb[2][1] = 0.0055386950515785164; Tb[2][2] = 0.99954853411673616; Tb[2][3] = 0; 
     Tb[3][0] = 0;   Tb[3][1] = 0;   Tb[3][2] = 0;    Tb[3][3] = 1; 


    aa = Tb*A1*A2; 
    cout.precision(25); 
    cout <<aa[0][0]<<' '<<aa[0][1]<<' '<<aa[0][2]<<' '<<aa[0][3]<<endl 
     <<aa[1][0]<<' '<<aa[1][1]<<' '<<aa[1][2]<<' '<<aa[1][3]<<endl 
     <<aa[2][0]<<' '<<aa[2][1]<<' '<<aa[2][2]<<' '<<aa[2][3]<<endl 
     <<aa[3][0]<<' '<<aa[3][1]<<' '<<aa[3][2]<<' '<<aa[3][3]<<endl; 
} 

и это определение operator*:

Matrix_nxn<T, N> res; 
size_t i, j, k; 
for (i = 0; i < N; ++i) 
{ 
    for (j = 0; j < N; ++j) 
    { 
    for (k = 0; k < N; ++k) 
    { 
     res[i][j] += m1[i][k] * m2[k][j]; 
    } 
    if (MVTools::isNearInf(res[i][j])) 
    { 
     if (MVTools::isNearPosInf(res[i][j])) 
     throw MVException(MVException::PosInfValue); 
     else 
     throw MVException(MVException::NegInfValue); 
    } 
    } 
} 
return res; 

Странная вещь я делаю те же матрицы с одинаковыми значениями внутри Matlab и I получить разные результаты. Вот код Matlab:

Tb = [0.99956387949834924,-0.00016363183229951183,-0.029530052943282908,0;0,0.99998464792303143,-0.0055411116439683869,0;0.029530506297888514,0.0055386950515785164,0.99954853411673616,0;0,0,0,1]; 
A1 = [0.99958087959447828,1.7725781974830023e-18,0.028949354900049871,0;-0.028949354900049871,6.1204654815537932e-17,0.99958087959447828,0;0,-1,6.1230317691118863e-17,0.21129000000000001;0,0,0,1]; 
A2 = [0.90634806393366396,-0.42253187690835708,0,0;0.42253187690835708,0.90634806393366396,0,0;0,0,1,0;0,0,0,1]; 
aa = Tb*A1*A2; 
aa - aaa 

ans = 

    1.0e-16 * 

       0 -0.555111512312578     0     0 
       0     0     0     0 
       0     0     0     0 
       0     0     0     0 

while aaa - это результат реализации C++. Я знаю, что ошибки так мало, но я хочу знать, что вызывает проблему! Я хочу отладить много кода, и мне нужна нулевая разница для хорошей отладки.

Answer 1

Я вижу, что вы ожидаете точности в 25 цифр от кода на C++. Это очень маловероятно, используя тип double. Вы можете получить более точную точность, используя long double, но, возможно, не так, как мускус, как 25 цифр.

См: What is the precision of long double in C++?

Answer 2

Причина другое значение (однако незначителен это может быть) в том, что алгоритмы, используемые MATLAB и вы не то же самое.

Ваш алгоритм простой O(N^3) матричное умножение. Существуют специальные алгоритмы для эффективного вычисления матриц малого размера, а также сложные алгоритмы, которые имеют лучшую асимптотику, чем O(N^3).

Если вы заинтересованы, см:

• Laderman's 3x3 matrix multiplication with only 23 multiplications, is it worth it?
• http://en.wikipedia.org/wiki/Strassen_algorithm