Perché questo codice non scala in modo lineare?

Question

Ho scritto questo codice SOR solver. Non preoccupatevi troppo di ciò che fa questo algoritmo, non è la preoccupazione qui. Ma solo per completezza: può risolvere un sistema lineare di equazioni, a seconda di quanto sia ben condizionato il sistema.

Lo eseguo con una matrice sparso di righe 2097152 mal condizionate (che non converge mai), con al massimo 7 colonne diverse da zero per riga.

Traducendo: l'esterno do-while ciclo eseguirà 10000 iterazioni (il valore passo come max_iters), mezzo for eseguirà 2097152 iterazioni, diviso in blocchi di work_line, suddivisi tra i fili OpenMP. Il ciclo più interno di for avrà 7 iterazioni, tranne in pochissimi casi (meno dell'1%) in cui può essere inferiore.

C'è una dipendenza tra i thread nei valori della matrice sol. Ogni iterazione del mezzo for aggiorna un elemento ma legge fino a 6 altri elementi dell'array. Poiché SOR non è un algoritmo esatto, durante la lettura, può avere qualsiasi valore precedente o attuale su quella posizione (se si ha familiarità con i risolutori, questo è un Gauss-Siedel che tollera il comportamento di Jacobi in alcuni punti per il bene di parallelismo).

typedef struct{ 
    size_t size; 

    unsigned int *col_buffer; 
    unsigned int *row_jumper; 
    real *elements; 
} Mat; 

int work_line; 

// Assumes there are no null elements on main diagonal 
unsigned int solve(const Mat* matrix, const real *rhs, real *sol, real sor_omega, unsigned int max_iters, real tolerance) 
{ 
    real *coefs = matrix->elements; 
    unsigned int *cols = matrix->col_buffer; 
    unsigned int *rows = matrix->row_jumper; 
    int size = matrix->size; 
    real compl_omega = 1.0 - sor_omega; 
    unsigned int count = 0; 
    bool done; 

    do { 
     done = true; 
     #pragma omp parallel shared(done) 
     { 
      bool tdone = true; 

      #pragma omp for nowait schedule(dynamic, work_line) 
      for(int i = 0; i < size; ++i) { 
       real new_val = rhs[i]; 
       real diagonal; 
       real residual; 
       unsigned int end = rows[i+1]; 

       for(int j = rows[i]; j < end; ++j) { 
        unsigned int col = cols[j]; 
        if(col != i) { 
         real tmp; 
         #pragma omp atomic read 
         tmp = sol[col]; 

         new_val -= coefs[j] * tmp; 
        } else { 
         diagonal = coefs[j]; 
        } 
       } 

       residual = fabs(new_val - diagonal * sol[i]); 
       if(residual > tolerance) { 
        tdone = false; 
       } 

       new_val = sor_omega * new_val/diagonal + compl_omega * sol[i]; 
       #pragma omp atomic write 
       sol[i] = new_val; 
      } 

      #pragma omp atomic update 
      done &= tdone; 
     } 
    } while(++count < max_iters && !done); 

    return count; 
} 

Come si può vedere, non v'è alcun blocco all'interno della regione in parallelo, così, per quello che ci insegnano sempre, è il tipo di problema al 100% in parallelo. Questo non è ciò che vedo in pratica.

Tutti i miei test sono stati eseguiti su CPU Intel (R) Xeon (R) E5-2670 v2 @ 2,50 GHz, 2 processori, 10 core ciascuno, abilitato per hyper-thread e somma di fino a 40 core logici.

Durante il mio primo set, work_line è stato corretto su 2048 e il numero di thread variava da 1 a 40 (40 corse in totale). Questo è il grafico del tempo di esecuzione di ogni ciclo (secondi x numero di thread):

