Creazione di piani FFTW tramite OpenMP

Question

Sto cercando di eseguire più FFT in parallelo. Sto usando FFTW e OpenMP. Ogni FFT è diversa, quindi non sto facendo affidamento sul multithreading incorporato di FFTW (che conosco utilizza OpenMP).

int m; 

// assume: 
// int numberOfColumns = 100; 
// int numberOfRows = 100; 

#pragma omp parallel for default(none) private(m) shared(numberOfColumns, numberOfRows)// num_threads(4) 
    for(m = 0; m < 36; m++){ 

     // create pointers 
     double   *inputTest; 
     fftw_complex *outputTest; 
     fftw_plan  testPlan; 

     // preallocate vectors for FFTW 
     outputTest = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*numberOfRows*numberOfColumns); 
     inputTest = (double *)fftw_malloc(sizeof(double)*numberOfRows*numberOfColumns); 

     // confirm that preallocation worked 
     if (inputTest == NULL || outputTest == NULL){ 
      logger_.log_error("\t\t FFTW memory not allocated on m = %i", m); 
     } 

     // EDIT: insert data into inputTest 
     inputTest = someDataSpecificToThisIteration(m); // same size for all m 

     // create FFTW plan 
     #pragma omp critical (make_plan) 
     { 
      testPlan = fftw_plan_dft_r2c_2d(numberOfRows, numberOfColumns, inputTest, outputTest, FFTW_ESTIMATE); 
     } 

     // confirm that plan was created correctly 
     if (testPlan == NULL){ 
      logger_.log_error("\t\t failed to create plan on m = %i", m); 
     } 

     // execute plan 
     fftw_execute(testPlan); 

     // clean up 
     fftw_free(inputTest); 
     fftw_free(outputTest); 
     fftw_destroy_plan(testPlan); 

    }// end parallelized for loop 

Tutto funziona correttamente. Tuttavia, se rimuovo il costrutto critico dalla creazione del piano (fftw_plan_dft_r2c_2d) il mio codice fallirà. Qualcuno può spiegare perché? fftw_plan_dft_r2c_2d non è davvero un "orfano", giusto? È perché due thread potrebbero entrambi tentare di raggiungere la posizione di memoria numeroOfRows o numberOfColumns allo stesso tempo?

Answer 1

E 'praticamente tutto scritto nella documentazione FFTW su thread safety:

... ma una certa cura deve essere presa perché le routine planner condividere i dati (ad esempio, la saggezza e tabelle trigonometriche) tra le chiamate e progetti.

Il risultato è che l'unica routine thread-safe (rientranti) in FFTW è fftw_execute (e le varianti del nuovo array). Tutte le altre routine (ad esempio il pianificatore) devono essere chiamate da un thread alla volta. Quindi, ad esempio, puoi avvolgere un blocco semaforo attorno a qualsiasi chiamata al planner; ancora più semplicemente, puoi semplicemente creare tutti i tuoi piani da un thread. Non pensiamo che questa dovrebbe essere una restrizione importante (FFTW è progettato per la situazione in cui l'unico codice sensibile alla performance è l'effettiva esecuzione della trasformazione) e i benefici dei dati condivisi tra i piani sono ottimi.

In una tipica applicazione di piani FFT vengono costruiti raramente, quindi non importa se è necessario sincronizzare la loro creazione. Nel tuo caso non è necessario creare un nuovo piano ad ogni iterazione, a meno che la dimensione dei dati non cambi. Si preferisce effettuare le seguenti operazioni:

#pragma omp parallel default(none) private(m) shared(numberOfColumns, numberOfRows) 
{ 
    // create pointers 
    double   *inputTest; 
    fftw_complex *outputTest; 
    fftw_plan  testPlan; 

    // preallocate vectors for FFTW 
    outputTest = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*numberOfRows*numberOfColumns); 
    inputTest = (double *)fftw_malloc(sizeof(double)*numberOfRows*numberOfColumns); 

    // confirm that preallocation worked 
    if (inputTest == NULL || outputTest == NULL){ 
     logger_.log_error("\t\t FFTW memory not allocated on m = %i", m); 
    } 

    // create FFTW plan 
    #pragma omp critical (make_plan) 
    testPlan = fftw_plan_dft_r2c_2d(numberOfRows, numberOfColumns, inputTest, outputTest, FFTW_ESTIMATE); 

    #pragma omp for 
    for (m = 0; m < 36; m++) { 
     // execute plan 
     fftw_execute(testPlan); 
    } 

    // clean up 
    fftw_free(inputTest); 
    fftw_free(outputTest); 
    fftw_destroy_plan(testPlan); 
} 

Ora i piani vengono creati solo una volta in ogni filo e l'overhead di serializzazione diminuirebbe con ogni esecuzione di fftw_execute(). Se si esegue su un sistema NUMA (ad esempio un sistema multi-socket AMD64 o Intel (post-) Nehalem), è necessario abilitare il binding del thread per ottenere le massime prestazioni.