Multithreaded e SIMD hanno vettorializzato Mandelbrot in R usando Rcpp & OpenMP

Come test delle prestazioni OpenMP e Rcpp , volevo verificare quanto velocemente avrei potuto calcolare il set Mandelbrot in R usando l’implementazione Rcpp + OpenMP più semplice e semplice. Attualmente quello che ho fatto è stato:

 #include  #include  // [[Rcpp::plugins(openmp)]] using namespace Rcpp; // [[Rcpp::export]] Rcpp::NumericMatrix mandelRcpp(const double x_min, const double x_max, const double y_min, const double y_max, const int res_x, const int res_y, const int nb_iter) { Rcpp::NumericMatrix ret(res_x, res_y); double x_step = (x_max - x_min) / res_x; double y_step = (y_max - y_min) / res_y; int r,c; #pragma omp parallel for default(shared) private(c) schedule(dynamic,1) for (r = 0; r < res_y; r++) { for (c = 0; c < res_x; c++) { double zx = 0.0, zy = 0.0, new_zx; double cx = x_min + c*x_step, cy = y_min + r*y_step; int n = 0; for (n=0; (zx*zx + zy*zy < 4.0 ) && ( n < nb_iter ); n++ ) { new_zx = zx*zx - zy*zy + cx; zy = 2.0*zx*zy + cy; zx = new_zx; } ret(c,r) = n; } } return ret; } 

E poi in R:

 library(Rcpp) sourceCpp("mandelRcpp.cpp") xlims=c(-0.74877,-0.74872); ylims=c(0.065053,0.065103); x_res=y_res=1080L; nb_iter=10000L; system.time(m <- mandelRcpp(xlims[[1]], xlims[[2]], ylims[[1]], ylims[[2]], x_res, y_res, nb_iter)) # 0.92s rainbow=c(rgb(0.47,0.11,0.53),rgb(0.27,0.18,0.73),rgb(0.25,0.39,0.81),rgb(0.30,0.57,0.75),rgb(0.39,0.67,0.60),rgb(0.51,0.73,0.44),rgb(0.67,0.74,0.32),rgb(0.81,0.71,0.26),rgb(0.89,0.60,0.22),rgb(0.89,0.39,0.18),rgb(0.86,0.13,0.13)) cols=c(colorRampPalette(rainbow)(100),rev(colorRampPalette(rainbow)(100)),"black") # palette par(mar=c(0, 0, 0, 0)) system.time(image(m^(1/7), col=cols, asp=diff(ylims)/diff(xlims), axes=F, useRaster=T)) # 0.5s 

inserisci la descrizione dell'immagine qui

