[src] angefangen openMP einzubinden

diff --git a/src/mex_build_V.cpp b/src/mex_build_V.cpp
index 5a4f71357ac83b948461bf5001d0c2723b93d154..3ad38d1a4bedbec706f5c93a8fc3d515dbaba4b0 100644 (file)
--- 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)\r
  *\r
  * 3 : Semianalytisch wie (2) und Quadratur ueber eine Achse\r
- *             Quadratur bei : min{x[1],y1} < zeta[0] * dist(x[1],y1)\r
+ *             Quadratur bei : min{x[1],y[1]} < zeta[0] * dist(x[1],y[1])\r
  *\r
  * 4 : Semianalytisch wie (2) und Quadratur ueber ein Element\r
  *             Quadratur bei : min{dia(T1),dia(T2)} < zeta[0] * dist(T1,T2)\r
@@ -30,82 +30,31 @@
  //#include "gauss.hpp"\r
  */\r
 //#define M_EPS 10e-8\r
-#include "tbb/include/tbb/parallel_for.h"\r
-#include "tbb/include/tbb/blocked_range.h"\r
-\r
-using namespace tbb;\r
+//#include "tbb/include/tbb/parallel_for.h"\r
+//#include "tbb/include/tbb/blocked_range.h"\r
+//using namespace tbb;\r
+//#include <omp.h>\r
 \r
 #include "slpRectangle.hpp"\r
 //using namespace slpR;\r
 \r
-// Gibt vom Vektor vec die Dimension != 0 zurueck\r
-int inline dimOfVec(double* vec) {\r
-       if (vec[2] != 0)\r
-               return 2;\r
-       else if (vec[1] != 0)\r
-               return 1;\r
-       else if (vec[0] != 0)\r
-               return 0;\r
-       else {\r
-               mexErrMsgTxt("length of Site is zero");\r
-               return -1;\r
-       }\r
-}\r
-\r
-\r
-// Gibt von [0 1 2] die Fehlende Zahl zurueck\r
-int inline dimOfThird(int a, int b) {\r
-       return ((-(a + b) % 3) + 3) % 3;\r
-}\r
-\r
-void inline add(double* a, double* b) {\r
-       for (int i = 0; i < 3; ++i)\r
-               a[i] += b[i];\r
-}\r
-void inline add(double* a, double* b, double* c) {\r
-       for (int i = 0; i < 3; ++i)\r
-               a[i] = b[i] + c[i];\r
-}\r
-\r
-void inline sub(double* a, double* b) {\r
-       for (int i = 0; i < 3; ++i)\r
-               a[i] -= b[i];\r
-}\r
-void inline sub(double* a, double* b, double* c) {\r
-       for (int i = 0; i < 3; ++i)\r
-               a[i] = b[i] - c[i];\r
-}\r
-\r
-void inline set(double* a, double* b) {\r
-       for (int i = 0; i < 3; ++i)\r
-               a[i] = b[i];\r
-}\r
+void evalInte(double* A, double* C, double * E, int cm, int em, int j, int k,\r
+               double (*ctypeP)(double, double, double, double, double, double, double,\r
+                               double*),\r
+               double (*ctypeO)(double, double, double, double, double, double, double,\r
+                               double*), double *zeta) {\r
 \r
-void inline set(double* a, double b) {\r
-       for (int i = 0; i < 3; ++i)\r
-               a[i] = b;\r
-}\r
 \r
-//gibt die kleinste Ecke zurueck\r
-void inline getSCorner(double* d, double x[4][3]) {\r
-       set(d,x[0]);\r
-       for (int i = 1; i < 4; ++i) {\r
-               for (int j = 0; j < 3; ++j) {\r
-                       if (d[j] > x[i][j])\r
-                               set(d, x[i]);\r
-                       else if (d[j] < x[i][j])\r
-                               break;\r
-               }\r
-       }\r
 }\r
 \r
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
+       mexPrintf("%d",omp_get_max_threads());\r
+       omp_set_num_threads(3);\r
+       mexPrintf("%d",omp_get_num_threads());\r
 \r
-       int i, j, k; //Schleifenindizes\r
-       double tmp; //Zwischenspeicherung der Einzelnen Werte\r
+       int j, k; //Schleifenindizes\r
        int count;\r
-\r
-       double x[4][3] = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};\r
+       double tmp;\r
 \r
 //Sicherheitsabfragen zu Datengroessen\r
        if ((nrhs != 4))\r
