Matlab et produit matriciel C++: même entrée, différentes sorties

Question

J'utilise Matlab 2010 dans mon ubuntu 12.04 sur un ordinateur x86, et g ++ 4.6.3. voilà comment je fais la production et les entrées:

#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; 
} 

et c'est la définition de 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; 

La chose étrange est que je fais les mêmes matrices avec les mêmes valeurs à l'intérieur Matlab et je obtenir des résultats différents. Voici le code 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 

alors que aaa est la sortie de l'implémentation C++. Je sais que l'erreur est si petite mais je veux savoir ce qui cause le problème! Je veux déboguer beaucoup de code et j'ai besoin de zéro différences pour un bon débogage.

Answer 1

Je vois que vous attendez une précision de 25 chiffres du code C++. Ceci est très improbable en utilisant le type double. Vous pouvez obtenir une meilleure précision en utilisant long double mais peut-être pas aussi musc que 25 chiffres.

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

Answer 2

La raison de la valeur différente (mais il pourrait être insignifiante) est que les algorithmes utilisés par Matlab et vous êtes pas même.

Votre algorithme est simple O(N^3) multiplication matricielle. Il existe des algorithmes spéciaux pour les matrices de petite taille à calculer efficacement, ainsi que des algorithmes complexes qui ont un meilleur comportement asymptotique que O(N^3).

Si vous êtes intéressé, voir:

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