La regola d'oro di OpenMP è che tutte le variabili (con alcune esclusioni), che sono definite in un ambito esterno, sono condivise di default nella regione parallela. Poiché in Fortran prima del 2008 non ci sono ambiti locali (cioè non c'è lo BLOCK ... END BLOCK nelle versioni precedenti), tutte le variabili (eccetto quelle threadprivate) sono condivise, il che è molto naturale per me (a differenza di Ian Bush, non sono un grande fan dell'uso default(none) e quindi ridichiarare la visibilità di tutte le oltre 100 variabili locali in vari codici scientifici complessi).
Ecco come determinare la classe di condivisione di ogni variabile:
N - condivisi, perché dovrebbe essere la stessa in tutte le discussioni e hanno solo leggere il suo valore.ii - è il contatore del ciclo, soggetto a una direttiva di condivisione del lavoro, quindi la sua classe di condivisione è predeterminata per essere private. Non fa male dichiararlo esplicitamente in una clausola PRIVATE, ma non è davvero necessario.jj - contatore di loop di un loop, che non è soggetto a una direttiva di condivisione del lavoro, quindi jj deve essere private.X - condiviso, perché tutti i thread fanno riferimento e leggono solo da esso.distance_vector - ovviamente dovrebbe essere private poiché ogni thread funziona su diverse coppie di particelle.distance, distance2 e coff - idem.M - deve essere condiviso per gli stessi motivi di X.PE - funge da variabile di accumulo (immagino che questa sia l'energia potenziale del sistema) e dovrebbe essere oggetto di un'operazione di riduzione, vale a dire dovrebbe essere inserito in una clausola REDUCTION(+:....).A - questo è difficile. Può essere condiviso o aggiornato a A(jj,:) protetto con i costrutti di sincronizzazione, oppure è possibile utilizzare la riduzione (OpenMP consente riduzioni su variabili di array in Fortran a differenza di C/C++). A(ii,:) non viene mai modificato da più di un thread, quindi non richiede un trattamento speciale.
Con riduzione nel A in atto, ogni thread otterrebbe la sua copia privata di A e questo potrebbe essere un ingordo di memoria, anche se dubito si può usare questo O diretta (N) codice di simulazione per calcolare i sistemi con un numero molto grande di particelle. C'è anche un certo overhead associato all'implementazione della riduzione. In questo caso è sufficiente aggiungere A all'elenco della clausola REDUCTION(+:...).
Con i costrutti di sincronizzazione sono disponibili due opzioni. È possibile utilizzare il costrutto ATOMIC o il costrutto CRITICAL. Come ATOMIC è applicabile solo a contesti scalari, si dovrà "unvectorise" il ciclo di assegnazione e applicare ATOMIC ad ogni dichiarazione a parte, ad esempio:
!$OMP ATOMIC UPDATE
A(jj,1)=A(jj,1)+(M(ii)/coff)*(distance_vector(1))
!$OMP ATOMIC UPDATE
A(jj,2)=A(jj,2)+(M(ii)/coff)*(distance_vector(2))
!$OMP ATOMIC UPDATE
A(jj,3)=A(jj,3)+(M(ii)/coff)*(distance_vector(3))
Si può anche riscrivere questo come un loop - non dimenticate di dichiarare il contatore di loop private.
Con CRITICAL non v'è alcuna necessità di unvectorise loop:
!$OMP CRITICAL (forceloop)
A(jj,:)=A(jj,:)+(M(ii)/coff)*(distance_vector)
!$OMP END CRITICAL (forceloop)
Naming regioni critiche è opzionale e un po 'inutile, in questo caso particolare, ma in generale si permette di separare regioni critiche indipendenti.
Quale è più veloce? Srotolato con ATOMIC o CRITICAL? Dipende da molte cose. Solitamente lo CRITICAL è molto più lento poiché spesso implica chiamate di funzione al runtime di OpenMP mentre gli incrementi atomici, almeno su x86, sono implementati con istruzioni di aggiunta bloccate. Come spesso dicono, YMMV.
ricapitolare, una versione funzionante del ciclo dovrebbe essere qualcosa di simile:
!$OMP PARALLEL DO PRIVATE(jj,kk,distance_vector,distance2,distance,coff) &
!$OMP& REDUCTION(+:PE)
do ii=1,N-1
do jj=ii+1,N
distance_vector=X(ii,:)-X(jj,:)
distance2=sum(distance_vector*distance_vector)
distance=DSQRT(distance2)
coff=distance*distance*distance
PE=PE-M(II)*M(JJ)/distance
do kk=1,3
!$OMP ATOMIC UPDATE
A(jj,kk)=A(jj,kk)+(M(ii)/coff)*(distance_vector(kk))
end do
A(ii,:)=A(ii,:)-(M(jj)/coff)*(distance_vector)
end do
end do
!$OMP END PARALLEL DO
ho pensato che il vostro sistema è 3-dimensionale.
Con tutto questo detto, io secondo Ian Bush che è necessario ripensare il modo di posizione e le matrici di accelerazione sono disposti in memoria. L'utilizzo corretto della cache potrebbe aumentare il codice e consentirebbe anche alcune operazioni, ad es. X(:,ii)-X(:,jj) da vettorizzare, cioè implementato usando le istruzioni vettoriali SIMD.
Grazie ancora Hristo, di nuovo. Solo una domanda molto semplice: perché deve essere così complicata per cicli più semplici. Voglio dire, vuoi solo che il ciclo più esterno sia diviso in modi nthread e popolare una matrice (almeno nelle cose che faccio). Posso capire perché hai bisogno di privati, ecc. Per casi complessi e quindi è un vaso di 'ingegneria extra', ma perché non lo stato! $ OMP PARALLELL DO SHARED (MATRIX) quindi divide solo N/nthreads e consegna. Tutto in ogni thread sarà fatto in sequenza. Immagino di star pensando: perché non rendere tutto privato di default, non condiviso, quindi basta impostare la condivisione. – Griff
@Griff, in C/C++ queste variabili possono essere dichiarate all'interno del ciclo più interno e quindi sono automaticamente private. Il fatto che siano condivisi deriva dal desiderio di far coincidere la semantica di OpenMP con quella del codice seriale. A proposito, puoi sbarazzarti degli atomici se ti dimentichi per un secondo della terza legge di Newton ed esegui entrambi i cicli su '1, N' con assegnamento a' A (ii, :) 'solo (comunque inefficace). –
C'è comunque una possibilità per mostrarti un nuovo codice senza doverlo postare come domanda da controllare? È un po 'più semplice, ma più lungo e voglio solo controllarlo per assicurarmi di averlo capito correttamente a parole. Grazie. – Griff