 	// memset:
 	CUDA_CHECK_RETURN( cudaMemset(dev_oPLOT, 0,  (size_t)oPLOT_bytes) 		);
 
+	/*
+	 * 		****3-D aggregation****
+	 */
 	// kernel size:
 	unsigned int 	BLOCKSIZE, GRIDSIZE;
 	BLOCKSIZE	= 32;//floor(sqrt( devProp.maxThreadsPerBlock ));

		CUDA_CHECK_RETURN( cudaMemcpy(oGRID, dev_oGRID, iMap_bytes, cudaMemcpyDeviceToHost) );
		end_t = clock();
		break;
	}
	// elapsed time [ms]:
	elapsed_time += (int)( (double)(end_t  - start_t ) / (double)CLOCKS_PER_SEC * 1000 );

	// call 2d reduction to computing mean & std for scatterplot:
	dim3 block_2d( BLOCKSIZE*BLOCKSIZE,1,1);
	dim3 grid_2d ( Nmaps,1,1);

		}
	}
}
__global__ void reduction_3d_mean( const double *lin_maps, unsigned int map_len, unsigned int Nmaps, double *sum_map ){
	/*
	 * 		lin_maps:	|------|------|------|	...	|------|
	 * 					   1st	  2nd	 3rd	...	   Nth
	 *
	 * 		block:		blockDim.Y(=1) * blockDim.X(=32^2, other);
	 * 		grid:		gridDim.Y(=1)  * gridDim.X(=map_len / blockDim.X);
	 */
	unsigned int r 			= threadIdx.y;
	unsigned int c 			= threadIdx.x;
	unsigned int bdx		= blockDim.x;
	unsigned int bdy		= blockDim.y;
	unsigned int bix		= blockIdx.x;
	unsigned int biy		= blockIdx.y;
	unsigned int gdx		= gridDim.x;
	//unsigned int gdy		= gridDim.y;
	unsigned long int tid 	= (  c + bdx*( r + bdy * (bix + gdx*biy) )  );
/* 			  				     |_________|
 * 								 block-offset
 * 										     + |_____grid-offset_____|
 */
	unsigned int ii 		= 0;
	for( ii=0; ii<Nmaps-(Nmaps % 2); ii+=2 ){
		if(tid < map_len){
			sum_map[tid] = sum_map[tid] + lin_maps[tid + ii*map_len] + lin_maps[(ii+1)*map_len + tid];
		}
	}
	if(Nmaps % 2){
		if(tid < map_len){
			sum_map[tid] = sum_map[tid] + lin_maps[(Nmaps-1)*map_len + tid];
		}
	}
	if(tid < map_len){
		sum_map[tid] = sum_map[tid] / Nmaps;
	}
}

__global__ void reduction_2d_sum( const double *lin_maps, unsigned int map_len, unsigned int Nmaps, double *sum_map_2d ){
	/*
	 * 		lin_maps:	|------|------|------|	...	|------|
