Description: upstream: parallelization: surface
 Attempt to parallelize surface drawing (surface without cuts for now);
 implemented on behalf of imaginary.org by Jerome Benoit <calculus@rezozer.net>.
Origin: imaginary.org
Forwarded: submitted
Author:
 Jerome Benoit <calculus@rezozer.net>
Last-Update: 2016-08-02

--- a/yaccsrc/Script.cc
+++ b/yaccsrc/Script.cc
@@ -35,6 +35,13 @@
  **************************************************************************/
 
 
+#define ENABLE_PARALLELIZATION_SURFACE
+
+#ifdef ENABLE_PARALLELIZATION_SURFACE
+#define ENABLE_PARALLELIZATION_SURFACE_PTHREAD
+#endif
+
+
 #include <assert.h>
 #include <stdio.h>
 #include <errno.h>
@@ -44,6 +51,10 @@
 #include <sys/stat.h>
 #include <unistd.h>
 
+#ifdef ENABLE_PARALLELIZATION_SURFACE_PTHREAD
+#include <pthread.h>
+#endif
+
 #include "FileWriter.h"
 #include "TreePolynom.h"
 #include "Misc.h"
@@ -80,6 +91,12 @@
 // #define DEBUG
 #include "debug.h"
 
+#ifdef ENABLE_PARALLELIZATION_SURFACE
+#	ifndef SURFALGGEO_ENVIRONMENT_NPROCESSORS
+#	define SURFALGGEO_ENVIRONMENT_NPROCESSORS "SURFALGGEO_NPROCESSORS"
+#	endif
+#endif
+
 using namespace std;
 
 TSDrawingArea	*Script::display	= 0;
@@ -104,7 +121,7 @@
 	} else {
 		DMESS("adding new command " << name);
 	}
-	
+
 	symtab_add_surface_name (name, SYM_COMMAND, 1, (void*)func);
 }
 
@@ -396,21 +413,204 @@
 	checkVariables();
 
 	BEGIN("Script::setSize");
-	if(buffer->getWidth() != main_width_data ||
-	   buffer->getHeight() != main_height_data) {
+	if(buffer->getWidth() != main_width_data || buffer->getHeight() != main_height_data) {
 		buffer->Realloc(main_width_data, main_height_data);
 		if (display) {
 			display->setSize(main_width_data, main_height_data);
 		}
-		buffer->Realloc(main_width_data, main_height_data);
+		//buffer->Realloc(main_width_data, main_height_data);
 	}
 }
 
 
 extern double Y_AXIS_LR_ROTATE;
 
