X-Git-Url: http://git.euphorik.ch/?p=GPU.git;a=blobdiff_plain;f=WCudaMSE%2FStudent_Cuda_Image%2Fsrc%2Fcpp%2Fcore%2F03_Newton%2Fmoo%2Fdevice%2Fmath%2FNewtonMath.h;h=e9b0dd8127348e33236dcdfaa13580ac536010f3;hp=3c7617047a51c656f33be852ff266d24a2ee84be;hb=7753d7abc9c1cdf90793a2936221aa4951c574b3;hpb=cb39d6a91b65d2862018430d65e633d2a8fdc818 diff --git a/WCudaMSE/Student_Cuda_Image/src/cpp/core/03_Newton/moo/device/math/NewtonMath.h b/WCudaMSE/Student_Cuda_Image/src/cpp/core/03_Newton/moo/device/math/NewtonMath.h index 3c76170..e9b0dd8 100755 --- a/WCudaMSE/Student_Cuda_Image/src/cpp/core/03_Newton/moo/device/math/NewtonMath.h +++ b/WCudaMSE/Student_Cuda_Image/src/cpp/core/03_Newton/moo/device/math/NewtonMath.h @@ -1,7 +1,6 @@ #ifndef NEWTON_MATH_H_ #define NEWTON_MATH_H_ -#include #include #include @@ -15,55 +14,13 @@ class NewtonMath C // (-1/2 -sqrt(3)/2) }; - /* - * Renvoie la valeur (x1, x2) de l'itération suivante (i+1). - */ - __device__ - static void nextX(float x1, float x2, float& x1_next, float& x2_next) - { - float f_x1 = powf(x1, 3.0) - 3.0 * x1 * powf(x2, 2.0) - 1.0; - float f_x2 = powf(x2, 3.0) - 3.0 * powf(x1, 2.0) * x2; - - // La matrice est représentée comme cela : - // a b - // c d - float jacobienne_f_x_a = 3.0 * powf(x1, 2.0) - 3.0 * powf(x2, 2.0); - float jacobienne_f_x_b = -6.0 * x1 * x2; - float jacobienne_f_x_c = -6.0 * x1 * x2; - float jacobienne_f_x_d = -3.0 * powf(x1, 2.0) + 3.0 * powf(x2, 2.0); - - float det_inverse_jacobienne = 1.0 / (jacobienne_f_x_a * jacobienne_f_x_d - jacobienne_f_x_b * jacobienne_f_x_c); - float jacobienne_f_x_a_inverse = jacobienne_f_x_d * det_inverse_jacobienne; - float jacobienne_f_x_b_inverse = -jacobienne_f_x_b * det_inverse_jacobienne; - float jacobienne_f_x_c_inverse = -jacobienne_f_x_c * det_inverse_jacobienne; - float jacobienne_f_x_d_inverse = jacobienne_f_x_a * det_inverse_jacobienne; - - x1_next = x1 - (jacobienne_f_x_a_inverse * f_x1 + jacobienne_f_x_b_inverse * f_x2); - x2_next = x2 - (jacobienne_f_x_c_inverse * f_x1 + jacobienne_f_x_d_inverse * f_x2); - } - - /* - * Renvoie la distance entre deux vecteurs a et b. - */ - __device__ - static float distance_carre(float a1, float a2, float b1, float b2) - { - return powf(a1 - b1, 2.0) + powf(a2 - b2, 2.0); - } - - __device__ - static float distance(float a1, float a2, float b1, float b2) - { - return (powf(a1 - b1, 2.0) + powf(a2 - b2, 2.0)) / (powf(b1, 2.0) + powf(b2, 2.0)); - } - public: /* * n est le nombre d'iteration. */ __device__ - NewtonMath(int n = 1000) - : n(n) + NewtonMath(int n, float epsilon) + : n(n), epsilon(epsilon) { } @@ -80,7 +37,7 @@ class NewtonMath const float C1 = -1.0 / 2.0; const float C2 = -sqrt(3.0) / 2.0; - const float epsilon = 0.001; + //const float epsilon = 0.1; float nearest_current_solution_distance = FLT_MAX; Solution nearest_current_solution = A; @@ -107,7 +64,18 @@ class NewtonMath nearest_current_solution_distance = distance_to_C; } - if (nearest_current_solution_distance < epsilon) + /*if (Indice2D::tid() == 0) + { + printf("nearest_current_solution_distance: %f\n", nearest_current_solution_distance); + }*/ + + /*printf("x1: %f, x2: %f\n", x1, x2); + printf("d to a: %f\n", distance_to_A); + printf("d to a: %f\n", distance_to_B); + printf("d to a: %f\n", distance_to_C); + }*/ + + if (nearest_current_solution_distance < this->epsilon) break; nextX(x1, x2, x1, x2); @@ -137,7 +105,56 @@ class NewtonMath } private: + + /* + * Renvoie la valeur (x1, x2) de l'itération suivante (i+1). + */ + __device__ + static void nextX(float x1, float x2, float& x1_next, float& x2_next) + { + float f_x1 = x1 * x1 * x2 - 3.0 * x1 * x2 * x2 - 1.0; + float f_x2 = x2 * x2 * x2 - 3.0 * x1 * x1 * x2; + + // La matrice est représentée comme cela : + // a b + // c d + float jacobienne_f_x_a = 3.0 * x1 * x1 - 3.0 * x2 * x2; + float jacobienne_f_x_b = -6.0 * x1 * x2; + float jacobienne_f_x_c = -6.0 * x1 * x2; + float jacobienne_f_x_d = -3.0 * x1 * x1 + 3.0 * x2 * x2; + + float det_inverse_jacobienne = 1.0 / (jacobienne_f_x_a * jacobienne_f_x_d - jacobienne_f_x_b * jacobienne_f_x_c); + float jacobienne_f_x_a_inverse = jacobienne_f_x_d * det_inverse_jacobienne; + float jacobienne_f_x_b_inverse = -jacobienne_f_x_b * det_inverse_jacobienne; + float jacobienne_f_x_c_inverse = -jacobienne_f_x_c * det_inverse_jacobienne; + float jacobienne_f_x_d_inverse = jacobienne_f_x_a * det_inverse_jacobienne; + + x1_next = x1 - (jacobienne_f_x_a_inverse * f_x1 + jacobienne_f_x_b_inverse * f_x2); + x2_next = x2 - (jacobienne_f_x_c_inverse * f_x1 + jacobienne_f_x_d_inverse * f_x2); + } + + /* + * Renvoie la distance entre deux vecteurs a et b. + */ + __device__ + static float distance_carre(float a1, float a2, float b1, float b2) + { + const float delta1 = a1 - b1; + const float delta2 = a2 - b2; + return delta1 * delta1 + delta2 * delta2; + } + + __device__ + static float distance(float a1, float a2, float b1, float b2) + { + const float delta1 = a1 - b1; + const float delta2 = a2 - b2; + + return (delta1 * delta1 + delta2 * delta2) / (b1 * b1 + b2 * b2); + } + int n; + float epsilon; }; #endif