[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 (executable)
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 6f5d9dc32e9f3710f83e6c365690e640942150bc..91d6d96f476eb89d8f553dd5522131272eda329a 100644 (file)
--- a/src/mex_build_V.cpp
+++ b/src/mex_build_V.cpp
@@ -33,31 +33,91 @@
 //#include "tbb/include/tbb/parallel_for.h"\r
 //#include "tbb/include/tbb/blocked_range.h"\r
 //using namespace tbb;\r
-//#include <omp.h>\r
+#include <omp.h>\r
 \r
 #include "slpRectangle.hpp"\r
 //using namespace slpR;\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
+// 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
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
        /*\r
-       mexPrintf("%d",omp_get_max_threads());\r
-       omp_set_num_threads(3);\r
-       mexPrintf("%d",omp_get_num_threads());\r
-       */\r
+        mexPrintf("%d",omp_get_max_threads());\r
+        omp_set_num_threads(3);\r
+        mexPrintf("%d",omp_get_num_threads());\r
+        */\r
 \r
        int j, k; //Schleifenindizes\r
        int count;\r
        double tmp;\r
 \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
 //Sicherheitsabfragen zu Datengroessen\r
        if ((nrhs != 4))\r
                mexErrMsgTxt(\r
@@ -115,18 +175,107 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
        }\r
 \r
        count = 0;\r
+       omp_set_dynamic(1);\r
 \r
-\r
+       //LageInformationen\r
+       int rx, rxa, rxb, ry, rya, ryb;\r
 \r
 //Ausrechnen\r
-//#pragma omp parallel private(j,k,tmp)\r
+#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)\r
        {\r
-//#pragma omp parallel for\r
+#pragma omp parallel for\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
+                       // 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
                        for (k = 0; k <= j; ++k) {\r
-                               tmp = evalInt(C, E, cm, em, j, k, ctypeP, ctypeO, zeta);\r
 \r
-//#pragma omp critical\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
+                                               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
+                               }\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

diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp
index c88a0c6e1905aa282c3cc79f2156b7df153f7800..be3928d287ab6008e0894c9a2335fca194d8d21c 100644 (file)
--- 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);\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
+return 0;\r
 }\r

diff --git a/src/test_sol.m b/src/test_sol.m
index a741ec74aa544fee3831ed01734cdaa3c802f6ad..82836a1545621e37661c111f9393f1cda2a04a72 100644 (file)
--- 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];