+#ifdef ENABLE_PARALLELIZATION_SURFACE
+static size_t Script_Tile_sys_numberof_processors(void) {
+	char * env_nprocessors=getenv(SURFALGGEO_ENVIRONMENT_NPROCESSORS);
+	size_t numberof_processor=1;
+
+	if (env_nprocessors) {
+		long int long_numberof_processor=strtol(env_nprocessors,NULL,0);
+		if (long_numberof_processor < 0) {
+			long_numberof_processor=-long_numberof_processor;
+		}
+		else {
+			const long int long_max_numberof_processor=4L*sysconf(_SC_NPROCESSORS_CONF);
+			if (long_max_numberof_processor < long_numberof_processor) {
+				long_numberof_processor=long_max_numberof_processor;
+			}
+		}
+		numberof_processor=(long_numberof_processor<2)?1:(size_t)(long_numberof_processor);
+	}
+	else {
+		long int long_numberof_processor=sysconf(_SC_NPROCESSORS_ONLN);
+		double nloadavg[3]={0.0,0.0,0.0};
+
+		if (getloadavg(nloadavg,3) == 3) {
+			double one_over_numberof_processor=1.0/(double)(long_numberof_processor);
+			nloadavg[0] *= one_over_numberof_processor;
+			nloadavg[1] *= one_over_numberof_processor;
+			nloadavg[2] *= one_over_numberof_processor;
+			if (2.5 < (nloadavg[1])) {
+				long_numberof_processor = 1L;
+				}
+			else if ((1.75 < nloadavg[1]) || (2.5 < nloadavg[0])) {
+				long_numberof_processor /= 2L;
+				}
+		}
+		numberof_processor=(long_numberof_processor<2)?1:(size_t)(long_numberof_processor);
+	}
+
+	return numberof_processor;
+}
+
+typedef struct {
+	size_t cardinality;
+	size_t tilingWidth;
+	size_t tilingHeight;
+	int tileWidth;
+	int tileHeight;
+	int areaHorizontalLowerBound;
+	int areaVerticalLowerBound;
+	int areaHorizontalUpperBound;
+	int areaVerticalUpperBound;
+	} Script_Tile_tiling_tuple;
+
+static bool Script_Tile_tiling_compute(size_t cardinality, int xmin, int ymin, int xmax, int ymax, Script_Tile_tiling_tuple *tuple)
+{
+	bool answer = false;
+	if ((tuple) /* && (xmin<xmax) && (ymin<ymax) */) {
+		if (1<cardinality) {
+			const int tile_delta = 67;
+			const int tile_4_3_delta = (4*tile_delta)/3; // 1/3 satisfies (1+epsilon)/2 = (1-epsilon)
+			int tile_cardinal = INT_MAX;
+			int tile_xcardinal = INT_MAX;
+			int tile_ycardinal = INT_MAX;
+			int tile_xdiff = INT_MAX;
+			int tile_ydiff = INT_MAX;
+			int tile_xtile = 0;
+			int tile_ytile = 0;
+			int tile_xmin = INT_MIN;
+			int tile_ymin = INT_MIN;
+			int tile_xmax = INT_MAX;
+			int tile_ymax = INT_MAX;
+			int tile_xstep = SHRT_MAX;
+			int tile_ystep = SHRT_MAX;
+			int xtile, ytile;
+
+			tile_cardinal = (int)(cardinality);
+			tile_xdiff = xmax - xmin;
+			tile_ydiff = ymax - ymin;
+
+			if (tile_4_3_delta < tile_xdiff) {
+				if (tile_4_3_delta < tile_ydiff) { // quadratic tiling
+					tile_xcardinal = (int) ( floor ( sqrt (tile_cardinal) ) );
+					tile_ycardinal = tile_cardinal / tile_xcardinal;
+					// (tile_xcardinal <= tile_ycardinal)
+					if (tile_xdiff <= tile_ydiff) { // square or vertical rectangle
+						tile_xtile = MAX(2, tile_xdiff / tile_delta);
+						tile_ytile = MAX(2, tile_ydiff / tile_delta);
+						if (tile_xcardinal <= tile_xtile)
+							tile_xtile = tile_xcardinal;
+						else {
+							tile_ycardinal = tile_cardinal / tile_xtile;
+							//tile_xcardinal = tile_xtile;
+						}
+						if (tile_ycardinal < tile_ytile)
+							tile_ytile = tile_ycardinal;
+					}
+					else { // (tile_ydiff < tile_xdiff) // horizontal rectangle
+						if (tile_xcardinal < tile_ycardinal) {
+							xtile = tile_ycardinal; tile_ycardinal = tile_xcardinal; tile_xcardinal = xtile;
+						} // (tile_ytile <= tile_xtile)
+						tile_ytile = MAX(2, tile_ydiff / tile_delta);
+						tile_xtile = MAX(2, tile_xdiff / tile_delta);
+						if (tile_ycardinal <= tile_ytile)
+							tile_ytile = tile_ycardinal;
+						else {
+							tile_xcardinal = tile_cardinal / tile_ytile;
+							//tile_ycardinal = tile_ytile;
+						}
+						if (tile_xcardinal < tile_xtile)
+							tile_xtile = tile_xcardinal;
+					}
+				}
+				else { // horizontal linear tiling
+					tile_ytile = 1;
+					tile_xtile = MAX(2, tile_xdiff / tile_delta);
+					if (tile_cardinal < tile_xtile) tile_xtile = tile_cardinal;
+				}
+			}
+			else if (tile_4_3_delta < tile_ydiff) { // vertical linear tiling
+				tile_xtile = 1;
+				tile_ytile = MAX(2, tile_ydiff / tile_delta);
+				if (tile_cardinal < tile_ytile) tile_ytile = tile_cardinal;
+			}
+			else { // scalar tiling
+				tile_ytile = tile_xtile = 1;
+			}
+
+			tile_cardinal = tile_xtile * tile_ytile;
+			tile_xstep = tile_xdiff / tile_xtile;
+			tile_ystep = tile_ydiff / tile_ytile;
+
+			tuple->cardinality = (size_t)(tile_cardinal);
+			tuple->tilingWidth = (size_t)(tile_xtile);
+			tuple->tilingHeight = (size_t)(tile_ytile);
+			tuple->tileWidth = tile_xstep;
+			tuple->tileHeight = tile_ystep;
+			tuple->areaHorizontalLowerBound = xmin;
+			tuple->areaVerticalLowerBound = ymin;
+			tuple->areaHorizontalUpperBound = xmax;
+			tuple->areaVerticalUpperBound = ymax;
+
+			answer = (1<tile_cardinal)?true:false;
+		}
+		else {
+			tuple->cardinality = 1;
+			tuple->tilingWidth = 1;
+			tuple->tilingHeight = 1;
+			tuple->tileWidth = xmax - xmin;
+			tuple->tileHeight = ymax -ymin;
+			tuple->areaHorizontalLowerBound = xmin;
+			tuple->areaVerticalLowerBound = ymin;
+			tuple->areaHorizontalUpperBound = xmax;
+			tuple->areaVerticalUpperBound = ymax;
+			//answer = false;
+		}
+	}
+	return answer;
+}
+
+typedef struct {
+	size_t tdx;
+	SurfaceCalc *sc;
+	int xmin;
+	int ymin;
+	int xmax;
+	int ymax;
+	RgbBuffer *buffer;
+	bool wipedOff;
+	} Script_Tile_tile_tuple;
+
+static void * Script_Tile_tile_tile(void * tuple)
+{
+	Script_Tile_tile_tuple * tile_tuple = (Script_Tile_tile_tuple *)(tuple);
+	SurfaceCalc * sc = tile_tuple->sc;
+	if (tile_tuple->wipedOff) sc->wipeOff(); else tile_tuple->wipedOff = false;
+	(sc->surface_calculate)( tile_tuple->xmin, tile_tuple->ymin, tile_tuple->xmax, tile_tuple->ymax, *(tile_tuple->buffer));
+	pthread_exit(NULL);
+	return NULL;
+}
+
+#endif
+
 void Script::drawSurface()
 {
+#ifdef ENABLE_PARALLELIZATION_SURFACE
+	bool flag_multithreading = false;
+#endif
 	int xmin, ymin;
 	int xmax, ymax;
 	checkVariables();
@@ -437,12 +637,99 @@
 
 	Y_AXIS_LR_ROTATE = 0.0;
 
+#ifdef ENABLE_PARALLELIZATION_SURFACE
+
+	Script_Tile_tiling_tuple tiling_tuple;
+
+	flag_multithreading = Script_Tile_tiling_compute(Script_Tile_sys_numberof_processors(), xmin, ymin, xmax, ymax, &tiling_tuple);
+
+#ifdef ENABLE_PARALLELIZATION_SURFACE_PTHREAD
+	pthread_t listof_thread[tiling_tuple.cardinality];
+#endif
+	SurfaceCalc *arrayof_sc = new SurfaceCalc [tiling_tuple.cardinality];
+	Script_Tile_tile_tuple listof_tile_tuple[tiling_tuple.cardinality];
+	Script_Tile_tile_tuple * tile_tuple = NULL;
+	size_t tdx = 0;
+
+//	#pragma omp parallel shared(arrayof_sc)
+//	#pragma omp parallel shared(listof_tile_tuple)
+
+	if (flag_multithreading) {
+
+#define PROCNAME_SURFACECALCULATE(XTile,YTile,TileXmin,TileYmin,TileXmax,TileYmax) \
+	/* fprintf(stderr, \
+			"::::: DEBUG:%s:%s:%d:\t<%2zu/%zu>\t[%2zu/%zu,%2zu/%zu] (%4i,%4i) (%4i,%4i)\n", \
+		__FILE__,__func__,__LINE__, \
+		(tdx + 1), tiling_tuple.cardinality, \
+		XTile, tile_xtile, YTile, tile_ytile, \
+		TileXmin, TileYmin, TileXmax, TileYmax \
+		); */ \
+	tile_tuple = listof_tile_tuple + tdx; \
+	tile_tuple->tdx = tdx; \
+	tile_tuple->sc = arrayof_sc + tdx; \
+	(tile_tuple->sc)->setDisplay (getDisplay()); \
+	(tile_tuple->sc)->setPreview (getPreview()); \
+	tile_tuple->xmin = TileXmin; tile_tuple->ymin = TileYmin; \
+	tile_tuple->xmax = TileXmax; tile_tuple->ymax = TileYmax; \
+	tile_tuple->buffer = intensity; \
+	tile_tuple->wipedOff = true; \
+	++tdx;
+
+		{ // WIP: might be macrofied
+			const size_t tile_xtile = tiling_tuple.tilingWidth;
+			const size_t tile_ytile = tiling_tuple.tilingHeight;
+			const int tile_xstep = tiling_tuple.tileWidth;
+			const int tile_ystep = tiling_tuple.tileHeight;
+			int tile_xmin = INT_MAX;
+			int tile_ymin = INT_MAX;
+			int tile_xmax = INT_MIN;
+			int tile_ymax = INT_MIN;
+			size_t xtile = 0;
+			size_t ytile = 0;
+			for( ytile=1, tile_ymin=ymin, tile_ymax=tile_ymin+tile_ystep; ytile<tile_ytile; tile_ymin=tile_ymax,tile_ymax+=tile_ystep,++ytile) {
+				for( xtile=1, tile_xmin=xmin, tile_xmax=tile_xmin+tile_xstep; xtile<tile_xtile; tile_xmin=tile_xmax,tile_xmax+=tile_xstep,++xtile) {
+					PROCNAME_SURFACECALCULATE(xtile, ytile, tile_xmin, tile_ymin, tile_xmax, tile_ymax);
+				}
+				PROCNAME_SURFACECALCULATE(xtile, ytile, tile_xmin, tile_ymin, xmax, tile_ymax);
+			}
+			for( xtile=1, tile_xmin=xmin, tile_xmax=tile_xmin+tile_xstep; xtile<tile_xtile; tile_xmin=tile_xmax,tile_xmax+=tile_xstep,++xtile) {
+				PROCNAME_SURFACECALCULATE(xtile, ytile, tile_xmin, tile_ymin, tile_xmax, ymax);
+			}
+			PROCNAME_SURFACECALCULATE(xtile, ytile, tile_xmin, tile_ymin, xmax, ymax);
+		}
+
+#ifdef ENABLE_PARALLELIZATION_SURFACE_PTHREAD
+		for (tdx=0,tile_tuple=listof_tile_tuple;tdx<tiling_tuple.cardinality;++tdx,++tile_tuple) {
+			pthread_create(&listof_thread[tdx],NULL,Script_Tile_tile_tile,(void *)(tile_tuple));
+		}
+		for (tdx=0;tdx<tiling_tuple.cardinality;++tdx) {
+			pthread_join(listof_thread[tdx],NULL);
+		}
+#else
+		for (tdx=0,tile_tuple=listof_tile_tuple;tdx<tiling_tuple.cardinality;++tdx,++tile_tuple) {
+			(tile_tuple->sc->surface_calculate)( tile_tuple->xmin, tile_tuple->ymin, tile_tuple->xmax, tile_tuple->ymax, *intensity);
+		}
+#endif
+
+#undef PROCNAME_SURFACECALCULATE
+	}
+	else {
+		(arrayof_sc->setDisplay)(getDisplay());
+		(arrayof_sc->setPreview)(getPreview());
+		(arrayof_sc->surface_calculate)(xmin, ymin, xmax, ymax, *intensity);
+	}
+
+#else
 	SurfaceCalc sc;
 	sc.setDisplay (getDisplay());
 	sc.setPreview (getPreview());
+#endif
 
+#ifdef ENABLE_PARALLELIZATION_SURFACE
+#else
 	//sc.wipeOff();
 	sc.surface_calculate(xmin, ymin, xmax, ymax, *intensity);
+#endif
 
 	if( display_numeric.stereo_eye ) {
 
@@ -450,8 +737,28 @@
 		*getZBuffer3d() = -10.0;
 		intensity->StereoLeft();
 
+#ifdef ENABLE_PARALLELIZATION_SURFACE
+	if (flag_multithreading) {
+#ifdef ENABLE_PARALLELIZATION_SURFACE_PTHREAD
+		for (tdx=0,tile_tuple=listof_tile_tuple;tdx<tiling_tuple.cardinality;++tdx,++tile_tuple) {
+			pthread_create(&listof_thread[tdx],NULL,Script_Tile_tile_tile,(void *)(tile_tuple));
+			}
+		for (tdx=0;tdx<tiling_tuple.cardinality;++tdx) {
+			pthread_join(listof_thread[tdx],NULL);
+		}
+#else
+		for (tdx=0,tile_tuple=listof_tile_tuple;tdx<tiling_tuple.cardinality;++tdx,++tile_tuple) {
+			(tile_tuple->sc->surface_calculate)( tile_tuple->xmin, tile_tuple->ymin, tile_tuple->xmax, tile_tuple->ymax, *(tile_tuple->buffer));
+		}
+#endif
+	}
+	else {
+		(arrayof_sc->surface_calculate)(xmin, ymin, xmax, ymax, *intensity);
+	}
+#else
 		sc.wipeOff();
 		sc.surface_calculate(xmin, ymin, xmax, ymax, *intensity);
+#endif
 
 		int     back =(int)( 0.299*((float)(color_background_data[RED]))
 				     +0.587*((float)(color_background_data[GREEN]))
@@ -469,6 +776,10 @@
 			getDisplay()->drawRgbBuffer (*intensity);
 	}
 
+#ifdef ENABLE_PARALLELIZATION_SURFACE
+	delete[] arrayof_sc;
+#endif
+
 }
 
 
--- a/src/cube.h
+++ b/src/cube.h
@@ -50,7 +50,7 @@
 	double    equations[6][4];
 } cube;
 
