<commit_before>/** @file simplex_local.h
 *  @brief simplex feasible set for petscvector with local problems
 *
 *  @author <NAME>
 */

#ifndef PASC_SIMPLEXFEASIBLESET_LOCAL_H
#define	PASC_SIMPLEXFEASIBLESET_LOCAL_H

/* for debugging, if >= 100, then print info about ach called function */
extern int DEBUG_MODE;

#include "pascinference.h"

#ifndef USE_PETSCVECTOR
 #error 'SIMPLEXFEASIBLESET_LOCAL is for PETSCVECTOR'
#endif

typedef petscvector::PetscVector PetscVector;

#ifdef USE_CUDA
    #include <../src/vec/vec/impls/seq/seqcuda/cudavecimpl.h>
#endif

namespace pascinference {


/** \class SimplexFeasibleSet_Local
 *  \brief class for manipulation with simplex set
 *
 *  Provides structure for manipulation with simplex feasible set, i.e.
 * \f[
 * \Omega = 
 *  \left\lbrace x \in R^{KT}: 
 *    \forall t = 0,\dots,T-1: \sum\limits_{k=0}^{K-1} x_{t+kT} = 1,  
 *    x \geq 0
 *  \right\rbrace
 *	\f]
 * 
*/
class SimplexFeasibleSet_Local: public GeneralFeasibleSet<PetscVector> {
	private:
		
		/** @brief sort array using bubble sort
		 * 
		 * @param x array with values
		 * @param n size of array
		*/ 		
		void sort_bubble(double *x, int n);

		/** @brief compute projection to one simplex subset
		 * 
		 * @param x_sub values of subvector in array
		 * @param n size of subvector
		*/ 		
		void project_sub(int t, double *x, int T, int K);

		#ifdef USE_CUDA
			double *x_sorted;
		#endif
		
	public:
		/** @brief default constructor
		*/ 	
		SimplexFeasibleSet_Local(int T, int K_local);
		
		/** @brief default destructor
		 */ 
		~SimplexFeasibleSet_Local();

		/** @brief print properties of this feasible set
		 * 
		 * @param output where to print
		 */ 
		void print(ConsoleOutput &output) const;

		/** @brief get name of this feasible set
		 */
		virtual std::string get_name() const;

		/* variables */
		int T; /**< number of disjoint simplex subsets */
		int K_local; /**< size of each simplex subset */
		
		/** @brief compute projection onto feasible set
		 * 
		 * @param x point which will be projected
		 */		
		void project(GeneralVector<PetscVector> &x);


		
		
};

/* cuda kernels cannot be a member of class */
#ifdef USE_CUDA
__device__ void device_sort_bubble(double *x, double *x_sorted, int t, int T, int K);
__global__ void project_kernel(double *x, double *x_sorted, int T, int K);
#endif

} // end of namespace

/* ------------- implementation ----------- */
//TODO: move to impls

