Problema di N-Body: Parallelizzazione efficiente del ciclo double per

Question

Un problema molto comune in un problema con N-body è l'uso di un doppio ciclo per calcolare le interazioni tra le particelle. Considerando un problema corpo N con n particelle, il ciclo può essere scritta trovi

for (i = 0, i < n; i++) 
    for (j = i+1, j < n; j++) 
     // calculate interaction 

mia domanda riguarda come può essere ciclo parallelizzata utilizzando fili differenti. L'obiettivo è che ogni thread "idealmente" debba calcolare lo stesso numero di interazioni.

La mia idea era di separare il ciclo esterno, il ciclo i, su intervalli diversi, diciamo a_k = a (k), dove k = 1,2, ..., p dove p è il numero di fili che voglio dividere il problema in

Così, il ciclo potrebbe essere scritta come

for (k = 1, k < p; k++) 
    for (i = a(k), i < a(k+1); i++) 
     for (j = i+1, j < n; j++) 
      // calculate interaction 

Dove il ciclo più esterno, il k-ciclo, è quello di essere parallelizzati.

Poiché il numero di interazioni del ciclo più interno, il j-ciclo, è N- (i + 1), il numero di interazioni calcolati da ciascun filo è

\ sum_ {i = a (k)}^{a (k + 1)} n - (i + 1)

Ciò significa che si vorrebbe trovare la funzione a_k discreta tale che il funzionale

f [a_k] = \ sum_ {i = a (k)}^{a (k + 1)} n - (i + 1)

con le condizioni al contorno a (1) = 0 e a (p) = n è una costante funzionale, forzando così t il numero di interazioni su ogni thread è lo stesso.

Ho provato a utilizzare "euristica" diversa (ad esempio a_k polinomiale, esponenziale, log), e finora nessuno mi ha dato una risposta soddisfacente. Una soluzione diretta di questo problema non è evidente per me.

Per piccolo p, questo problema può essere messo sul "problemi di minimizzazione del sacco" in cui sostanzialmente ogni a_k è una variabile per minimizzare la funzione

f (a_1, a_2, A_3, ...) = somma (| f [a_k] - n/p |^2)

Ma è possibile indovinare, questo non è efficiente (o addirittura converge) per valori più elevati di p.

Qualcuno ha idea di come potrebbe essere affrontato questo problema?

Answer 1

(Scusa se questo non è espresso chiaramente, ha senso nella mia testa).

Quando sommando tutti i numeri da 1 a N, si può notare che N + 1 = (N - 1) + 2 = (N - 2) + 3, ecc

Allora, che cosa se ogni thread usava un io piccolo e uno grande, tale che le somme venivano sempre sommate?

Oppure, supponiamo di voler utilizzare sempre 5 thread. Thread 1 farebbe il primo 10% e l'ultimo 10%, il thread 2 farebbe il secondo 10% e il penultimo 10%, e così via. Ogni abbinamento di una sezione "in anticipo" e "in ritardo" equivarrebbe allo stesso numero totale di interazioni.

EDIT:

Rubare un diagramma da un altro post ...

0 1 2 3 4 5 6 7 8 

0 - A B C D D C B A 
1 - B C D D C B A 
2  - C D D C B A 
3  - D D C B A 
4   - D C B A 
5   - C B A 
6    - B A 
7    - A 
8     - 

Questo mostra più chiaramente cosa intendo?

Answer 2

È possibile dividere gli oggetti in k gruppi di circa N/k corpi, e usare questo per sezionare il vostro triangolo iniziale di interazioni in k*(k + 1)/2 pezzi:

0 1 2 3 4 5 6 7 8 
         -- N=9; k=3; N/k=3 
0 - A A B B B C C C 
1 - A B B B C C C -- diagonal pieces: A, D, F 
2  - B B B C C C 
3  - D D E E E -- non-diagonal pieces: B, C, E 
4   - D E E E 
5   - E E E 
6    - F F 
7    - F 
8     - 

Questo punto di vista è complicato dal fatto che ci sono due tipi di pezzi: quelli lungo la diagonale (che sono triangoli con elementi (N/k)*(N/k - 1)/2) e quelli che non lo sono (che sono quadrati con elementi (N/k)*(N/k)). Tuttavia, poiché i pezzi diagonali sono circa la metà delle dimensioni dei pezzi quadrati, è possibile assegnare due a ciascun filo per bilanciare il carico - per un totale di k*k/2 attività approssimativamente uguali.

Un vantaggio di questo metodo è che ogni attività ha solo bisogno di accedere ai dati per i corpi , il che potrebbe rendere notevolmente più cache-friendly.

Answer 3

Oggi ho appena trovato la soluzione. Non sto ad accettare finché qualcuno lo conferma

Al fine di f [a_k] essere una funzione costante rispetto a k, poi

f [a_ {k + 1}] - f [a_k] = 0

deve essere vero per k = 1,2,3, ..., p-1.

Possiamo espandere questa equazione usando le definizioni che ho postato sulla domanda, e arriviamo ad un sistema di "p" 2º di equazioni algebriche di ordine rispetto ad a_k, k = 1,2,3, ... , p. Non vedo una soluzione chiusa a un p arbitrario, ma può essere risolta analiticamente per ogni p.

ho confermato che:

1. la somma, quando si utilizza l'a_k Ho calcolato era n (n-1)/2, il numero totale di interazioni di questo problema.

2. il numero di interazioni per thread è effettivamente costante per p = 2,3,4,5 e 10 (dove p = 10 ha richiesto del tempo per calcolare su mathematica®).

EDIT

Dopo l'ispezione dettagliata delle soluzioni per diversi valori di p, ho raggiunto al generale chiusa soluzione

a_k = 1/(2 p) (-p + 2 pn - sqrt [p^2 + 4 p (p + 1 - k) (n - 1) n])

valido per ogni p> = 2, n> 1.

Questo completa la risposta.

Answer 4

Supponendo vostro compilatore supporta OpenMP, perché non si può semplicemente cercare di fare

#pragma omp parallel for schedule(dynamic) // or: schedule(guided) 
for (i = 0; i < n; i++) 
    for (j = i+1; j < n; j++) 
     // calculate interaction 

o anche (avrete bisogno di punto di riferimento per capire quale funziona meglio)

#pragma omp parallel 
const int stride = omp_get_num_threads() + 1; 
for (i = omp_get_thread_num(); i < n; i += stride) 
    for (j = i+1; j < n; j++) 
     // calculate interaction