-extern  cube    clip_cube;
+extern  cube    GLOBAL_clip_cube;
 
 void    init_cube( cube* );
 
--- a/src/cube.cc
+++ b/src/cube.cc
@@ -50,7 +50,7 @@
 //  global data
 // -------------
 
-cube    clip_cube;
+cube    GLOBAL_clip_cube;
 
 // -----------------------------------
 //  rotate the cube around the x-axis
--- a/src/SurfaceDataStruct.h
+++ b/src/SurfaceDataStruct.h
@@ -84,7 +84,6 @@
 	int    root_n_all;  /* sk Anzahl aller vorhandenen NUS aller Flchen*/        
 	bool   init_all;    /* sk wurde berhaupt eine Flche initialisiert*/
 
-
 };
 
 #endif
--- a/src/SurfaceDataStruct.cc
+++ b/src/SurfaceDataStruct.cc
@@ -73,7 +73,7 @@
 	//  init clipping_cube
 	// --------------------
 	
-	init_cube( &clip_cube );
+//	init_cube( &GLOBAL_clip_cube );
 
 	surf_n=k;
 	init_all = k > 0;
--- a/draw/SurfaceCalc.cc
+++ b/draw/SurfaceCalc.cc
@@ -182,7 +182,8 @@
 							zero_of_polyx( &sf_ds.formulas[i].hpxyz->pZ,
 								       zmax,zmin,root[i],			
 								       root_last_pixel[i],root_last_pixel_n[i],
-                                                                       clipper );
+                       clipper,
+											 flag_numeric_multi_root );
 						sf_ds.root_n_all += sf_ds.surfaces[i].root_n;
 						root_last_pixel_n[i] = sf_ds.surfaces[i].root_n;     
 					} 