@@ -136,28 +85,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
        int type = (int) *(mxGetPr(prhs[3]));\r
        double * zeta = mxGetPr(prhs[2]);\r
 \r
-//Initialisieren der Hilfsvektoren\r
-\r
-\r
-/*\r
-       double * x[0] = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-       double * x[1] = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-       double * x[2] = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-       double * x[3] = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-*/\r
-       double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-       double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
-       double * y0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-       double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-       double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-       double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-       double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-       double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
-       double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-       double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));\r
-\r
 // Welche Funktion soll verwendet werden\r
        double (*ctypeP)(double, double, double, double, double, double, double,\r
                        double*);\r
@@ -175,7 +102,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                ctypeP = cParO2;\r
                ctypeO = cOrtO1;\r
                break;\r
-       case 3: // Analytisch oder Qy1x2\r
+       case 3: // Analytisch oder Qy[1]x2\r
                ctypeP = cParO3;\r
                ctypeO = cOrtO1;\r
                break;\r
@@ -185,135 +112,25 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                break;\r
        }\r
 \r
-//LageInformationen\r
-       int rx, rxa, rxb, ry, rya, ryb;\r
-\r
        count = 0;\r
-//Ausrechnen\r
-       for (j = 0; j < em; ++j) {\r
-               x[0][0] = C[(int) E[j] - 1];\r
-               x[0][1] = C[cm + (int) E[j] - 1];\r
-               x[0][2] = C[2 * cm + (int) E[j] - 1];\r
-\r
-               x[1][0] = C[(int) E[em + j] - 1];\r
-               x[1][1] = C[cm + (int) E[em + j] - 1];\r
-               x[1][2] = C[2 * cm + (int) E[em + j] - 1];\r
-\r
-               x[2][0] = C[(int) E[2 * em + j] - 1];\r
-               x[2][1] = C[cm + (int) E[2 * em + j] - 1];\r
-               x[2][2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
-\r
-               x[3][0] = C[(int) E[3 * em + j] - 1];\r
-               x[3][1] = C[cm + (int) E[3 * em + j] - 1];\r
-               x[3][2] = C[2 * cm + (int) E[3 * em + j] - 1];\r
-\r
-               sub(xa, x[1], x[0]);\r
-               sub(xb, x[3], x[0]);\r
-\r
-               // Lageeigenschaften des Flaechenstuecks\r
-\r
-               rxa = dimOfVec(xa);\r
-               rxb = dimOfVec(xb);\r
-               rx = dimOfThird(rxa, rxb);\r
-\r
-               //kleinste Ecke finden und fuer \delta verwenden\r
-               getSCorner(dt, x);\r
 \r
-               for (k = 0; k <= j; ++k) {\r
-//             for (k = 0; k < em; ++k) {\r
-                       y0[0] = C[(int) E[k] - 1];\r
-                       y0[1] = C[cm + (int) E[k] - 1];\r
-                       y0[2] = C[2 * cm + (int) E[k] - 1];\r
 \r
-                       y1[0] = C[(int) E[em + k] - 1];\r
-                       y1[1] = C[cm + (int) E[em + k] - 1];\r
-                       y1[2] = C[2 * cm + (int) E[em + k] - 1];\r
-\r
-                       y2[0] = C[(int) E[2 * em + k] - 1];\r
-                       y2[1] = C[cm + (int) E[2 * em + k] - 1];\r
-                       y2[2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
-\r
-                       y3[0] = C[(int) E[3 * em + k] - 1];\r
-                       y3[1] = C[cm + (int) E[3 * em + k] - 1];\r
-                       y3[2] = C[2 * cm + (int) E[3 * em + k] - 1];\r
-\r
-                       // Seiten Vektoren aufbauen\r
-                       sub(ya, y1, y0);\r
-                       sub(yb, y3, y0);\r
-\r
-                       // Lageeigenschaften des Flaechenstuecks\r
-                       rya = dimOfVec(ya);\r
-                       ryb = dimOfVec(yb);\r
-                       ry = dimOfThird(rya, ryb);\r
-\r
-                       //kleinste Ecke finden und fuer \delta verwenden\r
-                       if (ya[rya] > 0) {\r
-                               if (yb[ryb] > 0) {\r
-                                       for (i = 0; i < 3; ++i)\r
-                                               d[i] = -dt[i] + y0[i];\r
-                               } else {\r
-                                       for (i = 0; i < 3; ++i)\r
-                                               d[i] = -dt[i] + y3[i];\r
-                               }\r
-                       } else {\r
-                               if (yb[ryb] > 0) {\r
-                                       for (i = 0; i < 3; ++i)\r
-                                               d[i] = -dt[i] + y1[i];\r
-                               } else {\r
-                                       for (i = 0; i < 3; ++i)\r
-                                               d[i] = -dt[i] + y2[i];\r
-                               }\r
-                       }\r
-//            mexPrintf("(%f,%f)(%f,%f)+(%f,%f,%f)\n",fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]),\r
-//                                                     fabs(yb[rxb]), d[rxa], d[rxb], d[rx]);\r
-                       if (rx == ry) {\r
-                               if (rxa == rya) {\r
-                                       tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]),\r
-                                                       fabs(yb[rxb]), d[rxa], d[rxb], d[rx], zeta);\r
-\r
-                               } else {\r
-                                       tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[rxa]),\r
-                                                       fabs(ya[rxb]), d[rxa], d[rxb], d[rx], zeta);\r
-                               }\r
-\r
-                       } else {\r
-\r
-                               if (rxa == rya) {\r
-                                       tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), fabs(ya[rya]),\r
-                                                       fabs(yb[ryb]), d[rxb], d[rxa], d[rx], zeta);\r
-                               } else if (rxa == ryb) {\r
-                                       tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), fabs(yb[ryb]),\r
-                                                       fabs(ya[rya]), d[rxb], d[rxa], d[rx], zeta);\r
-                               } else if (rxb == rya) {\r
-                                       tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rya]),\r
-                                                       fabs(yb[ryb]), d[rxa], d[rxb], d[rx], zeta);\r
-                               } else {\r
-                                       tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[ryb]),\r
-                                                       fabs(ya[rya]), d[rxa], d[rxb], d[rx], zeta);\r
-                               }\r
 \r
