From: Peter Schaefer Date: Fri, 3 Aug 2012 06:33:41 +0000 (+0200) Subject: [src] evalInt wieder entfernt X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=a873ec23f35063f16ad16fbb99b2c28fdaa7120c;p=bacc.git [src] evalInt wieder entfernt [src] openMP läuft ... ist aber leider langsamer --- diff --git a/res/OpenMP3.1-CCard.pdf b/res/OpenMP3.1-CCard.pdf new file mode 100755 index 0000000..5921533 Binary files /dev/null and b/res/OpenMP3.1-CCard.pdf differ diff --git a/src/mex_build_V.cpp b/src/mex_build_V.cpp index 6f5d9dc..91d6d96 100644 --- a/src/mex_build_V.cpp +++ b/src/mex_build_V.cpp @@ -33,31 +33,91 @@ //#include "tbb/include/tbb/parallel_for.h" //#include "tbb/include/tbb/blocked_range.h" //using namespace tbb; -//#include +#include #include "slpRectangle.hpp" //using namespace slpR; -void evalInte(double* A, double* C, double * E, int cm, int em, int j, int k, - double (*ctypeP)(double, double, double, double, double, double, double, - double*), - double (*ctypeO)(double, double, double, double, double, double, double, - double*), double *zeta) { +// Gibt vom Vektor vec die Dimension != 0 zurueck +int inline dimOfVec(double* vec) { + if (vec[2] != 0) + return 2; + else if (vec[1] != 0) + return 1; + else if (vec[0] != 0) + return 0; + else { +// mexErrMsgTxt("length of Site is zero"); + return -1; + } +} + +// Gibt von [0 1 2] die Fehlende Zahl zurueck +int inline dimOfThird(int a, int b) { + return ((-(a + b) % 3) + 3) % 3; +} + +void inline add(double* a, double* b) { + for (int i = 0; i < 3; ++i) + a[i] += b[i]; +} +void inline add(double* a, double* b, double* c) { + for (int i = 0; i < 3; ++i) + a[i] = b[i] + c[i]; +} +void inline sub(double* a, double* b) { + for (int i = 0; i < 3; ++i) + a[i] -= b[i]; +} +void inline sub(double* a, double* b, double* c) { + for (int i = 0; i < 3; ++i) + a[i] = b[i] - c[i]; +} +void inline set(double* a, double* b) { + for (int i = 0; i < 3; ++i) + a[i] = b[i]; +} + +void inline set(double* a, double b) { + for (int i = 0; i < 3; ++i) + a[i] = b; +} + +//gibt die kleinste Ecke zurueck +void inline getSCorner(double* d, double x[4][3]) { + set(d, x[0]); + for (int i = 1; i < 4; ++i) { + for (int j = 0; j < 3; ++j) { + if (d[j] > x[i][j]) + set(d, x[i]); + else if (d[j] < x[i][j]) + break; + } + } } void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* - mexPrintf("%d",omp_get_max_threads()); - omp_set_num_threads(3); - mexPrintf("%d",omp_get_num_threads()); - */ + mexPrintf("%d",omp_get_max_threads()); + omp_set_num_threads(3); + mexPrintf("%d",omp_get_num_threads()); + */ int j, k; //Schleifenindizes int count; double tmp; + double x[4][3] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; + double y[4][3] = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } }; + double xa[3] = { 0, 0, 0 }; + double xb[3] = { 0, 0, 0 }; + double ya[3] = { 0, 0, 0 }; + double yb[3] = { 0, 0, 0 }; + double d[3] = { 0, 0, 0 }; + double dt[3] = { 0, 0, 0 }; + //Sicherheitsabfragen zu Datengroessen if ((nrhs != 4)) mexErrMsgTxt( @@ -115,18 +175,107 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { } count = 0; + omp_set_dynamic(1); - + //LageInformationen + int rx, rxa, rxb, ry, rya, ryb; //Ausrechnen -//#pragma omp parallel private(j,k,tmp) +#pragma omp parallel private(j,k,tmp,x,y,xa,xb,ya,yb,d,dt,rx, rxa, rxb, ry, rya, ryb) shared(C,E,ctypeO,ctypeP) { -//#pragma omp parallel for +#pragma omp parallel for for (j = 0; j < em; ++j) { + 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); + for (k = 0; k <= j; ++k) { - tmp = evalInt(C, E, cm, em, j, k, ctypeP, ctypeO, zeta); -//#pragma omp critical + 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); + + // mexPrintf("(%f,%f)(%f,%f)+(%f,%f,%f)\n",fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]), + // fabs(yb[rxb]), d[rxa], d[rxb], d[rx]); + 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); + } + + } + A[(k * em) + j] = 1. / (4 * M_PI) * tmp; if (k != j) A[(j * em) + k] = 1. / (4 * M_PI) * tmp; diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp index c88a0c6..be3928d 100644 --- a/src/slpRectangle.cpp +++ b/src/slpRectangle.cpp @@ -567,169 +567,13 @@ double cParO4(double b, double d, double t, double v, double d1, double d2, return calcParIntA(b, d, t, v, d1, d2, d3); } -// Gibt vom Vektor vec die Dimension != 0 zurueck -int inline dimOfVec(double* vec) { - if (vec[2] != 0) - return 2; - else if (vec[1] != 0) - return 1; - else if (vec[0] != 0) - return 0; - else { -// mexErrMsgTxt("length of Site is zero"); - return -1; - } -} - -// Gibt von [0 1 2] die Fehlende Zahl zurueck -int inline dimOfThird(int a, int b) { - return ((-(a + b) % 3) + 3) % 3; -} - -void inline add(double* a, double* b) { - for (int i = 0; i < 3; ++i) - a[i] += b[i]; -} -void inline add(double* a, double* b, double* c) { - for (int i = 0; i < 3; ++i) - a[i] = b[i] + c[i]; -} - -void inline sub(double* a, double* b) { - for (int i = 0; i < 3; ++i) - a[i] -= b[i]; -} -void inline sub(double* a, double* b, double* c) { - for (int i = 0; i < 3; ++i) - a[i] = b[i] - c[i]; -} - -void inline set(double* a, double* b) { - for (int i = 0; i < 3; ++i) - a[i] = b[i]; -} - -void inline set(double* a, double b) { - for (int i = 0; i < 3; ++i) - a[i] = b; -} - -//gibt die kleinste Ecke zurueck -void inline getSCorner(double* d, double x[4][3]) { - set(d, x[0]); - for (int i = 1; i < 4; ++i) { - for (int j = 0; j < 3; ++j) { - if (d[j] > x[i][j]) - set(d, x[i]); - else if (d[j] < x[i][j]) - break; - } - } -} - double evalInt(double* C, double * E, int cm, int em, int j, int k, double (*ctypeP)(double, double, double, double, double, double, double, double*), double (*ctypeO)(double, double, double, double, double, double, double, double*), double *zeta) { - double x[4][3] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; - double y[4][3] = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } }; - double xa[3] = { 0, 0, 0 }; - double xb[3] = { 0, 0, 0 }; - double ya[3] = { 0, 0, 0 }; - double yb[3] = { 0, 0, 0 }; - double d[3] = { 0, 0, 0 }; - double dt[3] = { 0, 0, 0 }; - - //LageInformationen - int rx, rxa, rxb, ry, rya, ryb; - - 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); - -// mexPrintf("(%f,%f)(%f,%f)+(%f,%f,%f)\n",fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]), -// fabs(yb[rxb]), d[rxa], d[rxb], d[rx]); - if (rx == ry) { - if (rxa == rya) { - return ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]), - fabs(yb[rxb]), d[rxa], d[rxb], d[rx], zeta); - } else { - return ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[rxa]), - fabs(ya[rxb]), d[rxa], d[rxb], d[rx], zeta); - } - - } else { - if (rxa == rya) { - return ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), fabs(ya[rya]), - fabs(yb[ryb]), d[rxb], d[rxa], d[rx], zeta); - } else if (rxa == ryb) { - return ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), fabs(yb[ryb]), - fabs(ya[rya]), d[rxb], d[rxa], d[rx], zeta); - } else if (rxb == rya) { - return ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rya]), - fabs(yb[ryb]), d[rxa], d[rxb], d[rx], zeta); - } else { - return ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[ryb]), - fabs(ya[rya]), d[rxa], d[rxb], d[rx], zeta); - } - } +return 0; } diff --git a/src/test_sol.m b/src/test_sol.m index a741ec7..82836a1 100644 --- a/src/test_sol.m +++ b/src/test_sol.m @@ -2,11 +2,9 @@ format longG % Matrix MEX Funktion neu Compilieren -% mex mex_build_V.cpp slpRectangle.cpp CFLAGS="\$CFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp" -% mex mex_build_V.cpp slpRectangle.cpp COMPFLAGS="/openmp \$COMPFLAGS" -% mex COPTIMFLAGS="\$COPTIMFLAGS -fopenmp" LDOPTIMFLAGS="\$LDOPTIMFLAGS -fopenmp" -lgomp mex_build_V.cpp slpRectangle.cpp +mex -v mex_build_V.cpp slpRectangle.cpp CFLAGS="\$CFLAGS -fopenmp" CXXFLAGS="\$CXXFLAGS -fopenmp" LDFLAGS="\$LDFLAGS -fopenmp" -mex mex_build_V.cpp slpRectangle.cpp +% mex mex_build_V.cpp slpRectangle.cpp % Test ausführen @@ -14,7 +12,7 @@ mex mex_build_V.cpp slpRectangle.cpp steps = 15; %Art der Berechnungen -type = [1 2]; +type = [1]; zeta = [1];