Non ero sicuro se ci fossero altri evidenti miglioramenti della velocità che potrei sfruttare oltre al multithreading di OpenMP, ad esempio tramite la vettorizzazione simd ? (usando le opzioni di simd nel file #pragma non sembrava fare nulla)

PS all’inizio il mio codice stava andando in crash, ma in seguito ho trovato che questo era stato risolto sostituendo ret[r,c] = n; con ret(r,c) = n; Usare le lezioni di Armadillo come suggerito nella risposta qui sotto rende le cose leggermente più veloci, sebbene i tempi siano quasi gli stessi. Inoltre, ruotando attorno a y si ottiene un orientamento corretto quando si stampa con l’ image() . L’uso della velocità di 8 thread è ca. 350 volte più veloce della semplice versione R Mandelbrot vettorizzata qui e anche circa 7,3 volte più veloce rispetto alla versione Python / Numba (senza multithreading) qui (simile alle velocità PyCUDA o PyOpenCL), quindi abbastanza felice con questo … Rasterizzare / visualizzare ora sembra il collo di bottiglia in R ….

Non utilizzare OpenMP con gli oggetti *Vector o *Matrix SEXP mentre mascherano SEXP funzioni SEXP / allocazioni di memoria a thread singolo. OpenMP è un approccio multi-thread .

Questo è il motivo per cui il codice si sta bloccando.

Un modo per aggirare questa limitazione è usare una struttura di dati non R per memorizzare i risultati. Uno dei seguenti sarà sufficiente: arma::mat o Eigen::MatrixXd o std::vector … Come preferisco l’armadillo, cambierò la matrice res in arma::mat da Rcpp::NumericMatrix . Pertanto, quanto segue eseguirà il tuo codice in parallelo:

 #include  // Note the changed include and new attribute // [[Rcpp::depends(RcppArmadillo)]] // Avoid including header if openmp not on system #ifdef _OPENMP #include  #endif // [[Rcpp::plugins(openmp)]] // Note the changed return type // [[Rcpp::export]] arma::mat mandelRcpp(const double x_min, const double x_max, const double y_min, const double y_max, const int res_x, const int res_y, const int nb_iter) { arma::mat ret(res_x, res_y); // note change double x_step = (x_max - x_min) / res_x; double y_step = (y_max - y_min) / res_y; unsigned r,c; #pragma omp parallel for shared(res) for (r = 0; r < res_y; r++) { for (c = 0; c < res_x; c++) { double zx = 0.0, zy = 0.0, new_zx; double cx = x_min + c*x_step, cy = y_min + r*y_step; unsigned n = 0; for (; (zx*zx + zy*zy < 4.0 ) && ( n < nb_iter ); n++ ) { new_zx = zx*zx - zy*zy + cx; zy = 2.0*zx*zy + cy; zx = new_zx; } if(n == nb_iter) { n = 0; } ret(r, c) = n; } } return ret; } 

Con il codice di test (nota y e x non sono stati definiti, quindi ho assunto y = ylims e x = xlims ) abbiamo:

 xlims = ylims = c(-2.0, 2.0) x_res = y_res = 400L nb_iter = 256L system.time(m <- mandelRcpp(xlims[[1]], xlims[[2]], ylims[[1]], ylims[[2]], x_res, y_res, nb_iter)) rainbow = c( rgb(0.47, 0.11, 0.53), rgb(0.27, 0.18, 0.73), rgb(0.25, 0.39, 0.81), rgb(0.30, 0.57, 0.75), rgb(0.39, 0.67, 0.60), rgb(0.51, 0.73, 0.44), rgb(0.67, 0.74, 0.32), rgb(0.81, 0.71, 0.26), rgb(0.89, 0.60, 0.22), rgb(0.89, 0.39, 0.18), rgb(0.86, 0.13, 0.13) ) cols = c(colorRampPalette(rainbow)(100), rev(colorRampPalette(rainbow)(100)), "black") # palette par(mar = c(0, 0, 0, 0)) image(m, col = cols, asp = diff(range(ylims)) / diff(range(xlims)), axes = F) 

Per:

inserisci la descrizione dell'immagine qui

Sono andato avanti e ho vettorializzato il codice dell’OP usando le estensioni vettoriali di GCC e Clang. Prima di mostrare come l’ho fatto, permettimi di mostrare le prestazioni con il seguente hardware:

 Skylake (SKL) at 3.1 GHz with 4 cores Knights Landing (KNL) at 1.5 GHz with 68 cores ARMv8 Cortex-A57 arch64 (Nvidia Jetson TX1) 4 cores at ? GHz nb_iter = 1000000 GCC Clang SKL_scalar 6m5,422s SKL_SSE41 3m18,058s SKL_AVX2 1m37,843s 1m39,943s SKL_scalar_omp 0m52,237s SKL_SSE41_omp 0m29,624s 0m31,356s SKL_AVX2_omp 0m14,156s 0m16,783s ARM_scalar 15m28.285s ARM_vector 9m26.384s ARM_scalar_omp 3m54.242s ARM_vector_omp 2m21.780s KNL_scalar 19m34.121s KNL_SSE41 11m30.280s KNL_AVX2 5m0.005s 6m39.568s KNL_AVX512 2m40.934s 6m20.061s KNL_scalar_omp 0m9.108s KNL_SSE41_omp 0m6.666s 0m6.992s KNL_AVX2_omp 0m2.973s 0m3.988s KNL_AVX512_omp 0m1.761s 0m3.335s 

L’accelerazione teorica di KNL vs. SKL è

 (68 cores/4 cores)*(1.5 GHz/3.1 Ghz)* (8 doubles per lane/4 doubles per lane) = 16.45 

Sono andato in dettaglio sulle funzionalità di estensioni vettoriali di GCC e Clang qui . Per vettorizzare il codice dell’OP qui ci sono tre operazioni vettoriali aggiuntive che dobbiamo definire.

1. Trasmissione

Per un vettore v e uno scalare s GCC non può fare v = s ma Clang può. Ma ho trovato una bella soluzione che funziona per GCC e Clang qui . Per esempio

 vsi v = s - (vsi){}; 

2. Una funzione any() come in OpenCL o come in R.

Il meglio che ho trovato è una funzione generica

 static bool any(vli const & x) { for(int i=0; i 

Clang genera in realtà un codice relativamente efficiente per questo usando l'istruzione ptest (ma non per AVX512 ) ma GCC no.

3. Compressione

I calcoli vengono eseguiti come doppi a 64 bit, ma il risultato viene scritto come numeri interi a 32 bit. Quindi due calcoli vengono eseguiti utilizzando numeri interi a 64 bit e quindi i due calcoli vengono compressi in un vettore di numeri interi a 32 bit. Ho trovato una soluzione generica con cui Clang fa un buon lavoro

 static vsi compress(vli const & lo, vli const & hi) { vsi lo2 = (vsi)lo, hi2 = (vsi)hi, z; for(int i=0; i 

La soluzione seguente funziona meglio per GCC ma non è migliore per Clang . Ma dal momento che questa funzione non è critica, io uso solo la versione generica.

 static vsi compress(vli const & low, vli const & high) { #if defined(__clang__) return __builtin_shufflevector((vsi)low, (vsi)high, MASK); #else return __builtin_shuffle((vsi)low, (vsi)high, (vsi){MASK}); #endif } 

Queste definizioni non si basano su qualcosa di specifico x86 e il codice (definito sotto) compila anche per i processori ARM con GCC e Clang.


Ora che questi sono definiti qui è il codice

 #include  #include  #include  using namespace Rcpp; #ifdef _OPENMP #include  #endif // [[Rcpp::plugins(openmp)]] // [[Rcpp::plugins(cpp14)]] #if defined ( __AVX512F__ ) || defined ( __AVX512__ ) static const int SIMD_SIZE = 64; #elif defined ( __AVX2__ ) static const int SIMD_SIZE = 32; #else static const int SIMD_SIZE = 16; #endif static const int VSI_SIZE = SIMD_SIZE/sizeof(int32_t); static const int VLI_SIZE = SIMD_SIZE/sizeof(int64_t); static const int VDF_SIZE = SIMD_SIZE/sizeof(double); #if defined(__clang__) typedef int32_t vsi __attribute__ ((ext_vector_type(VSI_SIZE))); typedef int64_t vli __attribute__ ((ext_vector_type(VLI_SIZE))); typedef double vdf __attribute__ ((ext_vector_type(VDF_SIZE))); #else typedef int32_t vsi __attribute__ ((vector_size (SIMD_SIZE))); typedef int64_t vli __attribute__ ((vector_size (SIMD_SIZE))); typedef double vdf __attribute__ ((vector_size (SIMD_SIZE))); #endif static bool any(vli const & x) { for(int i=0; i 

Il codice R è quasi lo stesso del codice dell'OP

 library(Rcpp) sourceCpp("frac.cpp", verbose=TRUE, rebuild=TRUE) xlims=c(-0.74877,-0.74872); ylims=c(0.065053,0.065103); x_res=y_res=1080L; nb_iter=100000L; t = system.time(m <- frac(xlims[[1]], xlims[[2]], ylims[[1]], ylims[[2]], x_res, y_res, nb_iter)) print(t) m2 = matrix(m, ncol = x_res) rainbow = c( rgb(0.47, 0.11, 0.53), rgb(0.27, 0.18, 0.73), rgb(0.25, 0.39, 0.81), rgb(0.30, 0.57, 0.75), rgb(0.39, 0.67, 0.60), rgb(0.51, 0.73, 0.44), rgb(0.67, 0.74, 0.32), rgb(0.81, 0.71, 0.26), rgb(0.89, 0.60, 0.22), rgb(0.89, 0.39, 0.18), rgb(0.86, 0.13, 0.13) ) cols = c(colorRampPalette(rainbow)(100), rev(colorRampPalette(rainbow)(100)),"black") # palette par(mar = c(0, 0, 0, 0)) image(m2^(1/7), col=cols, asp=diff(ylims)/diff(xlims), axes=F, useRaster=T) 

Per compilare GCC o Clang cambiare il file ~/.R/Makevars in

 CXXFLAGS= -Wall -std=c++14 -O3 -march=native -ffp-contract=fast -fopenmp #uncomment the following two lines for clang #CXX=clang-5.0 #LDFLAGS= -lomp 

Se hai problemi a far funzionare OpenMP per Clang, guarda questo .


Il codice produce più o meno la stessa immagine. inserisci la descrizione dell'immagine qui