@@ -442,7 +443,7 @@
 										&sf_ds.formulas[i].hpxyz->pZ,
 										zmax,zmin,root[i],
 										root_last_pixel[i],
-										root_last_pixel_n[i],clipper );
+										root_last_pixel_n[i],clipper,flag_numeric_multi_root );
 								    sf_ds.root_n_all+=sf_ds.surfaces[i].root_n;
 								    root_last_pixel_n[i] = sf_ds.surfaces[i].root_n;
 								} 
@@ -836,8 +837,7 @@
 					   double stepwidth,
 					   colorrgb& color_intens ) 		
 {
-	// FIXME
-	static    int     inside;         // true if normal points into monitor
+	int     inside;         // true if normal points into monitor
 
 	Vector  P;              // our point on the surface
 	Vector  N;              // the surface normal in P
@@ -1148,7 +1148,7 @@
 							&sf_ds.formulas[i].hpxyz->pZ,	
 							zmax,zmin,root[i],			
 							root_last_pixel[i],
-							root_last_pixel_n[i],clipper );
+							root_last_pixel_n[i],clipper,flag_numeric_multi_root );
 						sf_ds.root_n_all += sf_ds.surfaces[i].root_n; 	
 						root_last_pixel_n[i] = sf_ds.surfaces[i].root_n;
 					}