+//Ausrechnen\r
+//#pragma omp parallel private(j,k,tmp)\r
+       {\r
+//#pragma omp parallel for\r
+               for (j = 0; j < em; ++j) {\r
+                       for (k = 0; k <= j; ++k) {\r
+                               tmp = evalInt(C, E, cm, em, j, k, ctypeP, ctypeO, zeta);\r
+\r
+//#pragma omp critical\r
+                               A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
+                               if (k != j)\r
+                                       A[(j * em) + k] = 1. / (4 * M_PI) * tmp;\r
                        }\r
-                       A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
-                       if (k != j)\r
-                               A[(j * em) + k] = 1. / (4 * M_PI) * tmp;\r
 \r
                }\r
-\r
        }\r
-//     cout << endl;\r
-//Rueckgabe (eventuell zurueck schreiben)\r
-       /*      mxFree(x[0]);\r
-        mxFree(x[1]);\r
-        mxFree(x[3]);\r
-        mxFree(xa);\r
-        mxFree(xb);\r
-        mxFree(y0);\r
-        mxFree(y1);\r
-        mxFree(y3);\r
-        mxFree(ya);\r
-        mxFree(yb);\r
-        mxFree(d);\r
-        mxFree(dt);\r
-        */\r
        return;\r
 }\r

diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp
index aad61cbb4c51d9f9de503e465563255b1b64b863..c88a0c6e1905aa282c3cc79f2156b7df153f7800 100644 (file)
--- 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,
 \r
        return calcParIntA(b, d, t, v, d1, d2, d3);\r
 }\r
+\r
+// Gibt vom Vektor vec die Dimension != 0 zurueck\r
+int inline dimOfVec(double* vec) {\r
+       if (vec[2] != 0)\r
+               return 2;\r
+       else if (vec[1] != 0)\r
+               return 1;\r
+       else if (vec[0] != 0)\r
+               return 0;\r
+       else {\r
+//             mexErrMsgTxt("length of Site is zero");\r
+               return -1;\r
+       }\r
+}\r
+\r
+// Gibt von [0 1 2] die Fehlende Zahl zurueck\r
+int inline dimOfThird(int a, int b) {\r
+       return ((-(a + b) % 3) + 3) % 3;\r
+}\r
+\r
+void inline add(double* a, double* b) {\r
+       for (int i = 0; i < 3; ++i)\r
+               a[i] += b[i];\r
+}\r
+void inline add(double* a, double* b, double* c) {\r
+       for (int i = 0; i < 3; ++i)\r
+               a[i] = b[i] + c[i];\r
+}\r
+\r
+void inline sub(double* a, double* b) {\r
+       for (int i = 0; i < 3; ++i)\r
+               a[i] -= b[i];\r
+}\r
+void inline sub(double* a, double* b, double* c) {\r
+       for (int i = 0; i < 3; ++i)\r
+               a[i] = b[i] - c[i];\r
+}\r
+\r
+void inline set(double* a, double* b) {\r
+       for (int i = 0; i < 3; ++i)\r
+               a[i] = b[i];\r
+}\r
+\r
+void inline set(double* a, double b) {\r
+       for (int i = 0; i < 3; ++i)\r
+               a[i] = b;\r
+}\r
+\r
+//gibt die kleinste Ecke zurueck\r
+void inline getSCorner(double* d, double x[4][3]) {\r
+       set(d, x[0]);\r
+       for (int i = 1; i < 4; ++i) {\r
+               for (int j = 0; j < 3; ++j) {\r
+                       if (d[j] > x[i][j])\r
+                               set(d, x[i]);\r
+                       else if (d[j] < x[i][j])\r
+                               break;\r
+               }\r
+       }\r
+}\r
+\r
+double evalInt(double* C, double * E, int cm, int em, int j, int k,\r
+               double (*ctypeP)(double, double, double, double, double, double, double,\r
+                               double*),\r
+               double (*ctypeO)(double, double, double, double, double, double, double,\r
+                               double*), double *zeta) {\r
+       double x[4][3] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };\r
+       double y[4][3] = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } };\r
+       double xa[3] = { 0, 0, 0 };\r
+       double xb[3] = { 0, 0, 0 };\r
+       double ya[3] = { 0, 0, 0 };\r
+       double yb[3] = { 0, 0, 0 };\r
+       double d[3] = { 0, 0, 0 };\r
+       double dt[3] = { 0, 0, 0 };\r
+\r
+       //LageInformationen\r
+       int rx, rxa, rxb, ry, rya, ryb;\r
+\r
+       x[0][0] = C[(int) E[j] - 1];\r
+       x[0][1] = C[cm + (int) E[j] - 1];\r
+       x[0][2] = C[2 * cm + (int) E[j] - 1];\r
+\r
+       x[1][0] = C[(int) E[em + j] - 1];\r
+       x[1][1] = C[cm + (int) E[em + j] - 1];\r
+       x[1][2] = C[2 * cm + (int) E[em + j] - 1];\r
+\r
+       x[2][0] = C[(int) E[2 * em + j] - 1];\r
+       x[2][1] = C[cm + (int) E[2 * em + j] - 1];\r
+       x[2][2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
+\r
+       x[3][0] = C[(int) E[3 * em + j] - 1];\r
+       x[3][1] = C[cm + (int) E[3 * em + j] - 1];\r
+       x[3][2] = C[2 * cm + (int) E[3 * em + j] - 1];\r
+\r
+       // Seitenvektoren aufbauen\r
+       sub(xa, x[1], x[0]);\r
+       sub(xb, x[3], x[0]);\r
+\r
+       // Lageeigenschaften des Flaechenstuecks\r
+       rxa = dimOfVec(xa);\r
+       rxb = dimOfVec(xb);\r
+       rx = dimOfThird(rxa, rxb);\r
+\r
+       //kleinste Ecke finden und fuer \delta verwenden\r
+       getSCorner(dt, x);\r
+\r
+       y[0][0] = C[(int) E[k] - 1];\r
+       y[0][1] = C[cm + (int) E[k] - 1];\r
+       y[0][2] = C[2 * cm + (int) E[k] - 1];\r
+\r
+       y[1][0] = C[(int) E[em + k] - 1];\r
+       y[1][1] = C[cm + (int) E[em + k] - 1];\r
+       y[1][2] = C[2 * cm + (int) E[em + k] - 1];\r
+\r
+       y[2][0] = C[(int) E[2 * em + k] - 1];\r
+       y[2][1] = C[cm + (int) E[2 * em + k] - 1];\r
+       y[2][2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
+\r
+       y[3][0] = C[(int) E[3 * em + k] - 1];\r
+       y[3][1] = C[cm + (int) E[3 * em + k] - 1];\r
+       y[3][2] = C[2 * cm + (int) E[3 * em + k] - 1];\r
+\r
+       // Seiten Vektoren aufbauen\r
+       sub(ya, y[1], y[0]);\r
+       sub(yb, y[3], y[0]);\r
+\r
+       // Lageeigenschaften des Flaechenstuecks\r
+       rya = dimOfVec(ya);\r
+       ryb = dimOfVec(yb);\r
+       ry = dimOfThird(rya, ryb);\r
+\r
+       //kleinste Ecke finden\r
+       getSCorner(d, y);\r
+\r
+       //\delta berechnen\r
+       sub(d, dt);\r
+\r
+//            mexPrintf("(%f,%f)(%f,%f)+(%f,%f,%f)\n",fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]),\r
+//                                                     fabs(yb[rxb]), d[rxa], d[rxb], d[rx]);\r
+       if (rx == ry) {\r
+               if (rxa == rya) {\r
+                       return ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]),\r
+                                       fabs(yb[rxb]), d[rxa], d[rxb], d[rx], zeta);\r
+\r
+               } else {\r
+                       return ctypeP(fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[rxa]),\r
+                                       fabs(ya[rxb]), d[rxa], d[rxb], d[rx], zeta);\r
+               }\r
+\r
+       } else {\r
+\r
+               if (rxa == rya) {\r
+                       return ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), fabs(ya[rya]),\r
+                                       fabs(yb[ryb]), d[rxb], d[rxa], d[rx], zeta);\r
+               } else if (rxa == ryb) {\r
+                       return ctypeO(fabs(xb[rxb]), fabs(xa[rxa]), fabs(yb[ryb]),\r
+                                       fabs(ya[rya]), d[rxb], d[rxa], d[rx], zeta);\r
+               } else if (rxb == rya) {\r
+                       return ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rya]),\r
+                                       fabs(yb[ryb]), d[rxa], d[rxb], d[rx], zeta);\r
+               } else {\r
+                       return ctypeO(fabs(xa[rxa]), fabs(xb[rxb]), fabs(yb[ryb]),\r
+                                       fabs(ya[rya]), d[rxa], d[rxb], d[rx], zeta);\r
+               }\r
+\r
+       }\r
+}\r

diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp
index 1bfd66b9261331884f421482ec133c5d70cc0671..4c3a22d62ebd8965ad6986c4befe42e98e65ed70 100644 (file)
--- 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 46ed8a630ea0d14863e15b8b7ebba258509effc7..d0b629c830d85d0d4e0974300e14cac959836923 100644 (file)
--- 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