mex Funktion beschleunigen (Bottleneck unklar)

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

Vladi

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 03.12.2015, 12:41 Titel: mex Funktion beschleunigen (Bottleneck unklar)

Hallo Forum,

ich habe ein Problem mit einer mex Funktion, die ich geschrieben habe. Sie berechnet den Wert eines vierfach Integrals. Da das in Matlab eher lange dauert, habe ich das als mex Funktion (Test.cpp) geschrieben.

Code:

#include <mex.h>
#include<time.h>

float getTime() {
clock_t start = clock();
return (float)start / CLOCKS_PER_SEC;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs,
const mxArray *prhs[]) {

/* For matrix A (m x n) the element A[a, b] is accessed
* lineraly in Matlab by A[(b - 1) * m + a]. Since indices
* in C++ start with zero, the same element corresponds
* to A[b * m + a].
*
* delta : double scalar
* path : (1 x nSteps + 1) double vector
* mat_A : (nSteps x nSteps + 1) double matrix
* mat_B : (nSteps x nSteps + 1) double matrix
*/
double delta = *mxGetPr(prhs[0]);
double *path = mxGetPr(prhs[1]);
double *mat_A = mxGetPr(prhs[2]);
double *mat_B = mxGetPr(prhs[3]);

int nSteps = mxGetM(prhs[2]);

float t = 0.0;
float temp = 0.0;
float t1 = 0.0;
float temp1 = 0.0;

double I_ds = 0.0;
for(int s = 0; s < nSteps; s++){
double I_dt = 0.0;

temp1 = getTime();
for(int t = s; t < nSteps; t++){

temp = getTime();
double I_dv = 0.0;
for(int v = t; v < nSteps; v++){
double I_dWr = 0.0;
for(int r = t; r < v + 1; r++){
I_dWr += mat_A[r * nSteps + s];
}
I_dv += mat_B[v * nSteps + t] * I_dWr;
}
I_dv *= delta;
I_dt += path[t] * I_dv;
t += getTime() - temp;

}
t1 += getTime() - temp1;

I_ds += path[s] * I_dt * delta;
}
I_ds *= delta;

printf("\n\ntime %f2s\n", t);
printf("time1 %f2s\n\n", t1);

/* define output */
plhs[0] = mxCreateDoubleScalar(I_ds);

}

Funktion ohne Link?

Nach dem Kompilieren, kann ich die Funktion mittels

Code:

nSteps = 500;
delta = 0.01;
path = 0.8 * ones(1, nSteps + 1);
mat_A = 0.25 * ones(nSteps, nSteps + 1);
mat_B = 0.18 * ones(nSteps, nSteps + 1);
Test(delta, path, mat_A, mat_B)

Funktion ohne Link?

starten.
Das komische ist nun die Laufzeit. 'time' ist verschwinded klein, 'time1' ist ca 6 Sekunden und der ganz Call von Test() auch in etwa 6 Sekunden dauert.

Warum liefert die innere Messung kaum Zeit, während die äußere erheblich länger dauert?
Hat jemand eine Idee, wie man das beschleunigen kann?

Vielen Dank und viele Grüße
Vladi

Harald

Forum-Meister


	Beiträge: 24.501

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 03.12.2015, 12:45 Titel:

Hallo,

meine Frage wäre erst mal: was hast du innerhalb von MATLAB versucht?

So wie ich das sehe, sollte sich ein äquivalenter MATLAB-Code durch Vektorisierung erheblich beschleunigen lassen.

Grüße,
Harald

Vladi

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 03.12.2015, 13:11 Titel:

In der Tat hatte ich die Schleifen auch bereits in Matlab geschrieben.

Das Problem ist, dass man (meiner Meinung nach) nicht gut vektorisieren kann. In der inneren r-Schleife wird (Matlab notation)

Code:

sum(mat_A(s, t : v))

Funktion ohne Link?

