[C/C++] Consigli ottimizzazione velocità

[cionci]

Originariamente inviato da Unrue

Purtroppo non conosco Visual Studio,ma comunque restrict è presente nello standard C99, prova a guardare se è già configurato per tale standard.

Il compilatore MS non è compatibile con C99.

[cionci]

Ci puoi fare vedere il codice ora ? Potremmo magari vedere se si può fare qualche trasformazione sul calcolo matriciale.

[Tommo]

Un ebook interessante: http://www.agner.org/optimize/optimizing_cpp.pdf

[Rossi88]

@Tommo
Grazie del link, leggerò sicuramente il libro

@Unrue
Sicuramente ho già ottenuto un più che ottimo risultato, però più ottimizzo più posso ottenere risultati accurati, insomma finchè si può... ottimizzo
Purtroppo il codice è parallelizzato invarie parti (elaborazione e salvataggio dei dati) facendo uso di funzioni specefiche per Windows, dovrei quindi riscriverlo, inoltre il codice dovrà girare anche su windows. Mi costerebbe un po' di fatica, comunque vedo che si può fare, magari riesco a sfruttare le MKL.
L'unroll del codice l'ho tolto in quanto mi da tempi analoghi (se non leggermente peggiori, forse il compilatore fa un unroll più efficiente, o il processore ha degli algoritmi predettivi efficaci per i salti nei loop).

@cionci
Ecco il codice:

void derivsU(double x , double y[], double dydx[],double lam)
{
	 int i;
	 int n=N;
	 int nquadro=n*n;

	 double Kmod=(1+fastcos(2*pi/lam*x));

	 double **xdbeta;
	 double **mcos, **msin;

	 xdbeta=matrixd(1,n,1,n);
	 mcos=matrixd(1,n,1,n);
	 msin=matrixd(1,n,1,n);

	 dcopy(nquadro,&beta[1][1],1,&xdbeta[1][1],1);

	 dscal(nquadro,x,&xdbeta[1][1],1);	 

	 vrda_sincos(nquadro,&xdbeta[1][1],&msin[1][1],&mcos[1][1]);

	 for (i=1;i<=nquadro;i++){
		 msin[1][i]*=intePerParam[1][i];
		 mcos[1][i]*=intePerParam[1][i];
	 }

	 dgemv('t',n,n,-Kmod,&mcos[1][1],n,&y[n+1],1,0,&dydx[1],1);
	 dgemv('t',n,n,-Kmod,&msin[1][1],n,&y[1],1,1,&dydx[1],1);

	 dgemv('t',n,n,Kmod,&mcos[1][1],n,&y[1],1,0,&dydx[1+n],1);
	 dgemv('t',n,n,-Kmod,&msin[1][1],n,&y[n+1],1,1,&dydx[n+1],1);

	 free_matrixd(xdbeta,1,n,1,n);
	 free_matrixd(mcos,1,n,1,n);
	 free_matrixd(msin,1,n,1,n);
}

Le matrici sono salvate in memoria in modo che tutte le righe siano memorizzate in spazi consecutivi di memoria, in questo modo posso trattarle indifferentemente come matrici o come vettori.

Essenzialmente il codice prevede la copia della matrice beta nella matrice xdbeta (dcopy) che poi moltiplico per x (dscal). poi calcolo il seno e il coseno di questa matrice (vrda_sincos), moltiplico elemento per elemento le due matrici ottenute per la matrice intePerParam (tramite il ciclo for).
Poi effettuo dei prodotti matriciali (o meglio matrice-vettore) e salvo il risultato in dydx.

Sicuramente può essere ottimizzato, magari evitando di creare e liberare le matrici ad ogni chiamata alla funzione, ma il punto lento è il calcolo del seno e del coseno. Ho utilizzato AMD CodeAnalyst e i risultati sono:

vrda_sincos 65% del tempo (funzione per il calcolo del seno e del coseno)
zupmtr_ circa 20% del tempo (dovrebbe essere una funzione che consente di salvare le matrici, io non faccio esplicita richiesta a questa funzione, sarà richiamata da un'altra funzione BLAS ma non so quale sia)
il mio programma occupa meno del 10% del tempo di cui il 13% nella creazione delle matrici, e il 70% nel fare il prodotto elemento per elemento. Qundi in parole pavore c'è da concentrarsi sul calcolo dei seni e dei coseni e su quella funzione il resto darebbe risultati pressochè trascurabili.

Sicuramente ho già ottenuto risultati ottimi, vedo solo se è possibile migliorare ulteriormente , devo informarmi meglio sulla possibilità di utilizzare le schede video che dovrebbero essere delle schegge nelle funzioni trigonometriche.


Per chi avesse bisogno di utilizzare le funzioni BLAS ottime guide sono anche le seguenti:

http://www.intel.com/software/produc.../docs/mklqref/

e una volta individuata la funzione altre informazioni (molto più dettagliate) sulla funzione e sui parametri possono essere trovati in

http://techpubs.sgi.com/library/tpl/...0650&pth=/cat3

[cionci]

Non conosco le grandezze in gioco. Ti ci mette tanto tempo perché devi eseguire tante volte quella funzione ?

[Unrue]

Originariamente inviato da Rossi88

Sicuramente ho già ottenuto risultati ottimi, vedo solo se è possibile migliorare ulteriormente , devo informarmi meglio sulla possibilità di utilizzare le schede video che dovrebbero essere delle schegge nelle funzioni trigonometriche.

Secondo me ti conviene calcolare i flop/s che ottieni con il tuo processore. In modo da valutare se sei arrivato al limite con il processore che hai.

Per fare questo, devi calcolare il numero di operazioni floating point che effettui e dividere per il tempo che impieghi. Poi lo confronti con flop/s teorici del tuo processore. Se raggiungi il 30/40% di questo valore, vuol dire che hai già raggiunto il massimo possibile con questo processore. Questo perché comunque intervengono altri fattori limitanti quali ad esempio la memoria.

Adesso però, avendo chiamate BLAS, è più complicato, dovresti vedere l'assembly generato e calcolarle da lì. Non è banalissimo, ma secondo me conviene farlo. Altrimenti rischi di continuare a cercare di ottimizzare quando poi scopri che più di lì è impossibile andare. Basta farlo facendo un ciclo di for. Altrimenti è un bordello.

[Rossi88]

Sono riuscito ad usare un trucco per eliminare il calcolo continuo dei seni e dei coseni ottenendo quindi un notevole boost prestazionale (precalcolo i seni e i coseni prima di dare in pasto al solutore che richiamerà milioni di volte la funzione, i seni e i coseni sono calcolati in modo esatto).
Il codice risulta quindi (ho inserito sia il codice con le librerie BLAS che quello senza)

void derivsU(double x , double y[], double dydx[],double lam)
{
#if !defined NOT_USE_ACML
	int n=N;
	int Indice=(int) floor(x*2/h+0.5)+1;
	double Kmod=(1+fastcos(2*pi/lam*x));
	
	dgemv('t',n,n,-Kmod,&ValoriPreCalcolatiCos[Indice][1][1],n,&y[n+1],1,0,&dydx[1],1);
	dgemv('t',n,n,-Kmod,&ValoriPreCalcolatiSin[Indice][1][1],n,&y[1],1,1,&dydx[1],1);

	dgemv('t',n,n,Kmod,&ValoriPreCalcolatiCos[Indice][1][1],n,&y[1],1,0,&dydx[1+n],1);
	dgemv('t',n,n,-Kmod,&ValoriPreCalcolatiSin[Indice][1][1],n,&y[n+1],1,1,&dydx[n+1],1);
#else
	int i,j;
	int n=N;
	int Indice=(int) floor(x*2/h+0.5)+1;
	double Kmod=(1+cos(2*pi/lam*x));
	double q1,q2;

	for(i=1;i<=n;i++){
		dydx[i]=0;
		dydx[i+n]=0;
		for(j=1;j<=n;j++)
		{
			q1=Kmod*ValoriPreCalcolatiSin[Indice][i][j];
			q2=Kmod*ValoriPreCalcolatiCos[Indice][i][j];

			dydx[i]-=q1*y[j]+q2*y[j+n];
			dydx[i+n]+=q2*y[j]-q1*y[j+n];
		}
	}
#endif
}

Ora per ottimizzare ulteriormente non mi resta che evitare anche di fare somme e moltiplicazioni ("dydx[i]-=q1*y[j]+q2*y[j+n];") e poi il programma sarà diventerà ancora più veloce .
Scherzo ovviamente, meglio di così non credo si possa fare per quella funzione.

Scrivo questo messaggio semplicemente per suggerire a chiunque avesse problemi di questo tipo di pensare in generale a tutto il programma nella sua interezza senza focalizzarsi eccessivamente sulla singola funzione. Nel mio caso ad esempio i sin e cos dipendevano sì dalla x e quindi variavano ad ogni chiamata ma gli stessi identici valori di x erano usati migliaia e migliaia di volte.

Per la cronaca il tutto è inserito in un ciclo che serve per risolvere un sistema di equazioni differenziali del prim'ordine, il solutore è un solutore Runge Kutta del 4° ordine RK4 a step costante, utilizzare un Runge kutta del 4° 5° RK45 a step variabile non risolve molto dato che il problema mi sembra si presenti moderatamente stiff, il passo di soluzione quindi resta piccolo (circa il doppio) e in più in questo caso dovrei calcolare sempre seni e coseni.
Se qualcuno avesse da suggerirmi un solutore più efficiente, mi sembra di aver capito che i solutori impliciti dovrebbero adattarsi meglio a quelli espliciti (come il RK) in queste situazioni, ben venga

Inoltre in determinati contesti solo una parte dell'elaborazione risulta essere stiff mentre la restante no, in questi casi quindi è opportuno cambiare solutore, nel mio caso passo dal RK4 al RK45 che consente di aumentare enormemente le prestazioni per la parte non stiff.


In definitiva quindi se non si può ottimizzare ulteriormente la funzione più onerosa occorre ingegnarsi per capire come è possibile ridurre il numero di chiamate a tale funzione.
Ad ogni modo mi posso ritenere più che soddisfatto.
Grazie a tutti.

[Unrue]

Ottimo!

Quanto è il boost rispetto al codice originale?

[Rossi88]

Originariamente inviato da Unrue

Ottimo!

Quanto è il boost rispetto al codice originale?

Con un test sintetico in un tipico scenario di utilizzo

senza BLAS e senza precalcolo 8'33''
senza BLAS e con precalcolo 45''
con BLAS e senza precalcolo 3'05''
con BLAS e con precalcolo 53''

il test l'ho fatto con tanti altri programmi avviati, che non facevano niente ma possono aver modificato leggermente i valori. Le librerie BLAS risultano quindi estremamente valide per velocizzare il calcolo dei sin e cos, per i prodotti matriciali penso che il compilatore di Visual Studio riesca ad ottimizzare maggiormente, anche se di poco.
Dando uno sguardo al codice assembly, da quel poco che ne capisco, usa funzioni tipo xmms o qualcosa del genere quindi sfrutta le istruzioni SSE che velocizzano, anche le librerie BLAS (per l'esattezza l'implementazione dell'AMD ossia ACML) sicuramente le sfruttano ma il dover richiamare ogni volta delle DLL esterne probabilmente causa quel leggero degrado prestazionale. Almeno credo.

[Unrue]

Originariamente inviato da Rossi88

Con un test sintetico in un tipico scenario di utilizzo

senza BLAS e senza precalcolo 8'33''
senza BLAS e con precalcolo 45''
con BLAS e senza precalcolo 3'05''
con BLAS e con precalcolo 53''

il test l'ho fatto con tanti altri programmi avviati, che non facevano niente ma possono aver modificato leggermente i valori. Le librerie BLAS risultano quindi estremamente valide per velocizzare il calcolo dei sin e cos, per i prodotti matriciali penso che il compilatore di Visual Studio riesca ad ottimizzare maggiormente, anche se di poco.
Dando uno sguardo al codice assembly, da quel poco che ne capisco, usa funzioni tipo xmms o qualcosa del genere quindi sfrutta le istruzioni SSE che velocizzano, anche le librerie BLAS (per l'esattezza l'implementazione dell'AMD ossia ACML) sicuramente le sfruttano ma il dover richiamare ogni volta delle DLL esterne probabilmente causa quel leggero degrado prestazionale. Almeno credo.

In realtà le BLAS sono molto famose per i prodotti matriciali, strano che vadano leggermente peggio. Forse aumentando la dimensionalità del test case andranno sicuramente meglio. Non ricordo se giò te lo consigliai, ma se sei su processori Intel, forse vale la pena provare anche limplementazione BLAS delle MKL, che sono ancora più spinte. ( le trovi insieme al compilatore Intel e la versione non pro è gratuita)

Hai controllato anche l'accuratezza dei risultati? Hai differenze tra gli output originali e quelli ottimizzati? Se si, di quanto?

PS:in realtà gli xmms sono registri vettoriali, non funzioni , utilizzati dalle istruzioni SSE.

[Rossi88]

i sin e cos precalcolati sono calcolati in maniera esatta, cioè con le normali funzioni sin e cos (o con vrda_sincos se uso ACML), quindi i risultati sono pressoché identici. La massima variazione che ho riscontrato nel test è di 1.2488e-008%, sarà dovuto a qualche errore di troncamento sulle ultime cifre.

Le MKL mi sembrano che siano gratis solo per Linux , e il mio programma fa uso di funzioni di parallelizzazioni scritte per windows.

Comunque il codice Assembly che Visual Studio genera per quella funzione è il seguente:

249: void derivsU(double x , double y[], double dydx[],double lam)
   250: {
000000013F3F1300 40 53                push        rbx  
000000013F3F1302 57                   push        rdi  
000000013F3F1303 41 55                push        r13  
000000013F3F1305 41 57                push        r15  
000000013F3F1307 48 83 EC 68          sub         rsp,68h  
000000013F3F130B 0F 29 7C 24 30       movaps      xmmword ptr [rsp+30h],xmm7  
000000013F3F1310 0F 29 74 24 40       movaps      xmmword ptr [rsp+40h],xmm6  
000000013F3F1315 66 0F 28 F0          movapd      xmm6,xmm0  
   251: #if !defined NOT_USE_ACML
   252: 	/* La parte iniziale del vettore contiene i dati legati alla parte
   253: 	* reale dell'ampiezza, mentre la seconda parte contiene la parte
   254: 	* immaginaria.
   255: 	*/
   256: 	int n=N;
   257: 	int Indice=(int) floor(x*2/h+0.5)+1;
   258: 	double Kmod=(1+fastcos(2*pi/lam*x));
   259: 	
   260: 	dgemv('t',n,n,-Kmod,&ValoriPreCalcolatiCos[Indice][1][1],n,&y[n+1],1,0,&dydx[1],1);
   261: 	dgemv('t',n,n,-Kmod,&ValoriPreCalcolatiSin[Indice][1][1],n,&y[1],1,1,&dydx[1],1);
   262: 
   263: 	dgemv('t',n,n,Kmod,&ValoriPreCalcolatiCos[Indice][1][1],n,&y[1],1,0,&dydx[1+n],1);
   264: 	dgemv('t',n,n,-Kmod,&ValoriPreCalcolatiSin[Indice][1][1],n,&y[n+1],1,1,&dydx[n+1],1);
   265: #else
   266: 	int i,j;
   267: 	int n=N;
   268: 	int Indice=(int) floor(x*2/h+0.5)+1;
   269: 	double Kmod=(1+cos(2*pi/lam*x));
000000013F3F1319 F2 0F 10 05 CF 79 03 00 movsd       xmm0,mmword ptr [__real@401921fb54442d18 (13F428CF0h)]  
000000013F3F1321 4D 8B F8             mov         r15,r8  
000000013F3F1324 48 8B FA             mov         rdi,rdx  
000000013F3F1327 F2 0F 5E C3          divsd       xmm0,xmm3  
000000013F3F132B F2 0F 59 C6          mulsd       xmm0,xmm6  
000000013F3F132F E8 1C 3D 04 00       call        cos (13F435050h)  
000000013F3F1334 F2 0F 59 35 44 79 03 00 mulsd       xmm6,mmword ptr [__real@4000000000000000 (13F428C80h)]  
000000013F3F133C 48 63 1D 3D F1 03 00 movsxd      rbx,dword ptr [N (13F430480h)]  
000000013F3F1343 F2 0F 5E 35 7D 05 04 00 divsd       xmm6,mmword ptr [h (13F4318C8h)]  
000000013F3F134B 66 0F 28 F8          movapd      xmm7,xmm0  
000000013F3F134F F2 0F 58 3D B9 D0 03 00 addsd       xmm7,mmword ptr [VersioneFileSalvataggio+20h (13F42E410h)]  
000000013F3F1357 F2 0F 58 35 81 79 03 00 addsd       xmm6,mmword ptr [__real@3fe0000000000000 (13F428CE0h)]  
000000013F3F135F 66 0F 28 C6          movapd      xmm0,xmm6  
000000013F3F1363 E8 F8 86 02 00       call        floor (13F419A60h)  
000000013F3F1368 0F 28 74 24 40       movaps      xmm6,xmmword ptr [rsp+40h]  
000000013F3F136D F2 0F 2C C0          cvttsd2si   eax,xmm0  
000000013F3F1371 FF C0                inc         eax  
000000013F3F1373 4C 63 E8             movsxd      r13,eax  
   270: 	double q1,q2;
   271: 
   272: 	for(i=1;i<=n;i++){
000000013F3F1376 48 83 FB 01          cmp         rbx,1  
000000013F3F137A 0F 8C 34 03 00 00    jl          derivsU+3B4h (13F3F16B4h)  
   251: #if !defined NOT_USE_ACML
   252: 	/* La parte iniziale del vettore contiene i dati legati alla parte
   253: 	* reale dell'ampiezza, mentre la seconda parte contiene la parte
   254: 	* immaginaria.
   255: 	*/
   256: 	int n=N;
   257: 	int Indice=(int) floor(x*2/h+0.5)+1;
   258: 	double Kmod=(1+fastcos(2*pi/lam*x));
   259: 	
   260: 	dgemv('t',n,n,-Kmod,&ValoriPreCalcolatiCos[Indice][1][1],n,&y[n+1],1,0,&dydx[1],1);
   261: 	dgemv('t',n,n,-Kmod,&ValoriPreCalcolatiSin[Indice][1][1],n,&y[1],1,1,&dydx[1],1);
   262: 
   263: 	dgemv('t',n,n,Kmod,&ValoriPreCalcolatiCos[Indice][1][1],n,&y[1],1,0,&dydx[1+n],1);
   264: 	dgemv('t',n,n,-Kmod,&ValoriPreCalcolatiSin[Indice][1][1],n,&y[n+1],1,1,&dydx[n+1],1);
   265: #else
   266: 	int i,j;
   267: 	int n=N;
   268: 	int Indice=(int) floor(x*2/h+0.5)+1;
   269: 	double Kmod=(1+cos(2*pi/lam*x));
000000013F3F1380 48 89 AC 24 98 00 00 00 mov         qword ptr [y],rbp  
000000013F3F1388 48 89 74 24 60       mov         qword ptr [rsp+60h],rsi  
000000013F3F138D 4D 8D 4F 08          lea         r9,[r15+8]  
000000013F3F1391 4C 89 64 24 58       mov         qword ptr [rsp+58h],r12  
   273: 		dydx[i]=0;
000000013F3F1396 49 F7 DF             neg         r15  
000000013F3F1399 4C 8D 14 DD 00 00 00 00 lea         r10,[rbx*8]  
000000013F3F13A1 4C 89 74 24 50       mov         qword ptr [rsp+50h],r14  
000000013F3F13A6 48 8B EB             mov         rbp,rbx  
000000013F3F13A9 48 89 9C 24 90 00 00 00 mov         qword ptr [x],rbx  
000000013F3F13B1 33 C0                xor         eax,eax  
000000013F3F13B3 66 66 66 66 66 0F 1F 84 00 00 00 00 00 nop         word ptr [rax+rax]  
000000013F3F13C0 49 89 01             mov         qword ptr [r9],rax  
   274: 		dydx[i+n]=0;
   275: 		for(j=1;j<=n;j++)
000000013F3F13C3 BE 01 00 00 00       mov         esi,1  
000000013F3F13C8 4B 89 04 0A          mov         qword ptr [r10+r9],rax  
000000013F3F13CC 48 83 FB 04          cmp         rbx,4  
000000013F3F13D0 0F 8C 00 02 00 00    jl          derivsU+2D6h (13F3F15D6h)  
   270: 	double q1,q2;
   271: 
   272: 	for(i=1;i<=n;i++){
000000013F3F13D6 48 8B 05 33 06 04 00 mov         rax,qword ptr [ValoriPreCalcolatiSin (13F431A10h)]  
000000013F3F13DD 4B 8D 0C 0F          lea         rcx,[r15+r9]  
000000013F3F13E1 4C 8B F3             mov         r14,rbx  
000000013F3F13E4 4A 8B 04 E8          mov         rax,qword ptr [rax+r13*8]  
000000013F3F13E8 4D 8B E2             mov         r12,r10  
000000013F3F13EB 49 F7 DE             neg         r14  
000000013F3F13EE 4C 8B 04 01          mov         r8,qword ptr [rcx+rax]  
000000013F3F13F2 48 8B 05 1F 06 04 00 mov         rax,qword ptr [ValoriPreCalcolatiCos (13F431A18h)]  
000000013F3F13F9 49 F7 DC             neg         r12  
000000013F3F13FC 4A 8B 04 E8          mov         rax,qword ptr [rax+r13*8]  
   273: 		dydx[i]=0;
000000013F3F1400 48 8D 6C DF 10       lea         rbp,[rdi+rbx*8+10h]  
000000013F3F1405 48 8B 14 01          mov         rdx,qword ptr [rcx+rax]  
000000013F3F1409 49 8B C8             mov         rcx,r8  
000000013F3F140C B8 02 00 00 00       mov         eax,2  
000000013F3F1411 48 2B C3             sub         rax,rbx  
000000013F3F1414 4C 8B DA             mov         r11,rdx  
000000013F3F1417 49 2B CA             sub         rcx,r10  
000000013F3F141A 48 C1 E0 03          shl         rax,3  
000000013F3F141E 4D 2B DA             sub         r11,r10  
000000013F3F1421 48 2B CF             sub         rcx,rdi  
000000013F3F1424 48 89 84 24 A8 00 00 00 mov         qword ptr [lam],rax  
000000013F3F142C 48 2B C7             sub         rax,rdi  
000000013F3F142F 4C 2B DF             sub         r11,rdi  
000000013F3F1432 4C 03 C0             add         r8,rax  
000000013F3F1435 48 03 C2             add         rax,rdx  
000000013F3F1438 48 8B 94 24 A8 00 00 00 mov         rdx,qword ptr [lam]  
000000013F3F1440 48 89 44 24 20       mov         qword ptr [rsp+20h],rax  
000000013F3F1445 48 8B C3             mov         rax,rbx  
000000013F3F1448 48 8B 5C 24 20       mov         rbx,qword ptr [rsp+20h]  
000000013F3F144D 48 C1 E8 02          shr         rax,2  
000000013F3F1451 48 8D 34 85 01 00 00 00 lea         rsi,[rax*4+1]  
000000013F3F1459 0F 1F 80 00 00 00 00 nop         dword ptr [rax]  
   276: 		{
   277: 			q1=Kmod*ValoriPreCalcolatiSin[Indice][i][j];
000000013F3F1460 66 0F 28 D7          movapd      xmm2,xmm7  
   278: 			q2=Kmod*ValoriPreCalcolatiCos[Indice][i][j];
000000013F3F1464 66 0F 28 DF          movapd      xmm3,xmm7  
   281: 			dydx[i+n]+=q2*y[j]-q1*y[j+n];
000000013F3F1468 48 83 C5 20          add         rbp,20h  
000000013F3F146C 48 FF C8             dec         rax  
000000013F3F146F F2 41 0F 59 5C 2B D8 mulsd       xmm3,mmword ptr [r11+rbp-28h]  
000000013F3F1476 F2 0F 59 54 29 D8    mulsd       xmm2,mmword ptr [rcx+rbp-28h]  
000000013F3F147C 66 0F 28 CB          movapd      xmm1,xmm3  
000000013F3F1480 66 0F 28 C2          movapd      xmm0,xmm2  
000000013F3F1484 F2 41 0F 59 44 2C D8 mulsd       xmm0,mmword ptr [r12+rbp-28h]  
000000013F3F148B F2 0F 59 4D D8       mulsd       xmm1,mmword ptr [rbp-28h]  
000000013F3F1490 F2 0F 58 C8          addsd       xmm1,xmm0  
000000013F3F1494 F2 41 0F 10 01       movsd       xmm0,mmword ptr [r9]  
000000013F3F1499 F2 0F 5C C1          subsd       xmm0,xmm1  
000000013F3F149D F2 41 0F 11 01       movsd       mmword ptr [r9],xmm0  
000000013F3F14A2 F2 41 0F 59 5C 2C D8 mulsd       xmm3,mmword ptr [r12+rbp-28h]  
000000013F3F14A9 F2 0F 59 55 D8       mulsd       xmm2,mmword ptr [rbp-28h]  
000000013F3F14AE F2 0F 5C DA          subsd       xmm3,xmm2  
000000013F3F14B2 66 0F 28 D7          movapd      xmm2,xmm7  
000000013F3F14B6 F2 43 0F 58 1C 0A    addsd       xmm3,mmword ptr [r10+r9]  
000000013F3F14BC F2 43 0F 11 1C 0A    movsd       mmword ptr [r10+r9],xmm3  
000000013F3F14C2 66 0F 28 DF          movapd      xmm3,xmm7  
000000013F3F14C6 F2 0F 59 54 29 E0    mulsd       xmm2,mmword ptr [rcx+rbp-20h]  
000000013F3F14CC F2 41 0F 59 5C 2B E0 mulsd       xmm3,mmword ptr [r11+rbp-20h]  
000000013F3F14D3 66 0F 28 C2          movapd      xmm0,xmm2  
000000013F3F14D7 F2 42 0F 59 44 F5 E0 mulsd       xmm0,mmword ptr [rbp+r14*8-20h]  
000000013F3F14DE 66 0F 28 CB          movapd      xmm1,xmm3  
000000013F3F14E2 F2 0F 59 4D E0       mulsd       xmm1,mmword ptr [rbp-20h]  
000000013F3F14E7 F2 0F 58 C8          addsd       xmm1,xmm0  
000000013F3F14EB F2 41 0F 10 01       movsd       xmm0,mmword ptr [r9]  
000000013F3F14F0 F2 0F 5C C1          subsd       xmm0,xmm1  
000000013F3F14F4 F2 41 0F 11 01       movsd       mmword ptr [r9],xmm0  
000000013F3F14F9 F2 42 0F 59 5C F5 E0 mulsd       xmm3,mmword ptr [rbp+r14*8-20h]  
000000013F3F1500 F2 0F 59 55 E0       mulsd       xmm2,mmword ptr [rbp-20h]  
000000013F3F1505 F2 0F 5C DA          subsd       xmm3,xmm2  
000000013F3F1509 66 0F 28 D7          movapd      xmm2,xmm7  
000000013F3F150D F2 43 0F 58 1C 0A    addsd       xmm3,mmword ptr [r10+r9]  
000000013F3F1513 F2 43 0F 11 1C 0A    movsd       mmword ptr [r10+r9],xmm3  
000000013F3F1519 66 0F 28 DF          movapd      xmm3,xmm7  
000000013F3F151D F2 0F 59 54 29 E8    mulsd       xmm2,mmword ptr [rcx+rbp-18h]  
000000013F3F1523 F2 41 0F 59 5C 2B E8 mulsd       xmm3,mmword ptr [r11+rbp-18h]  
000000013F3F152A 66 0F 28 C2          movapd      xmm0,xmm2  
000000013F3F152E F2 41 0F 59 44 2C E8 mulsd       xmm0,mmword ptr [r12+rbp-18h]  
000000013F3F1535 66 0F 28 CB          movapd      xmm1,xmm3  
000000013F3F1539 F2 0F 59 4D E8       mulsd       xmm1,mmword ptr [rbp-18h]  
000000013F3F153E F2 0F 58 C8          addsd       xmm1,xmm0  
000000013F3F1542 F2 41 0F 10 01       movsd       xmm0,mmword ptr [r9]  
000000013F3F1547 F2 0F 5C C1          subsd       xmm0,xmm1  
000000013F3F154B F2 41 0F 11 01       movsd       mmword ptr [r9],xmm0  
000000013F3F1550 F2 41 0F 59 5C 2C E8 mulsd       xmm3,mmword ptr [r12+rbp-18h]  
000000013F3F1557 F2 0F 59 55 E8       mulsd       xmm2,mmword ptr [rbp-18h]  
000000013F3F155C F2 0F 5C DA          subsd       xmm3,xmm2  
000000013F3F1560 66 0F 28 D7          movapd      xmm2,xmm7  
000000013F3F1564 F2 43 0F 58 1C 0A    addsd       xmm3,mmword ptr [r10+r9]  
000000013F3F156A F2 43 0F 11 1C 0A    movsd       mmword ptr [r10+r9],xmm3  
000000013F3F1570 66 0F 28 DF          movapd      xmm3,xmm7  
000000013F3F1574 F2 41 0F 59 54 28 E0 mulsd       xmm2,mmword ptr [r8+rbp-20h]  
000000013F3F157B F2 0F 59 5C 2B E0    mulsd       xmm3,mmword ptr [rbx+rbp-20h]  
000000013F3F1581 66 0F 28 C2          movapd      xmm0,xmm2  
000000013F3F1585 F2 0F 59 44 2A E0    mulsd       xmm0,mmword ptr [rdx+rbp-20h]  
000000013F3F158B 66 0F 28 CB          movapd      xmm1,xmm3  
000000013F3F158F F2 0F 59 4D F0       mulsd       xmm1,mmword ptr [rbp-10h]  
000000013F3F1594 F2 0F 58 C8          addsd       xmm1,xmm0  
000000013F3F1598 F2 41 0F 10 01       movsd       xmm0,mmword ptr [r9]  
000000013F3F159D F2 0F 5C C1          subsd       xmm0,xmm1  
000000013F3F15A1 F2 41 0F 11 01       movsd       mmword ptr [r9],xmm0  
000000013F3F15A6 F2 0F 59 5C 2A E0    mulsd       xmm3,mmword ptr [rdx+rbp-20h]  
000000013F3F15AC F2 0F 59 55 F0       mulsd       xmm2,mmword ptr [rbp-10h]  
000000013F3F15B1 F2 0F 5C DA          subsd       xmm3,xmm2  
000000013F3F15B5 F2 43 0F 58 1C 0A    addsd       xmm3,mmword ptr [r10+r9]  
000000013F3F15BB F2 43 0F 11 1C 0A    movsd       mmword ptr [r10+r9],xmm3  
000000013F3F15C1 0F 85 99 FE FF FF    jne         derivsU+160h (13F3F1460h)  
000000013F3F15C7 48 63 1D B2 EE 03 00 movsxd      rbx,dword ptr [N (13F430480h)]  
000000013F3F15CE 48 8B AC 24 90 00 00 00 mov         rbp,qword ptr [x]  
   274: 		dydx[i+n]=0;
   275: 		for(j=1;j<=n;j++)
000000013F3F15D6 48 3B F3             cmp         rsi,rbx  
000000013F3F15D9 0F 8F A4 00 00 00    jg          derivsU+383h (13F3F1683h)  
000000013F3F15DF 4B 8D 0C 0F          lea         rcx,[r15+r9]  
000000013F3F15E3 48 8D 04 33          lea         rax,[rbx+rsi]  
000000013F3F15E7 48 8D 14 C7          lea         rdx,[rdi+rax*8]  
000000013F3F15EB 48 8B 05 1E 04 04 00 mov         rax,qword ptr [ValoriPreCalcolatiSin (13F431A10h)]  
000000013F3F15F2 4A 8B 04 E8          mov         rax,qword ptr [rax+r13*8]  
000000013F3F15F6 4C 8B 04 08          mov         r8,qword ptr [rax+rcx]  
000000013F3F15FA 48 8B 05 17 04 04 00 mov         rax,qword ptr [ValoriPreCalcolatiCos (13F431A18h)]  
000000013F3F1601 4A 8B 04 E8          mov         rax,qword ptr [rax+r13*8]  
000000013F3F1605 4D 2B C2             sub         r8,r10  
000000013F3F1608 4C 8B 1C 08          mov         r11,qword ptr [rax+rcx]  
000000013F3F160C 48 8B C3             mov         rax,rbx  
000000013F3F160F 48 8B CB             mov         rcx,rbx  
000000013F3F1612 4D 2B DA             sub         r11,r10  
000000013F3F1615 48 2B C6             sub         rax,rsi  
000000013F3F1618 4C 2B C7             sub         r8,rdi  
000000013F3F161B 4C 2B DF             sub         r11,rdi  
000000013F3F161E 48 F7 D9             neg         rcx  
000000013F3F1621 48 FF C0             inc         rax  
   276: 		{
   277: 			q1=Kmod*ValoriPreCalcolatiSin[Indice][i][j];
000000013F3F1624 66 0F 28 D7          movapd      xmm2,xmm7  
   278: 			q2=Kmod*ValoriPreCalcolatiCos[Indice][i][j];
000000013F3F1628 66 0F 28 DF          movapd      xmm3,xmm7  
000000013F3F162C 48 83 C2 08          add         rdx,8  
000000013F3F1630 48 FF C8             dec         rax  
000000013F3F1633 F2 41 0F 59 5C 13 F8 mulsd       xmm3,mmword ptr [r11+rdx-8]  
000000013F3F163A F2 41 0F 59 54 10 F8 mulsd       xmm2,mmword ptr [r8+rdx-8]  
   279: 
   280: 			dydx[i]-=q1*y[j]+q2*y[j+n];
000000013F3F1641 66 0F 28 C3          movapd      xmm0,xmm3  
000000013F3F1645 66 0F 28 CA          movapd      xmm1,xmm2  
000000013F3F1649 F2 0F 59 42 F8       mulsd       xmm0,mmword ptr [rdx-8]  
000000013F3F164E F2 0F 59 4C CA F8    mulsd       xmm1,mmword ptr [rdx+rcx*8-8]  
000000013F3F1654 F2 0F 58 C8          addsd       xmm1,xmm0  
000000013F3F1658 F2 41 0F 10 01       movsd       xmm0,mmword ptr [r9]  
000000013F3F165D F2 0F 5C C1          subsd       xmm0,xmm1  
000000013F3F1661 F2 41 0F 11 01       movsd       mmword ptr [r9],xmm0  
   281: 			dydx[i+n]+=q2*y[j]-q1*y[j+n];
000000013F3F1666 F2 0F 59 5C CA F8    mulsd       xmm3,mmword ptr [rdx+rcx*8-8]  
000000013F3F166C F2 0F 59 52 F8       mulsd       xmm2,mmword ptr [rdx-8]  
000000013F3F1671 F2 0F 5C DA          subsd       xmm3,xmm2  
000000013F3F1675 F2 43 0F 58 1C 0A    addsd       xmm3,mmword ptr [r10+r9]  
000000013F3F167B F2 43 0F 11 1C 0A    movsd       mmword ptr [r10+r9],xmm3  
000000013F3F1681 75 A1                jne         derivsU+324h (13F3F1624h)  
   270: 	double q1,q2;
   271: 
   272: 	for(i=1;i<=n;i++){
000000013F3F1683 49 83 C1 08          add         r9,8  
000000013F3F1687 48 FF CD             dec         rbp  
000000013F3F168A B8 00 00 00 00       mov         eax,0  
000000013F3F168F 48 89 AC 24 90 00 00 00 mov         qword ptr [x],rbp  
000000013F3F1697 0F 85 23 FD FF FF    jne         derivsU+0C0h (13F3F13C0h)  
000000013F3F169D 4C 8B 74 24 50       mov         r14,qword ptr [rsp+50h]  
000000013F3F16A2 4C 8B 64 24 58       mov         r12,qword ptr [rsp+58h]  
000000013F3F16A7 48 8B 74 24 60       mov         rsi,qword ptr [rsp+60h]  
000000013F3F16AC 48 8B AC 24 98 00 00 00 mov         rbp,qword ptr [y]  
   282: 		}
   283: 	}
   284: #endif
   285: }
000000013F3F16B4 0F 28 7C 24 30       movaps      xmm7,xmmword ptr [rsp+30h]  
000000013F3F16B9 48 83 C4 68          add         rsp,68h  
000000013F3F16BD 41 5F                pop         r15  
000000013F3F16BF 41 5D                pop         r13  
000000013F3F16C1 5F                   pop         rdi  
000000013F3F16C2 5B                   pop         rbx  
000000013F3F16C3 C3                   ret

Come si può notare usa un sacco di registri xmm, in particolare usa questi registri all'interno dei cicli, quindi immagino che faccia eseguire alla CPU più operazioni in simultanea per ogni thread (dovrebbero essere 2 trattandosi di dati di tipo double).

Comunque non pensavo che ci fossero tutte queste istruzioni per fare due prodotti e due somme..., vabbè che ottimizza e sfrutta SSE però...

[Unrue]

Originariamente inviato da Rossi88

Come si può notare usa un sacco di registri xmm, in particolare usa questi registri all'interno dei cicli, quindi immagino che faccia eseguire alla CPU più operazioni in simultanea per ogni thread (dovrebbero essere 2 trattandosi di dati di tipo double).

Comunque non pensavo che ci fossero tutte queste istruzioni per fare due prodotti e due somme..., vabbè che ottimizza e sfrutta SSE però...

Si, sono molte, ma considera che per ogni xmm ti risparmi 4 operazioni singole se sei in singola precisione o due in doppia. Insomma, senza le SSE, avresti almeno il doppio delle istruzioni, o il doppio delle iterazioni.

E con le prossime AVX arriveranno a 8 e 4 rispettivamente..

[Rossi88]

Originariamente inviato da Unrue

Si, sono molte, ma considera che per ogni xmm ti risparmi 4 operazioni singole se sei in singola precisione o due in doppia. Insomma, senza le SSE, avresti almeno il doppio delle istruzioni, o il doppio delle iterazioni.

E con le prossime AVX arriveranno a 8 e 4 rispettivamente..

Bè allora mi comprerò un bel Sandy Bridge o un Bulldozer da 8 core , così passando da 2 o 8 core e con le istruzioni AVX riduco il tempo di elaborazione ad 1/8 (o anche meno considerando l'aumento di IPC che c'è con due generazioni di mezzo di CPU) .

Comunque chiaramente ora il tempo di elaborazione è molto più che accettabile, sono pienamente soddisfatto dei progressi che ho fatto. Magari un giorno imparerò l'assembly per futuri programmi CPU-limited, anche se dubito che riuscierei ad ottenere un codice assembly più performante di quello generato dai compilatori odierni.

[TRF83]

Prima di precalcolarti tutti i seni e coseni, avrei provato anche ad usare il buon vecchio sviluppo di Taylor (limitandoti al terzo ordine che da risultati più che accettabili). Certo..se non hai problemi di memoria, una LUT è sempre la scelta ideale..in caso contrario c'è Taylor!

[Rossi88]

Originariamente inviato da TRF83

Prima di precalcolarti tutti i seni e coseni, avrei provato anche ad usare il buon vecchio sviluppo di Taylor (limitandoti al terzo ordine che da risultati più che accettabili). Certo..se non hai problemi di memoria, una LUT è sempre la scelta ideale..in caso contrario c'è Taylor!

Ma taylor presuppone uno sviluppo intorno ad un certo punto iniziale e man mano che ci si allontana dal punto iniziale l'approssimazione peggiora, nel mio caso non c'è un limite al valore all'argomento del seno coseno, o comunque se c'è è molto grande (anche svariate centinaia). Quindi dovrei riportare il valore in un range tramite la funzione fmod, ma in tal caso dato converrebbe forse creare un vettore con i valori tabellati nel range con la funzione sin e cos dato che il numero di valori tabellati è molto più esiguo.
O forse mi sbaglio?

Vorrei sottolineare, forse è sfuggito, che i valori precalcolati sono esattamente i valori che mi servirebbero e sono calcolati in maniera esatta tramite il semplice sin e cos. Inoltre si ottiene un notevole risparmio di tempo anche se i vettori hanno una dimensione molto ragguardevole (ho provato anche con vettori da 500MB, probabilmente ciò è dovuto alla cache da 6MB e al fatto che la funzione utilizzerà i valori della tabella in maniera sequenziale, quindi saranno caricati chessò 4-5MB della tabella nella cache e sfruttati tutti, poi altri 4-5MB e così via).


Aggiungo altre informazioni che potrebbero risultare interessanti ad altri utenti:

Ho provato la versione trial delle librerie MKL.
1) A differenza delle librerie ACML, è possibile con il compilatore della Microsoft sia utilizzare le librerie dinamiche (permesso anche dalle ACML) che quelle statiche. O forse sono io incapace, eppure nel manuale delle ACML non si fa menzione di utilizzo delle librerie statiche con il compilatore Microsoft.

2) Le librerie MKL risultano più efficienti, mi consentono di ottenere una ulteriore riduzione del tempo. Nello stesso scenario del precedente post, nel quale riuscivo ad ottenere come migliore risultato 45'', con le MKL ottengo 35'' (un'altra ragguardevole riduzione del tempo)

3) Il compilatore Microsoft secondo me fa comunque un ottimo lavoro! , considerando che deve interpretare il codice che ho scritto e cercare di sfruttare SSE nel miglior modo possibile senza sapere che alla fine il codice si tratta di prodotti vettore-matrice, mentre utilizzando le MKL ho compattato il tutto a 4 prodotti vettore-matrice.

[Torav]

Cmq anche a me è successo di dover ottimizzare in qualche modo il calcolo di seni e coseni e me la sono cavata con una look up table moooooolto piccola (100-500 valori per il seno e altrettanti per il coseno) + un'interpolazione cubica. La precisione, per quel che dovevo fare io, era ottimo (in valore assoluto la differenza tra il risultato "vero" e quello trovato in questa maniera era ~10^-6 - 10^-9)

[rеpne scasb]

■

[Tommo]

Mm analisi matematica per l'ottimizzazzione, questo thread è interessantissimo
Continuate così!

[cionci]

Originariamente inviato da rеpne scasb

Come vedi per calcolare il seno di k per k qualsiasi e' sufficiente conoscere solo il seno tra 0 e P greco mezzi.

Anche lui aveva parlato di fmod Nel tuo caso c'è una divisione FP, che è molto vicina alla fmod come tempo di esecuzione.
Ma fare una fmod è più veloce di una LUT ampia ?
L'unico caso probabilmente è in presenza di un cache miss della parte di LUT interessata, ma non sono comunque sicuro che sia più lento.