La fonction générale du programme

5.5.1.1.3 : La fonction générale du programme

Maintenant, implémentons la fonctions générale :

///Simulate the images
/**	@param nbRow : number of rows of the images to be created
 * 	@param nbCol : number of columns of the images to be created
 * 	@param nbImage : number of images to be created
 * 	@param nbExtraStep : number of extra steps to be computed between images
 * 	@param killRate : rate of the process which converts V into P
 * 	@param feedRate : rate of the process which feeds U and drains U, V and P
 * 	@param dt : time interval between two computation
 * 	@param outputFile : name of the file to be created
 * 	@param blockSizeRow : number of rows of the blocks
 * 	@param blockSizeCol : number of columns of the blocks
 * 	@param nbThread : number of thread to be used
 * 	@return true on succsess, false otherwise
*/
bool simulateImage(size_t nbRow, size_t nbCol, size_t nbImage, size_t nbExtraStep, float killRate, float feedRate, float dt, const std::string & outputFile,
		size_t blockSizeRow, size_t blockSizeCol, size_t nbThread)
{

Nous allons commencer par indiquer quelques paramètres utiles pour l'utilisateur :

1	std::cout << "simulateImage : nbImage = "<<nbImage<<", nbRow = " << nbRow << ", nbCol = " << nbCol << std::endl;

Maintenant, nous allons spécifier à TBB combien de thread on veut utiliser :

	int defaultNbThread = tbb::task_scheduler_init::default_num_threads();
	if(nbThread != 0lu){
		defaultNbThread = (int)nbThread;
	}
	tbb::task_scheduler_init init(defaultNbThread);

Cette implémentation est moche et va produire un warning car le tbb::task_scheduler_init est déprécié. Et il faudrait utilise explicitement TBB pour gérer cela, ce qui ne serait plus du C++20. Mais pour le moment, le C++20 ne permet pas de faire cela (en tout cas l'implémmentation courrante de GCC 9.3). Mais dans le futur on pourra spécifier le nombre de thread avec un executor et ce sera bien plus classe.

Ensuite, nous devons définir l'ensemble de nos images en HDF5 :

1
2
3

	MatrixHdf5 fullMat;
	fullMat.setAllDim(nbCol, nbRow);
	fullMat.resize(nbImage);

On créé et on initialise les PTensor et les pointeurs de nos images temporaires :

1
2
3

	PTensor<float> tmpInU, tmpInV, tmpOutU, tmpOutV;
	float *tmpU1 = NULL, *tmpU2 = NULL, *tmpV1 = NULL, *tmpV2 = NULL;
	allocate_temporary(tmpU1, tmpU2, tmpV1, tmpV2, tmpInU, tmpInV, tmpOutU, tmpOutV, nbRow, nbCol);

On définit notre stencil 3x3 (avec éventuellement des coéfficients plus sérieux mais qui donnent des résultats moins rigolos) :

	long nbStencilRow(3l), nbStencilCol(3l);
	
	float diffudionRateU(0.1f), diffusionRateV(0.05f);
	//This matrix of neigbour exchange is quite accurate but gives not so fun results
// 	float matDeltaSquare[] = 	{0.05f, 0.2f, 0.05f,
// 					0.2f, 0.0f, 0.2f,
// 					0.05f, 0.2f, 0.05f};
	float matDeltaSquare[] = 	{1.0f, 1.0f, 1.0f,
					1.0f, 1.0f, 1.0f,
					1.0f, 1.0f, 1.0f};

Cependant, nous devons maintenant convertir tout nos temporaires afin qu'ils soient complètement vectorisables :

	//Let's convert these temporaries into intrinsics temporaries
	PTensor<float> tmpVecInU(AllocMode::ALIGNED), tmpVecInV(AllocMode::ALIGNED), tmpVecOutU(AllocMode::ALIGNED), tmpVecOutV(AllocMode::ALIGNED);
	tmpVecInU.fromScalToVecNeigbhour(tmpInU, PLIB_VECTOR_SIZE_FLOAT);
	tmpVecOutU.fromScalToVecNeigbhour(tmpOutU, PLIB_VECTOR_SIZE_FLOAT);
	tmpVecInV.fromScalToVecNeigbhour(tmpInV, PLIB_VECTOR_SIZE_FLOAT);
	tmpVecOutV.fromScalToVecNeigbhour(tmpOutV, PLIB_VECTOR_SIZE_FLOAT);

La vectorisation de la matrice qui définit les poids du gradient :

1
2

	PTensor<float> vecMatDeltaSquare(AllocMode::ALIGNED, nbStencilRow, nbStencilCol*PLIB_VECTOR_SIZE_FLOAT);
	reshuffle_broadcastTensor(vecMatDeltaSquare.getData(), matDeltaSquare, nbStencilRow, nbStencilCol, 0lu, PLIB_VECTOR_SIZE_FLOAT);

On définit les pointeurs correspondants :

1	float * ptrVecMatStencil = vecMatDeltaSquare.getData();

Et le résultat de chaque image :

1	PTensor<float> tmpScalOutV(AllocMode::ALIGNED);

Nous devons également créer nos vecteur de blocs qui seront liés à leur tenseur respectif, ce qui nous économisera une copie de donnée :

	std::vector<PBlock<float> > vecBlockOutU, vecBlockOutV, vecBlockInU, vecBlockInV;
	tmpVecOutU.splitBlockLink(vecBlockOutU, blockSizeRow, blockSizeCol, 1lu);
	tmpVecOutV.splitBlockLink(vecBlockOutV, blockSizeRow, blockSizeCol, 1lu);
	tmpVecInU.splitBlockLink(vecBlockInU, blockSizeRow, blockSizeCol, 1lu);
	tmpVecInV.splitBlockLink(vecBlockInV, blockSizeRow, blockSizeCol, 1lu);

On défint la nouvelle taille des matrices :

1	size_t nbVecRow(tmpVecInV.getFullNbRow()), nbVecCol(tmpVecInV.getNbCol());

On créé notre barre de chargement :

1 2	ProgressTime progress(nbImage); progress.start();

On boucle sur toutes les images que l'on veut créer :

1	for(size_t i(0lu); i < nbImage; ++i){

Petite mise à jour de la barre de chargement :

1	progress.print();

On boucle sur les image supplémentaires à calculer entre deux images à sauvegarder :

1	for(size_t j(0lu); j < nbExtraStep/2lu; ++j){

On appelle notre calcul de Gray Scott :

			grayscott_propagation_link_block_parallel(
					vecBlockOutU, vecBlockOutV, vecBlockInU, vecBlockInV,
					ptrVecMatStencil, nbStencilRow, nbStencilCol,
					diffudionRateU, diffusionRateV, feedRate, killRate, dt);

Il ne faut pas oublier que notre calcul se fait que sur la moitier des valeurs duppliquées, il faut donc les mettre à jour :

1
2
3

			//Let's update the dupplicated values
			reshuffle_updateDupplicateVecNeighbour(tmpVecOutU.getData(), nbVecRow, nbVecCol, PLIB_VECTOR_SIZE_FLOAT);
			reshuffle_updateDupplicateVecNeighbour(tmpVecOutV.getData(), nbVecRow, nbVecCol, PLIB_VECTOR_SIZE_FLOAT);

Plutôt que d'intervertir les pointeurs input et output, il est plus simple de faire deux fois moins d'étapes doubles :

			///Let's swap the in and out tensors
			grayscott_propagation_link_block_parallel(
					vecBlockInU, vecBlockInV, vecBlockOutU, vecBlockOutV,
					ptrVecMatStencil, nbStencilRow, nbStencilCol,
					diffudionRateU, diffusionRateV, feedRate, killRate, dt);

Il ne faut pas oublier que notre calcul se fait que sur la moitier des valeurs duppliquées, il faut donc les mettre à jour :

			//Let's update the dupplicated values
			reshuffle_updateDupplicateVecNeighbour(tmpVecInU.getData(), nbVecRow, nbVecCol, PLIB_VECTOR_SIZE_FLOAT);
			reshuffle_updateDupplicateVecNeighbour(tmpVecInV.getData(), nbVecRow, nbVecCol, PLIB_VECTOR_SIZE_FLOAT);
		}

On récupère l'image que l'on veut sauvegarder (en faisant attention car on interverti les pointeurs régulièrement) :

1	tmpScalOutV.fromVecToScalNeigbhour(tmpVecInV); //The pointers were swaped

On sauvegarde l'image que l'on veut :

1 2	fullMat.setRow(i, tmpScalOutV.getData()); }

On termine la barre de chargement :

1 2	progress.finish(); std::cerr << "Done" << std::endl;

` On écrit le fichier HDF5 :

1 2	//Let's save the output file fullMat.write(outputFile);

On renvoie true car tout s'est bien passé et on finit la fonction :

1 2	return true; }