, dann, in der v-Schleife brauchen wir diese Summeund multiplizieren sie, etc, etc...

Meiner Meinung nach ist das so verschachtelt, dass man nicht vektorisieren kann. Natürlich bin ich aber gerne bereit, mich vom Gegenteil überzeigen zu lassen Smile

.

Viele Grüße
Vladi

Vladi

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 03.12.2015, 13:15 Titel:

Ach übrigens sind die konstanten Matritzen natürlich nur zum testn. Im eigentlichen Programm sind die Zahlen komplett unterschiedlich! Sorry dafür!!!

Harald

Forum-Meister


	Beiträge: 24.501

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 03.12.2015, 14:06 Titel:

Hallo,

kannst du deinen bisherigen MATLAB-Code zur Verfügung stellen? Vielleicht geht da noch mehr.

Grüße,
Harald

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 03.12.2015, 14:23 Titel: Re: mex Funktion beschleunigen (Bottleneck unklar)

Hallo Vladi,

Wenn Du mat_A und mat_B transponierst, ersparst Du dem Rechner eine Menge Multiplikationen und er kann auf benachbarte Speicher-Bereiche zugreifen, die in Blöcken in den Prozessor-Cache geladen werden.

Code:

for(int v = t; v < nSteps; v++){
double I_dWr = 0.0;
for(int r = t; r < v + 1; r++){
I_dWr += mat_A[r * nSteps + s];
}
I_dv += mat_B[v * nSteps + t] * I_dWr;
}

Funktion ohne Link?

Gruß, Jan

Tim

Forum-Century


	Beiträge: 140

	Anmeldedatum: 03.11.07

	Wohnort: Stuttgart

	Version: 2011b+aktuellstes Release

Verfasst am: 03.12.2015, 14:37 Titel:

Mit OpenMP kannst du bestimmt auch noch was rausholen. Zumindest wenn du nicht schon am Bandbreitenlimit hängst bzw. dir durch Parallelisierung deine Cachelines zerschiesst.

Vladi

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 04.12.2015, 11:10 Titel:

Hallo zusammen.

Schonmal vielen Dank für die verschiedenen Vorschläge.

Hier ist noch der zur mex äquivalente Matlab code.

Code:

tic
I_ds = 0;
for s = 1 : nSteps
I_dt = 0;
for t = s : nSteps
I_dv = 0;
for v = t : nSteps
I_dWr = 0;
for r = t : v
I_dWr = I_dWr + mat_A(s, r);
end
I_dv = I_dv + mat_B(t, v) * I_dWr;
end
I_dv = I_dv * delta;
I_dt = I_dt + path(t) * I_dv;
end
I_dt = I_dt * delta;
I_ds = I_ds + path(s) * I_dt;
end
I_ds = I_ds * delta;
toc

Funktion ohne Link?

Die Laufzeitunterschiede werden (logisch) mit größerem nSteps extremer.
Wie ich bereits geschrieben habe, denke ich, dass vektorisieren schwierig ist, bin aber bereits meine Meinung zu ändern Smile

Vielen Dank für eure Kommentare
Vladi

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 04.12.2015, 12:11 Titel:

Hallo Vladi,

Der Code lässt sich hervorragend vektorisieren.
Fange bei den inneren Schleifen an:

Code:

I_dv = 0;
for v = t : nSteps
I_dWr = 0;
for r = t : v
I_dWr = I_dWr + mat_A(s, r);
end
I_dv = I_dv + mat_B(t, v) * I_dWr;
end

% ==>
I_dv = 0;
for v = t : nSteps
I_dWr = sum(mat_A(s, t:v));
I_dv = I_dv + mat_B(t, v) * I_dWr;
end

Funktion ohne Link?

Und auch dies lässt sich noch verbessern, denn die Summe über mat_A enthält ja in jeder Iteration nur einen zusätzlichen Summanden.

Code:

I_dWr = 0;
I_dv = 0;
for v = t : nSteps
I_dWr = I_dWr + mat_A(s, v);
I_dv = I_dv + mat_B(t, v) * I_dWr;
end

Funktion ohne Link?

Und so weit ich sehen kann (ich habe kein Matlab zur Zeit um es auszuprobieren) ist das wiederum:

Code:

I_dv = mat_B(t, t:nSteps) * cumsum(mat_A(s, t:nSteps)).';

Funktion ohne Link?

Damit hätte man die beiden inneren Schleifen vektorisiert und das sollte sich in der Laufzeit bemerkbar machen.

Auch hier tritt der gleiche Effekt auf wie in einer Mex-Funktion: Der Zugriff auf das langsame RAM bremst stark, wenn die beiden Matrizen nicht mehr in den Prozessor-Cache passen. Ein Transponieren der Matrizen würde es erlauben, dass auf benachbarte Speicherzellen zugegriffen wird, was im Allgemeinen die Bearbeitungsgeschwindigkeit deutlich steigert.

Gruß, Jan

Vladi

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 04.12.2015, 14:35 Titel:

Hallo Jan,

deine Version funktioniert super! Danke hierfür. Stimmt, an cumsum hatte ich gar nicht gedacht Smile

Auch werde ich das mit dem Transponieren probieren. Vielleicht bringt das auch noch etwas...

An all die Matlab-Profis: Geht vielleicht noch mehr? Kan noch mehr zusammengefasst werden?

Vielen Dank schonmal an alle!

Vladi

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 04.12.2015, 15:08 Titel:

Ok, getestet: Transponieren führt nochmal zu einer Verbesserung. Der aktuelle Code (nur Matlab) sieht jetzt (dank Jan S) so aus:

Code:

nSteps = 600;
delta = 0.01;
path = 0.8 * ones(1, nSteps + 1);
mat_A = 0.25 * ones(nSteps, nSteps + 1);
mat_B = ones(nSteps, nSteps + 1);

% transpose mat_A and mat_B
mat_A_trans = mat_A';
mat_B_trans = mat_B';

tic
I_ds = 0;
for s = 1 : nSteps
I_dt = 0;
for t = s : nSteps
I_dv = mat_B_trans(t : nSteps, t)' * cumsum(mat_A_trans(t : nSteps, s)) * delta;
I_dt = I_dt + path(t) * I_dv;
end
I_dt = I_dt * delta;
I_ds = I_ds + path(s) * I_dt;
end
I_ds = I_ds * delta
toc

Funktion ohne Link?

Wie gesagt, vielleict geht noch mehr...

Viele Grüße
Vladi

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 05.12.2015, 17:09 Titel:

Hallo Vladi,

Ich habe ein wenig dran gebastelt unhd keine nennenswerten Verbeserungen mehr gefunden. Immerhin:

Matlab 2015b, 64, Win7, Core2Duo:
Original: 96.0 sec
Verbessert: 1.72 sec

Interessant ist der Unterschied zu 2011b:
Original: 127.3 sec
Verbessert: 1.40 sec

Ich mache noch ein paar Experimente mit Mex-files.

Gruß, Jan

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 05.12.2015, 18:30 Titel:

Hallo Vladi,

Die Multiplikationen mit delta kann man noch aus den Schleifen heraus ziehen, was aber kaum etwas bringt:

Code:

mat_A_T = mat_A';
mat_B_T = mat_B';

I_ds = 0;
for s = 1 : nSteps
I_dt = 0;
for t = s:nSteps
I_dv = mat_B_T(t:nSteps, t)' * cumsum(mat_A_T(t:nSteps, s));
I_dt = I_dt + path(t) * I_dv;
end
I_ds = I_ds + path(s) * I_dt;
end
I_ds = I_ds * delta^3;

Funktion ohne Link?

