From: Peter Schaefer Date: Thu, 2 Aug 2012 23:06:50 +0000 (+0200) Subject: [src] angefangen openMP einzubinden X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=2ee53c4a890c08a0a3303a0c590a64c0f3e72f2f;p=bacc.git [src] angefangen openMP einzubinden --- diff --git a/src/mex_build_V.cpp b/src/mex_build_V.cpp index 5a4f713..3ad38d1 100644 --- a/src/mex_build_V.cpp +++ b/src/mex_build_V.cpp @@ -14,7 +14,7 @@ * Quadratur bei : max{dia(T1),dia(T2)} < zeta[0] * dist(T1,T2) * * 3 : Semianalytisch wie (2) und Quadratur ueber eine Achse - * Quadratur bei : min{x[1],y1} < zeta[0] * dist(x[1],y1) + * Quadratur bei : min{x[1],y[1]} < zeta[0] * dist(x[1],y[1]) * * 4 : Semianalytisch wie (2) und Quadratur ueber ein Element * Quadratur bei : min{dia(T1),dia(T2)} < zeta[0] * dist(T1,T2) @@ -30,82 +30,31 @@ //#include "gauss.hpp" */ //#define M_EPS 10e-8 -#include "tbb/include/tbb/parallel_for.h" -#include "tbb/include/tbb/blocked_range.h" - -using namespace tbb; +//#include "tbb/include/tbb/parallel_for.h" +//#include "tbb/include/tbb/blocked_range.h" +//using namespace tbb; +//#include #include "slpRectangle.hpp" //using namespace slpR; -// 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 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) { -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()); - int i, j, k; //Schleifenindizes - double tmp; //Zwischenspeicherung der Einzelnen Werte + int j, k; //Schleifenindizes int count; - - double x[4][3] = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}; + double tmp; //Sicherheitsabfragen zu Datengroessen if ((nrhs != 4)) @@ -136,28 +85,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { int type = (int) *(mxGetPr(prhs[3])); double * zeta = mxGetPr(prhs[2]); -//Initialisieren der Hilfsvektoren - - -/* - double * x[0] = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - double * x[1] = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - double * x[2] = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - double * x[3] = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); -*/ - double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - - double * y0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - - double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL)); - // Welche Funktion soll verwendet werden double (*ctypeP)(double, double, double, double, double, double, double, double*); @@ -175,7 +102,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { ctypeP = cParO2; ctypeO = cOrtO1; break; - case 3: // Analytisch oder Qy1x2 + case 3: // Analytisch oder Qy[1]x2 ctypeP = cParO3; ctypeO = cOrtO1; break; @@ -185,135 +112,25 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { break; } -//LageInformationen - int rx, rxa, rxb, ry, rya, ryb; - count = 0; -//Ausrechnen - 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]; - - 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) { -// for (k = 0; k < em; ++k) { - y0[0] = C[(int) E[k] - 1]; - y0[1] = C[cm + (int) E[k] - 1]; - y0[2] = C[2 * cm + (int) E[k] - 1]; - y1[0] = C[(int) E[em + k] - 1]; - y1[1] = C[cm + (int) E[em + k] - 1]; - y1[2] = C[2 * cm + (int) E[em + k] - 1]; - - y2[0] = C[(int) E[2 * em + k] - 1]; - y2[1] = C[cm + (int) E[2 * em + k] - 1]; - y2[2] = C[2 * cm + (int) E[2 * em + k] - 1]; - - y3[0] = C[(int) E[3 * em + k] - 1]; - y3[1] = C[cm + (int) E[3 * em + k] - 1]; - y3[2] = C[2 * cm + (int) E[3 * em + k] - 1]; - - // Seiten Vektoren aufbauen - sub(ya, y1, y0); - sub(yb, y3, y0); - - // Lageeigenschaften des Flaechenstuecks - rya = dimOfVec(ya); - ryb = dimOfVec(yb); - ry = dimOfThird(rya, ryb); - - //kleinste Ecke finden und fuer \delta verwenden - if (ya[rya] > 0) { - if (yb[ryb] > 0) { - for (i = 0; i < 3; ++i) - d[i] = -dt[i] + y0[i]; - } else { - for (i = 0; i < 3; ++i) - d[i] = -dt[i] + y3[i]; - } - } else { - if (yb[ryb] > 0) { - for (i = 0; i < 3; ++i) - d[i] = -dt[i] + y1[i]; - } else { - for (i = 0; i < 3; ++i) - d[i] = -dt[i] + y2[i]; - } - } -// 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); - } +//Ausrechnen +//#pragma omp parallel private(j,k,tmp) + { +//#pragma omp parallel for + for (j = 0; j < em; ++j) { + for (k = 0; k <= j; ++k) { + tmp = evalInt(C, E, cm, em, j, k, ctypeP, ctypeO, 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; } - A[(k * em) + j] = 1. / (4 * M_PI) * tmp; - if (k != j) - A[(j * em) + k] = 1. / (4 * M_PI) * tmp; } - } -// cout << endl; -//Rueckgabe (eventuell zurueck schreiben) - /* mxFree(x[0]); - mxFree(x[1]); - mxFree(x[3]); - mxFree(xa); - mxFree(xb); - mxFree(y0); - mxFree(y1); - mxFree(y3); - mxFree(ya); - mxFree(yb); - mxFree(d); - mxFree(dt); - */ return; } diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp index aad61cb..c88a0c6 100644 --- a/src/slpRectangle.cpp +++ b/src/slpRectangle.cpp @@ -566,3 +566,170 @@ 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); + } + + } +} diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp index 1bfd66b..4c3a22d 100644 --- a/src/slpRectangle.hpp +++ b/src/slpRectangle.hpp @@ -50,5 +50,11 @@ double cParO3(double, double, double, double, double, double, double, double*); // A oder Q oder Qx1x2 double cParO4(double, double, double, double, double, double, double, double*); +double evalInt(double *, double *, int, int, int, int, + double (*)(double, double, double, double, double, double, double, + double*), + double (*)(double, double, double, double, double, double, double, + double*), double *); + //} #endif diff --git a/src/test_sol.m b/src/test_sol.m index 46ed8a6..d0b629c 100644 --- a/src/test_sol.m +++ b/src/test_sol.m @@ -2,6 +2,10 @@ 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 mex_build_V.cpp slpRectangle.cpp % Test ausführen