namespace pascinference {

/* constructor */
SimplexFeasibleSet_Local::SimplexFeasibleSet_Local(int Tnew, int Knew){
	LOG_FUNC_BEGIN

	/* set initial content */
	this->T = Tnew;
	this->K_local = Knew;

	#ifdef USE_CUDA
		/* allocate space for sorting */
		gpuErrchk( cudaMalloc((void **)&x_sorted,K_local*T*sizeof(double)) );
	#endif

	LOG_FUNC_END
}

/* general destructor */
SimplexFeasibleSet_Local::~SimplexFeasibleSet_Local(){
	LOG_FUNC_BEGIN
	
	#ifdef USE_CUDA
		/* destroy space for sorting */
		gpuErrchk( cudaFree(&x_sorted) );
	#endif	
	
	LOG_FUNC_END	
}

/* print info about feasible set */
void SimplexFeasibleSet_Local::print(ConsoleOutput &output) const {
	LOG_FUNC_BEGIN
	
	output <<  this->get_name() << std::endl;
	
	/* give information about presence of the data */
	output <<  " - T:       " << T << std::endl;
	output <<  " - K_local: " << K_local << std::endl;

	LOG_FUNC_END
}

std::string SimplexFeasibleSet_Local::get_name() const {
	return "SimplexFeasibleSet_Local";
}

void SimplexFeasibleSet_Local::project(GeneralVector<PetscVector> &x) {
	LOG_FUNC_BEGIN
	
	/* get local array */
	double *x_arr;
	
#ifdef USE_CUDA
	Timer mytimer;
	mytimer.restart();
	mytimer.start();

	TRY( VecCUDAGetArrayReadWrite(x.get_vector(),&x_arr) );

	/* use kernel to compute projection */
	//TODO: here should be actually the comparison of Vec type! not simple use_gpu
	int blockSize; /* The launch configurator returned block size */
	int minGridSize; /* The minimum grid size needed to achieve the maximum occupancy for a full device launch */
	int gridSize; /* The actual grid size needed, based on input size */
	
	gpuErrchk( cudaOccupancyMaxPotentialBlockSize( &minGridSize, &blockSize,project_kernel, 0, 0) );
	gridSize = (T + blockSize - 1)/ blockSize;
	
	coutMaster << "minGridSize:" << minGridSize << ", gridSize: " << gridSize << ", blockSize: " << blockSize << std::endl;
	
	project_kernel<<<gridSize, blockSize>>>(x_arr,x_sorted,T,K_local);
	gpuErrchk( cudaDeviceSynchronize() );

	TRY( VecCUDARestoreArrayReadWrite(x.get_vector(),&x_arr) );

	mytimer.stop();
	
#else
	TRY( VecGetArray(x.get_vector(),&x_arr) );
	
	/* use openmp */<commit_msg>writing cuda kernel<commit_after>/** @file simplex_local.h
 *  @brief simplex feasible set for petscvector with local problems
 *
 *  @author <NAME>
 */

#ifndef PASC_SIMPLEXFEASIBLESET_LOCAL_H
#define	PASC_SIMPLEXFEASIBLESET_LOCAL_H

/* for debugging, if >= 100, then print info about ach called function */
extern int DEBUG_MODE;

#include "pascinference.h"

#ifndef USE_PETSCVECTOR
 #error 'SIMPLEXFEASIBLESET_LOCAL is for PETSCVECTOR'
#endif

typedef petscvector::PetscVector PetscVector;

#ifdef USE_CUDA
    #include <../src/vec/vec/impls/seq/seqcuda/cudavecimpl.h>
#endif

namespace pascinference {


/** \class SimplexFeasibleSet_Local
 *  \brief class for manipulation with simplex set
 *
 *  Provides structure for manipulation with simplex feasible set, i.e.
 * \f[
 * \Omega = 
 *  \left\lbrace x \in R^{KT}: 
 *    \forall t = 0,\dots,T-1: \sum\limits_{k=0}^{K-1} x_{t+kT} = 1,  
 *    x \geq 0
 *  \right\rbrace
 *	\f]
 * 
*/
class SimplexFeasibleSet_Local: public GeneralFeasibleSet<PetscVector> {
	private:
		
		/** @brief sort array using bubble sort
		 * 
		 * @param x array with values
		 * @param n size of array
		*/ 		
		void sort_bubble(double *x, int n);

		/** @brief compute projection to one simplex subset
		 * 
		 * @param x_sub values of subvector in array
		 * @param n size of subvector
		*/ 		
		void project_sub(int t, double *x, int T, int K);

		#ifdef USE_CUDA
			double *x_sorted;
		#endif
		
	public:
		/** @brief default constructor
		*/ 	
		SimplexFeasibleSet_Local(int T, int K_local);
		
		/** @brief default destructor
		 */ 
		~SimplexFeasibleSet_Local();

		/** @brief print properties of this feasible set
		 * 
		 * @param output where to print
		 */ 
		void print(ConsoleOutput &output) const;

		/** @brief get name of this feasible set
		 */
		virtual std::string get_name() const;

		/* variables */
		int T; /**< number of disjoint simplex subsets */
		int K_local; /**< size of each simplex subset */
		
		/** @brief compute projection onto feasible set
		 * 
		 * @param x point which will be projected
		 */		
		void project(GeneralVector<PetscVector> &x);


		
		
};

/* cuda kernels cannot be a member of class */
#ifdef USE_CUDA
__device__ void device_sort_bubble(double *x, double *x_sorted, int t, int T, int K);
__global__ void project_kernel(double *x, double *x_sorted, int T, int K);
#endif

} // end of namespace

/* ------------- implementation ----------- */
//TODO: move to impls

namespace pascinference {

/* constructor */
SimplexFeasibleSet_Local::SimplexFeasibleSet_Local(int Tnew, int Knew){
	LOG_FUNC_BEGIN

	/* set initial content */
	this->T = Tnew;
	this->K_local = Knew;

	#ifdef USE_CUDA
		/* allocate space for sorting */
		gpuErrchk( cudaMalloc((void **)&x_sorted,K_local*T*sizeof(double)) );
	#endif

	LOG_FUNC_END
}

/* general destructor */
SimplexFeasibleSet_Local::~SimplexFeasibleSet_Local(){
	LOG_FUNC_BEGIN
	
	#ifdef USE_CUDA
		/* destroy space for sorting */
		gpuErrchk( cudaFree(&x_sorted) );
	#endif	
	
	LOG_FUNC_END	
}

/* print info about feasible set */
void SimplexFeasibleSet_Local::print(ConsoleOutput &output) const {
	LOG_FUNC_BEGIN
	
	output <<  this->get_name() << std::endl;
	
	/* give information about presence of the data */
	output <<  " - T:       " << T << std::endl;
	output <<  " - K_local: " << K_local << std::endl;

	LOG_FUNC_END
}

std::string SimplexFeasibleSet_Local::get_name() const {
	return "SimplexFeasibleSet_Local";
}

void SimplexFeasibleSet_Local::project(GeneralVector<PetscVector> &x) {
	LOG_FUNC_BEGIN
	
	/* get local array */
	double *x_arr;
	
#ifdef USE_CUDA
	Timer mytimer;
	mytimer.restart();
	mytimer.start();

	TRY( VecCUDAGetArrayReadWrite(x.get_vector(),&x_arr) );

	/* use kernel to compute projection */
	//TODO: here should be actually the comparison of Vec type! not simple use_gpu
	int blockSize; /* The launch configurator returned block size */
	int minGridSize; /* The minimum grid size needed to achieve the maximum occupancy for a full device launch */
	int gridSize; /* The actual grid size needed, based on input size */
	
	gpuErrchk( cudaOccupancyMaxPotentialBlockSize( &minGridSize, &blockSize,project_kernel, 0, 0) );
	gridSize = (T + blockSize - 1)/ blockSize;
	
//	coutMaster << "minGridSize:" << minGridSize << ", gridSize: " << gridSize << ", blockSize: " << blockSize << std::endl;
	
	project_kernel<<<gridSize, blockSize>>>(x_arr,x_sorted,T,K_local);
	gpuErrchk( cudaDeviceSynchronize() );

	TRY( VecCUDARestoreArrayReadWrite(x.get_vector(),&x_arr) );

	mytimer.stop();
	
#else
	TRY( VecGetArray(x.get_vector(),&x_arr) );
	
	/* use openmp */
{'fraction_non_alphanumeric': 0.12546205696890628, 'fraction_numerical': 0.002609262883235486, 'mean_word_length': 3.719856336582863, 'pattern_counts': {'":': 0, '<': 83, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4022741', 'n_tokens_mistral': 3419, 'n_tokens_neox': 3117, 'n_words': 870}