2010-05-15 9 views
2

Supposons que vous voulez calculer la somme des carrés des différences de points:calcul mathématique d'optimisation (multiplication et addition)

$ \ sum_ {i = 1}^{N-1} (x_i - x_ { i + 1})^2

$

le code le plus simple (l'entrée est std::vector<double> xs, l'ouput sum2) est:

double sum2 = 0.; 
double prev = xs[0]; 
for (vector::const_iterator i = xs.begin() + 1; 
i != xs.end(); ++i) 
{ 
sum2 += (prev - (*i)) * (prev - (*i)); // only 1 - with compiler optimization 
prev = (*i); 
} 

J'espère que le compilateur ne l'optimisation dans le commentaire ci-dessus. Si N est la longueur de xs vous avez N-1 multiplications et 2N-3 sommes (sommes signifie + ou -).

Supposons maintenant que vous connaissez cette variable:

x 1 $^2 + x_n^2 + 2 \ sum_ {i = 2}^{N-1} x_i^2

$

et appelez-sum . Extension du carré binomial:

$ sum_i^{N-1} (x_i-x_ {i + 1})^2 = sum - 2 \ somme_ {i = 1}^{N-1} x_i x_ {i +1} $

si le code devient:

double sum2 = 0.; 
double prev = xs[0]; 
for (vector::const_iterator i = xs.begin() + 1; 
i != xs.end(); ++i) 
{ 
sum2 += (*i) * prev; 
prev = (*i); 
} 
sum2 = -sum2 * 2. + sum; 

ici, j'ai N multiplications et les additions N-1. Dans mon cas N est environ 100.

Eh bien, la compilation avec g++ -O2 Je n'ai pas accéléré (j'essaie d'appeler la fonction inline 2M fois), pourquoi?

+1

Veuillez essayer de formater les éléments LaTex correctement. Le SO le soutient. –

+0

D'une certaine manière, je soupçonne qu'une opération matricielle pourrait aider. Essayez d'utiliser GPU ou vectoriser ceci :) –

+0

@ Hamish1: êtes-vous sûr? Comment? @ Hamish2: une vectorisation de l'opération devrait aider, mais maintenant je ne pense qu'à une simplification du calcul mathématique. –

Répondre

2

Les multiplications sont beaucoup plus coûteuses que les ajouts en termes de temps d'exécution. En outre, en fonction du processeur, les ajouts et les multiplications se feront en parallèle. C'est à dire. il commencera la multiplication suivante pendant qu'il fait l'addition (voir http://en.wikipedia.org/wiki/Out-of-order_execution). Donc, la réduction du nombre d'ajouts n'aidera pas beaucoup pour la performance.

Ce que vous pouvez faire est de faciliter la tâche du compilateur pour vectoriser votre code ou le vectoriser par vous-même. Pour faciliter la vectorisation du compilateur, j'utiliserais un tableau régulier de doubles, j'utiliserais des indices et non des pointeurs.

EDIT: N = 100 peut également être un petit nombre pour voir la différence de temps d'exécution. Essayez un N plus grand.

Code sale mais montrant l'amélioration de perf. Sortie:

1e+06 
59031558 
1e+06 
18710703 

L'accélération obtenue est ~ 3x.

#include <vector> 
#include <iostream> 

using namespace std; 

unsigned long long int rdtsc(void) 
{ 
    unsigned long long int x; 
    unsigned a, d; 

    __asm__ volatile("rdtsc" : "=a" (a), "=d" (d)); 

    return ((unsigned long long)a) | (((unsigned long long)d) << 32);; 
} 



double f(std::vector<double>& xs) 
{ 
    double sum2 = 0.; 
    double prev = xs[0]; 

    vector<double>::const_iterator iend = xs.end(); 
    for (vector<double>::const_iterator i = xs.begin() + 1; 
     i != iend; ++i) 
    { 
     sum2 += (prev - (*i)) * (prev - (*i)); // only 1 - with compiler optimization 
     prev = (*i); 
    } 

    return sum2; 
} 

double f2(double *xs, int N) 
{ 
    double sum2 = 0; 

    for(int i = 0; i < N - 1; i+=1) { 
    sum2 += (xs[i+1] - xs[i])*(xs[i+1] - xs[i]); 

    } 

    return sum2; 
} 

int main(int argc, char* argv[]) 
{ 
    int N = 1000001; 
    std::vector<double> xs; 
    for(int i=0; i<N; i++) { 
    xs.push_back(i); 
    } 

    unsigned long long int a, b; 
    a = rdtsc(); 
    std::cout << f(xs) << endl; 
    b = rdtsc(); 
    cout << b - a << endl; 

    a = rdtsc(); 
    std::cout << f2(&xs[0], N) << endl; 
    b = rdtsc(); 
    cout << b - a << endl; 
} 
+0

N est fixé par le problème –

+0

Droit. Mais si vous le pouvez, il est plus facile d'avoir N plus gros pour vérifier les accélérations. – Kamchatka

+0

a) Pourquoi utilisez-vous i + = 1 au lieu de ++ i? b) l'accélération est due au fait que vous utilisez un tableau simple à la place de std :: vector ou parce que vous n'utilisez pas la variable 'prev'? –

1

L'ajout peut être libre lorsque x + = a * b. Le compilateur devrait être capable de comprendre cela dans la première version, si l'architecture le supporte.

Le calcul est probablement effectué en parallèle avec *i, ce qui pourrait être plus lent.

N'appelez pas xs.end() à chaque itération de boucle, à moins que vous ne vous attendiez à ce que la valeur de retour change. Si le compilateur ne peut pas l'optimiser, cela éclipsera le reste de la boucle.