From: Peter Schaefer Date: Fri, 3 Aug 2012 12:18:21 +0000 (+0200) Subject: [src] openMP passt erst mal so X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=b25bad0d30e271fe75b06a17563868fc4db59e14;p=bacc.git [src] openMP passt erst mal so --- diff --git a/src/mex_build_V.cpp b/src/mex_build_V.cpp index 00e0d56..501b983 100644 --- a/src/mex_build_V.cpp +++ b/src/mex_build_V.cpp @@ -22,7 +22,7 @@ * Peter Schaefer */ -//#define DEBUG +#define DEBUG 0 #define PARALLEL //#include @@ -216,7 +216,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { int actualNumberOfThreads = omp_get_max_threads(); - if(MAX_WORKER= 3 mexPrintf("matrixSize (em) : %d\n" "targetSize : %d\n" "actualNumberOfThreads: %d\n", em, targetSize, @@ -240,137 +240,147 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { //Ausrechnen #ifdef PARALLEL -#pragma omp parallel private(i,j,k,tmp,x,y,xa,xb,ya,yb,d,dt,rx, rxa, rxb, ry, rya, ryb,lastRow,firstRow) shared(C,E,ctypeO,ctypeP,em,targetSize,actualNumberOfThreads) +//#pragma omp parallel p { -#pragma omp for -#endif +#pragma omp parallel for schedule(static,1) private(i,j,k,tmp,x,y,xa,xb,ya,yb,d,dt,rx, rxa, rxb, ry, rya, ryb,lastRow,firstRow) shared(C,E,ctypeO,ctypeP,em,targetSize,actualNumberOfThreads) - for (i = 0; i < actualNumberOfThreads; ++i) { - firstRow = (int) sqrt(2*i*targetSize); - if (i == actualNumberOfThreads - 1) { - lastRow = em - 1; - } else { - lastRow = (int) sqrt(2*(i+1)*targetSize)-1; - } -#ifdef DEBUG - mexPrintf("thread Nr. : %d\n" - "firstRow : %d\n" - "lastRow : %d\n", i, firstRow, lastRow); #endif + for (i = 0; i < actualNumberOfThreads; ++i) { + firstRow = (int) sqrt(2 * i * targetSize); + if (i == actualNumberOfThreads - 1) { + lastRow = em - 1; + } else { + lastRow = (int) sqrt(2 * (i + 1) * targetSize) - 1; + } - for (j = firstRow; j <= lastRow; ++j) { - for (k = 0; k <= j; ++k) { - x[0][0] = C[(int) E[j] - 1]; - x[0][1] = C[cm + (int) E[j] - 1]; - x[0][2] = C[2 * cm + (int) E[j] - 1]; - - x[1][0] = C[(int) E[em + j] - 1]; - x[1][1] = C[cm + (int) E[em + j] - 1]; - x[1][2] = C[2 * cm + (int) E[em + j] - 1]; - - x[2][0] = C[(int) E[2 * em + j] - 1]; - x[2][1] = C[cm + (int) E[2 * em + j] - 1]; - x[2][2] = C[2 * cm + (int) E[2 * em + j] - 1]; - - x[3][0] = C[(int) E[3 * em + j] - 1]; - x[3][1] = C[cm + (int) E[3 * em + j] - 1]; - x[3][2] = C[2 * cm + (int) E[3 * em + j] - 1]; - - // Seitenvektoren aufbauen - sub(xa, x[1], x[0]); - sub(xb, x[3], x[0]); - - // Lageeigenschaften des Flaechenstuecks - rxa = dimOfVec(xa); - rxb = dimOfVec(xb); - rx = dimOfThird(rxa, rxb); - - //kleinste Ecke finden und fuer \delta verwenden - getSCorner(dt, x); - - y[0][0] = C[(int) E[k] - 1]; - y[0][1] = C[cm + (int) E[k] - 1]; - y[0][2] = C[2 * cm + (int) E[k] - 1]; - - y[1][0] = C[(int) E[em + k] - 1]; - y[1][1] = C[cm + (int) E[em + k] - 1]; - y[1][2] = C[2 * cm + (int) E[em + k] - 1]; - - y[2][0] = C[(int) E[2 * em + k] - 1]; - y[2][1] = C[cm + (int) E[2 * em + k] - 1]; - y[2][2] = C[2 * cm + (int) E[2 * em + k] - 1]; - - y[3][0] = C[(int) E[3 * em + k] - 1]; - y[3][1] = C[cm + (int) E[3 * em + k] - 1]; - y[3][2] = C[2 * cm + (int) E[3 * em + k] - 1]; - - // Seiten Vektoren aufbauen - sub(ya, y[1], y[0]); - sub(yb, y[3], y[0]); - - // Lageeigenschaften des Flaechenstuecks - rya = dimOfVec(ya); - ryb = dimOfVec(yb); - ry = dimOfThird(rya, ryb); - - //kleinste Ecke finden - getSCorner(d, y); - - //\delta berechnen - sub(d, dt); -#ifdef DEBUG - mexPrintf("[%d %d](%f,%f)(%f,%f)+(%f,%f,%f)\n", j, k, - fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]), - fabs(yb[rxb]), d[rxa], d[rxb], d[rx]); -#endif - if (rx == ry) { - if (rxa == rya) { - tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), - fabs(ya[rxa]), fabs(yb[rxb]), d[rxa], d[rxb], - d[rx], zeta); - - } else { - tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), - fabs(yb[rxa]), fabs(ya[rxb]), d[rxa], d[rxb], - d[rx], zeta); - } - - } else { - - if (rxa == rya) { - tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), - fabs(ya[rya]), fabs(yb[ryb]), d[rxb], d[rxa], - d[rx], zeta); - } else if (rxa == ryb) { - tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), - fabs(yb[ryb]), fabs(ya[rya]), d[rxb], d[rxa], - d[rx], zeta); - } else if (rxb == rya) { - tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), - fabs(ya[rya]), fabs(yb[ryb]), d[rxa], d[rxb], - d[rx], zeta); - } else { - tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), - fabs(yb[ryb]), fabs(ya[rya]), d[rxa], d[rxb], - d[rx], zeta); - } +#if DEBUG >= 4 +#pragma omp critical + mexPrintf("thread Nr.: (%d) | firstRow/lastRow: [%d / %d] \n", - } #ifdef PARALLEL -#pragma omp critical + omp_get_thread_num() +#else + i #endif - A[(k * em) + j] = 1. / (4 * M_PI) * tmp; - if (k != j) - A[(j * em) + k] = 1. / (4 * M_PI) * tmp; -#ifdef PARALLEL + , firstRow, lastRow); #pragma omp end critical #endif + { + for (j = firstRow; j <= lastRow; ++j) { + for (k = 0; k <= j; ++k) { + x[0][0] = C[(int) E[j] - 1]; + x[0][1] = C[cm + (int) E[j] - 1]; + x[0][2] = C[2 * cm + (int) E[j] - 1]; + + x[1][0] = C[(int) E[em + j] - 1]; + x[1][1] = C[cm + (int) E[em + j] - 1]; + x[1][2] = C[2 * cm + (int) E[em + j] - 1]; + + x[2][0] = C[(int) E[2 * em + j] - 1]; + x[2][1] = C[cm + (int) E[2 * em + j] - 1]; + x[2][2] = C[2 * cm + (int) E[2 * em + j] - 1]; + + x[3][0] = C[(int) E[3 * em + j] - 1]; + x[3][1] = C[cm + (int) E[3 * em + j] - 1]; + x[3][2] = C[2 * cm + (int) E[3 * em + j] - 1]; + + // Seitenvektoren aufbauen + sub(xa, x[1], x[0]); + sub(xb, x[3], x[0]); + + // Lageeigenschaften des Flaechenstuecks + rxa = dimOfVec(xa); + rxb = dimOfVec(xb); + rx = dimOfThird(rxa, rxb); + + //kleinste Ecke finden und fuer \delta verwenden + getSCorner(dt, x); + + y[0][0] = C[(int) E[k] - 1]; + y[0][1] = C[cm + (int) E[k] - 1]; + y[0][2] = C[2 * cm + (int) E[k] - 1]; + + y[1][0] = C[(int) E[em + k] - 1]; + y[1][1] = C[cm + (int) E[em + k] - 1]; + y[1][2] = C[2 * cm + (int) E[em + k] - 1]; + + y[2][0] = C[(int) E[2 * em + k] - 1]; + y[2][1] = C[cm + (int) E[2 * em + k] - 1]; + y[2][2] = C[2 * cm + (int) E[2 * em + k] - 1]; + + y[3][0] = C[(int) E[3 * em + k] - 1]; + y[3][1] = C[cm + (int) E[3 * em + k] - 1]; + y[3][2] = C[2 * cm + (int) E[3 * em + k] - 1]; + + // Seiten Vektoren aufbauen + sub(ya, y[1], y[0]); + sub(yb, y[3], y[0]); + + // Lageeigenschaften des Flaechenstuecks + rya = dimOfVec(ya); + ryb = dimOfVec(yb); + ry = dimOfThird(rya, ryb); + + //kleinste Ecke finden + getSCorner(d, y); + + //\delta berechnen + sub(d, dt); +#if DEBUG >= 5 +#pragma omp critical + mexPrintf("[%d %d](%f,%f)(%f,%f)+(%f,%f,%f)\n", j, k, + fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]), + fabs(yb[rxb]), d[rxa], d[rxb], d[rx]); +#pragma omp end critical +#endif + if (rx == ry) { + if (rxa == rya) { + tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), + fabs(ya[rxa]), fabs(yb[rxb]), d[rxa], + d[rxb], d[rx], zeta); + + } else { + tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), + fabs(yb[rxa]), fabs(ya[rxb]), d[rxa], + d[rxb], d[rx], zeta); + } + + } else { + + if (rxa == rya) { + tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), + fabs(ya[rya]), fabs(yb[ryb]), d[rxb], + d[rxa], d[rx], zeta); + } else if (rxa == ryb) { + tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), + fabs(yb[ryb]), fabs(ya[rya]), d[rxb], + d[rxa], d[rx], zeta); + } else if (rxb == rya) { + tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), + fabs(ya[rya]), fabs(yb[ryb]), d[rxa], + d[rxb], d[rx], zeta); + } else { + tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), + fabs(yb[ryb]), fabs(ya[rya]), d[rxa], + d[rxb], d[rx], zeta); + } + + } + +#pragma omp critical + + A[(k * em) + j] = 1. / (4 * M_PI) * tmp; + if (k != j) + A[(j * em) + k] = 1. / (4 * M_PI) * tmp; + +#pragma omp end critical + + } + } } } - } - #ifdef PARALLEL -} + } #endif return; diff --git a/src/test_sol.m b/src/test_sol.m index 2a8e973..9168d75 100644 --- a/src/test_sol.m +++ b/src/test_sol.m @@ -9,7 +9,7 @@ mex mex_build_V.cpp slpRectangle.cpp CFLAGS="\$CFLAGS -fopenmp" CXXFLAGS="\$CXXF % Test ausführen %Anzahl der Schritte -steps = 16; +steps = 2; %Art der Berechnungen type = [1];