La sorpresa era la curva logaritmica, così ho pensato che poiché la linea di lavoro era così grande, le cache condivise non erano molto ben utilizzati, quindi ho scavato questo file virtuale /sys/devices/system/cpu/cpu0/cache/index0/coherency_line_size che mi ha detto che la cache L1 di questo processore sincronizza gli aggiornamenti in gruppi di 64 byte (8 raddoppia nell'array sol). Così mi sono messo il work_line-8:

Poi ho pensato 8 era troppo basso per evitare stalli NUMA e impostare work_line-16:

Durante l'esecuzione di quanto sopra, ho pensato "Chi sono io per prevedere che cosa è work_line è buono? Basta vedere ...", e programmato per eseguire ogni work_line da 8 a 2048, passaggi di 8 (vale a dire ogni multiplo della linea della cache, da 1 a 256). I risultati per 20 e 40 fili (secondi x dimensioni della scissione della metà for ciclo, diviso tra i fili):

ritengo casi con bassa work_line soffre molto sincronizzazione cache, mentre grande work_line non offre alcun vantaggio oltre un certo numero di thread (presumo perché il percorso della memoria è il collo di bottiglia).È molto triste che un problema che sembra parallelo al 100% presenti un comportamento così negativo su una macchina reale. Quindi, prima di essere convinto che i sistemi multi-core sono una bugia molto ben venduta, vi chiedo innanzitutto:

Come posso rendere questo codice in scala lineare al numero di core? Cosa mi manca? C'è qualcosa nel problema che lo rende non buono come sembra in un primo momento?

Aggiornamento

Seguendo i suggerimenti, ho provato sia con static e dynamic programmazione, ma la rimozione delle atomiche di lettura/scrittura sulla matrice sol. Come riferimento, le linee blu e arancione sono le stesse del grafico precedente (fino a work_line = 248;). Le linee gialle e verdi sono le nuove. Per quello che ho potuto vedere: static fa una differenza significativa per il basso work_line, ma dopo 96 i benefici di dynamic superano il suo sovraccarico, rendendolo più veloce. Le operazioni atomiche non fanno alcuna differenza.

Answer 1

La moltiplicazione vettoriale a matrice sparsa è vincolata alla memoria (vedere here) e potrebbe essere visualizzata con un modello di linea semplice. I problemi legati alla memoria traggono vantaggio dalla maggiore larghezza di banda della memoria dei sistemi NUMA multisocket, ma solo se l'inizializzazione dei dati viene eseguita in modo tale che i dati vengano distribuiti tra i due domini NUMA. Ho alcune ragioni per credere che stai caricando la matrice in seriale e quindi tutta la sua memoria è allocata su un singolo nodo NUMA. In questo caso non sarà possibile beneficiare della larghezza di banda doppia memoria disponibile su un sistema dual-socket e 'davvero non importa se si utilizza schedule(dynamic) o schedule(static). Quello che potresti fare è abilitare la politica NUMA di interleaving della memoria per fare in modo che l'allocazione di memoria si diffonda tra i due nodi NUMA. In questo modo ogni thread finirà con il 50% di accesso alla memoria locale e il 50% di accesso alla memoria remota invece di avere tutti i thread sulla seconda CPU colpiti dal 100% dell'accesso remoto alla memoria. Il modo più semplice per attivare il criterio è quello di utilizzare numactl:

$ OMP_NUM_THREADS=... OMP_PROC_BIND=1 numactl --interleave=all ./program ... 

OMP_PROC_BIND=1 consente filo pinning e dovrebbe migliorare le prestazioni un po '.

Vorrei anche sottolineare che questo:

done = true; 
#pragma omp parallel shared(done) 
{ 
    bool tdone = true; 

    // ... 

    #pragma omp atomic update 
    done &= tdone; 
} 

è un probabilmente un non molto efficiente ri-attuazione:

done = true; 
#pragma omp parallel reduction(&:done) 
{ 
    // ... 
     if(residual > tolerance) { 
      done = false; 
     } 
    // ... 
} 

Non avrà una differenza di prestazioni notevole fra le due implementazioni a causa della quantità di lavoro svolto nel ciclo interno, ma non è comunque una buona idea reimplementare i primitivi OpenMP esistenti per motivi di portabilità e leggibilità.

Answer 2

vostro ciclo interno ha un omp atomic read, e il ciclo di mezzo ha un omp atomic write ad una posizione che potrebbe essere lo stesso letto da uno dei legge. OpenMP è obbligato a garantire che le scritture e le letture atomiche della stessa posizione siano serializzate, quindi in effetti probabilmente è necessario introdurre un blocco, anche se non ce n'è uno esplicito.

Potrebbe anche essere necessario bloccare l'intera matrice sol a meno che non sia in grado di capire in qualche modo quali letture potrebbero essere in conflitto con le scritture e in realtà i processori OpenMP non sono necessariamente così intelligenti.

Nessun codice scala in modo assolutamente lineare, ma state certi che ci sono molti codici che sono molto più vicini alla linearità rispetto ai vostri.

Answer 3

Sospetto si stiano riscontrando problemi di memorizzazione nella cache. Quando un thread aggiorna un valore nella matrice sol, esso invalidi le cache su altre CPU che archiviano stessa linea di cache. Ciò impone l'aggiornamento delle cache, che quindi porta allo stallo delle CPU.

Answer 4

Provare a eseguire l'IPCM (Intel Performance Counter Monitor). Puoi guardare la larghezza di banda della memoria e vedere se raggiunge il massimo con più core. La mia sensazione istintiva è che la larghezza di banda della memoria è limitata.

Come un rapido ritorno del calcolo dell'inviluppo, trovo che la larghezza di banda di lettura non collegata è di circa 10 GB/s su un Xeon. Se il tuo orologio è a 2.5 GHz, questa è una parola a 32 bit per ciclo di clock. Il tuo ciclo interno è fondamentalmente solo un'operazione di aggiunta multipla i cui cicli puoi contare su una mano, più alcuni cicli per l'overhead del loop. Non mi sorprende che dopo 10 thread, non ottieni alcun guadagno in termini di prestazioni.

Answer 5

Anche se non si dispone di un blocco mutex esplicito nel codice, si dispone di una risorsa condivisa tra i processi: la memoria e il relativo bus. Non lo vedi nel tuo codice perché è l'hardware che si occupa di gestire tutte le diverse richieste dalle CPU, ma tuttavia è una risorsa condivisa.

Così, ogni volta che uno dei processi scrive in memoria, quella posizione di memoria dovrà essere ricaricata dalla memoria principale da tutti gli altri processi che la usano, e tutti devono usare lo stesso bus di memoria per farlo. Il bus di memoria si satura e non si ottengono ulteriori guadagni in termini di prestazioni da core CPU aggiuntivi che servono solo a peggiorare la situazione.