Aber nun ist der Code schon so zusammen geschrumpft, dass der Profiler den Flaschenhals gut identifizieren kann. Und nun ist ein C-Mex-File sinnvoll:

Code:

% I_dv = mat_B_T(t:nSteps, t)' * cumsum(mat_A_T(t:nSteps, s));
% As C-Mex file:
I_dv = BtimesCumsumA(mat_B_T, mat_A_T, s, t, nSteps);

Funktion ohne Link?

Und:

Code:

#include "mex.h"
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
// BtimesCumsumA(mat_B_T, mat_A_T, s, t, nSteps)
double *A, *B, c, r;
mwIndex s, t, n;
mwSize sizeA1, sizeB1, count;

sizeB1 = mxGetM(prhs[0]);
sizeA1 = mxGetM(prhs[1]);

s = (mwIndex) mxGetScalar(prhs[2]);
t = (mwIndex) mxGetScalar(prhs[3]);
n = (mwIndex) mxGetScalar(prhs[4]);

// I_dv = mat_B_T(t:nSteps, t)' * cumsum(mat_A_T(t:nSteps, s));
B = mxGetPr(prhs[0]) + (t - 1) * sizeB1 + t - 1; // EDITED, Bug fix: += -> =
A = mxGetPr(prhs[1]) + (s - 1) * sizeA1 + t - 1; // EDITED, Bug fix: += -> =

count = n - t + 1;
r = 0.0;
c = 0.0;
while (count-- != 0) {
c += *A++;
r += *B++ * c;
}

plhs[0] = mxCreateDoubleScalar(r);

return;
}

Funktion ohne Link?

Damit sinkt die Laufzeit von 1.72 auf 0.55 Sekunden.
Die eigentlichen Berechnungen in der innersten Schleife sind nun sehr elementar. Das ist einer der seltenen Fälle, in denen der C-Code hübscher aussieht als der Matlab Code! :-)

Gruß, Jan

Zuletzt bearbeitet von Jan S am 13.12.2015, 01:28, insgesamt einmal bearbeitet

Vladi

Gast


	Beiträge: ---

	Anmeldedatum: ---

	Wohnort: ---

	Version: ---

Verfasst am: 06.12.2015, 01:09 Titel:

Hallo Jan,

Respekt!

Vielen Dank für die gute Analyse und deine Mühe, auch besonders für die mex Datei. Ich werde das testen sobald ich Matlab wieder griffbereit habe. Nach deinen Laufzeiten zu schließen bin ich bereits sehr optimistisch!

Das beste an der Sache... Insgesamt habe ich 11 von diesen 4-fach Integralen, so dass die Laufzeit des Codes erheblich reduziert wird.

Nochmals vielen Dank und viele Grüße
Vladi

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 06.12.2015, 15:36 Titel:

Hallo Vladi,

Gerne!
Das war ein sehr schönes Beispiel für Code-Optimierung:
1. Den Matlab-Code so weit entschlacken, wie möglich.
2. Dann wird der Bottleneck klar.
3. Nur diesen in ein C-Mex auslagern.

Wenn die hauptsächliche Arbeit nicht von einer BLAS/LAPACK-Funktion ausgeführt wird, ist der C-Code oft schneller.
Der Zugriff auf benachbarte Speicher-Zellen ist effizienter als durch das RAM zu springen.

Wenn Du nun noch die äußeren Schleifen dur PARFOR-Loops ersetzt, sollte der Rechner sein bestes gegen.

Anmerkungen: Da einem C-Mex-Code um die Ohren fliegt, wenn die Inputs falsch sind, lohnt es sich nich ein paar Micro-Sekunden in die Prüfung der Inputs zu investieren: Sitmmen Typ und Dimension?
"path" ist eine wichtige Matlab-Funktion. Den Namen für eine Variable zu verwenden, kann sehr verwirrende Auswirkungen haben, z.B. während des Debuggens.

Gruß, Jan


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.