@@ -1325,14 +1325,14 @@
 	packLights( ); 
 	
 	// ------------------------
-	//  illumiation components
+	//  illumination components
 	// ------------------------
 	flag_ambient      = light_illumination_data & light_illumination_ambient_data;
 	flag_diffuse      = light_illumination_data & light_illumination_diffuse_data;
 	flag_reflected    = light_illumination_data & light_illumination_reflected_data;
 	flag_transmitted  = light_illumination_data & light_illumination_transmitted_data;
 	
-	numeric_multi_root = flag_transmitted;
+	flag_numeric_multi_root = flag_transmitted;
 	
 	part_A= (float)(((double)(light_settings[0].ambient))/100.0);
 		
--- a/draw/SurfaceCalc.h
+++ b/draw/SurfaceCalc.h
@@ -137,11 +137,13 @@
 	int      flag_reflected;
 	int      flag_transmitted;
 
+	int      flag_numeric_multi_root;
+
 	float    ambient[3][3];   	// sk :1'ten 3 definiert RGB, ...
 	float    diffuse[3][3];   	//     2'ten 3 surfs,inside,balls
 	float    tocolor[3];            // sk :Farbe des Farbverlaufs
 
-	double   reflected;       			
+	double   reflected;
 	double   transmitted;
 	int      thickness;
 	double   transmission;
--- a/src/roots.h
+++ b/src/roots.h
@@ -41,9 +41,8 @@
 //  Function protoytpes
 // ----------------------------------------------------------------------------
 
-extern int numeric_multi_root;
 class polyx;
-int     zero_of_polyx( polyx*,double,double,double*,double*,int,NewClip* );
+int zero_of_polyx( polyx*,double,double,double*,double*,int,NewClip*,int );
 
 #endif  /* ROOTS_H */
 
--- a/src/roots.cc
+++ b/src/roots.cc
@@ -47,9 +47,6 @@
 #include "gui_config.h"
 #include "NewClip.h"
 
-//  using namespace globals;
-int numeric_multi_root = 0;
-
 // ----------------------------------------------------------------------------
 //  compute roots of f in [zmin,zmax] where g is positive
 //      f:           polynomial
@@ -63,11 +60,12 @@
 //      numeric_root_finder_data: chooses method
 // ----------------------------------------------------------------------------
 
-int     zero_of_polyx( polyx *f,
+int zero_of_polyx( polyx *f,
                 double zmax,double zmin,
                 double *root,
                 double *root_old,int old_n,
-                NewClip *clipper )
+                NewClip *clipper,
+								int numeric_multi_root)
 {