CUDA approfondi

Dans le dossier common, ouvrez les solutions cutil.sln et paramgl.sln. Compilez les configurations Debug et Release.

Ensuite, dans le dossier cudpp, ouvrez la solution cudpp.sln. Sélectionnez la configuration qui vous intéresse (un conseil : compilez au moins une version Debug et une version Release, pour faciliter le débogage de votre application), puis lancez la compilation, comme pour toute solution.

Ensuite, il faut installer ces librairies. Copiez les fichiers .lib du dossier lib dans votre dossier CUDA\lib, les .h du dossier inc dans CUDA\include et les fichiers .dll du dossier lib dans votre dossier CUDA\bin.

Si vous n'arrivez pas à compiler correctement ces librairies, téléchargez cette archive, qui les contient toutes (vous n'avez qu'à l'extraire dans le dossier Visual Studio 9.0\VC). Les fichiers compilés sont prévus pour Visual Studio 9 et CUDA 2.2, toute utilisation avec d'autres versions n'est pas garantie de fonctionnement.

L'utilisation des compilateurs Intel est déconseillée dans ce cas, car nvcc n'a pas encore été testé avec lui ! Le résultat de la compilation peut être très aléatoire !

Pendant la compilation, vous verrez beaucoup d'avertissements : ils ne sont que les signes de l'optimiseur du compilateur, qui optimise tout ce qu'il peut, surtout du côté GPU.

Dans un shell Linux (sh, par exemple), lancez ces quelques commandes.

( cd common ; make cuda-install=/usr/local/cuda )
( cd cudpp  ; make cuda-install=/usr/local/cuda )

Elles compilent et installent le tout automatiquement, dans le dossier de CUDA.

Tout d'abord, qu'est-il censé effectuer ? Un scan, aussi dénommé somme partielle, somme des préfixes ou réduction des préfixes. Mais qu'est-ce donc ?

Prenez une liste. Appliquez-lui un scan. Il en résultera une autre liste, dont chaque élément est la somme des précédents.

1
2
4
8
16
14
2004
0

1
3
7
15
31
45
2049
2049

 

Sélectionnez

2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.

CUT_DEVICE_INIT();

unsigned int numElements = 32768;
unsigned int memSize = sizeof(float) * numElements;


float* h_idata = (float *) malloc(memSize);

for (unsigned int i = 0; i < numElements; ++i) 
{
    h_idata[i] = (float) (rand() & 0xf);
}

CUTIL initialise un périphérique, on spécifie le nombre d'éléments à traiter puis la mémoire utilisée pour stockerces éléments. On initialise le tableau qui va les contenir et on le remplit de données aléatoires (entre 0 et 15).

 

Sélectionnez

17.
18.
19.
20.
21.
22.
23.
24.
25.

float* d_idata;
CUDA_SAFE_CALL( cudaMalloc( (void **) & d_idata, memSize));

CUDA_SAFE_CALL( cudaMemcpy( d_idata, h_idata, memSize,
                            cudaMemcpyHostToDevice) );


float* d_odata;
CUDA_SAFE_CALL( cudaMalloc( (void **) & d_odata, memSize));

On crée un tableau semblable sur le GPU et on y transfère les données. On prépare aussi un emplacement pour le résultat.

Maintenant, nous pouvons enfin commencer à parler de CUDPP.

Nous devons commencer par configurer CUDPP correctement pour scanner notre tableau. La configuration des algorithmes de CUDPP repose sur le concept des plans.

Un plan est, tout comme dans CUFFT, une structure qui contient les informations sur les algorithmes à employer pour le calcul.

Ici, nous allons effectuer un scan sur des flottants, exclusivement vers l'avant, d'un nombre d'éléments connu d'avance.

 

Sélectionnez

27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.

CUDPPConfiguration config;
config.op = CUDPP_ADD;
config.datatype = CUDPP_FLOAT;
config.algorithm = CUDPP_SCAN;
config.options = CUDPP_OPTION_FORWARD | CUDPP_OPTION_EXCLUSIVE;

CUDPPHandle scanplan = 0;
CUDPPResult result = cudppPlan(&scanplan, config, numElements, 1, 0);

if (CUDPP_SUCCESS != result)
{
    printf("Error creating CUDPPPlan\n");
    exit(-1);
}

 

Sélectionnez

43.

cudppScan(scanplan, d_odata, d_idata, numElements);

CUDPP effectue le scan, enfin.

 

Sélectionnez

46.
47.
48.
49.
50.
51.
52.
53.
54.
55.
56.

float* h_odata = (float*) malloc( memSize);  

CUDA_SAFE_CALL( cudaMemcpy( h_odata, d_odata, memSize,  
                            cudaMemcpyDeviceToHost) );  

float* reference = (float*) malloc( memSize);  
computeSumScanGold( reference, h_idata, numElements, config);  


CUTBoolean res = cutComparef( reference, h_odata, numElements);  
printf( "Test %s\n", (1 == res) ? "PASSED" : "FAILED");

On récupère le résultat du calcul, on effectue le même calcul en local (sur le CPU) puis on compare les résultats : s'ils sont identiques, le test est réussi.

 

Sélectionnez

58.
59.
60.
61.
62.
63.
64.
65.
66.
67.
68.
69.
70.

result = cudppDestroyPlan(scanplan);  
if (CUDPP_SUCCESS != result)  
{  
    printf("Error destroying CUDPPPlan\n");  
    exit(-1);  
}  

free(h_idata);  
free(h_odata);  
free(reference);  
CUDA_SAFE_CALL(cudaFree(d_idata));  
CUDA_SAFE_CALL(cudaFree(d_odata));  
}

On détruit le plan, puis on libère la mémoire allouée.

Maintenant, vous pouvez voir la simplicité de CUDPP : on alloue de la mémoire, on la remplit, on configure CUDPP en lui disant ce qu'il doit faire, puis il le fait et on récupère le résultat. C'est tout. Et c'est optimisé.

Premièrement, nous utilisons une variable de type CUDPPAlgorithm pour configurer un algorithme. Il s'agit, justement, de l'algorithme à employer.

• CUDPP_SCAN
• CUDPP_SEGMENTED_SCAN
• CUDPP_COMPACT
• CUDPP_REDUCE
• CUDPP_SORT_RADIX
• CUDPP_SORT_RADIX_GLOBAL
• CUDPP_SPMVMULT
• CUDPP_SORT_INVALID

Évidemment, cet algorithme devra correspondre à la fonction appelée plus tard avec cette configuration.

Deuxièmement, vous pouvez spécifier un opérateur, l'opération qui sera effectuée, uniquement lors d'un scan. Pour les autres algorithmes, laissez ceci vide.

Tous ces opérateurs sont binairement associables.

• CUDPP_ADD
• CUDPP_MULTIPLY
• CUDPP_MIN
• CUDPP_MAX

Troisièmement, vous devez définir le type des données en entrée et en sortie. Voici un tableau récapitulatif de ces valeurs, avec leur équivalent en C.

La configuration d'un algorithme se fait, finalement, aussi grâce à un entier. Cet entier est constitué d'une combinaison de paramètres, grâce à l'opérateur pipe (|).

CUDPPResult cudppScan(CUDPPHandle   planHandle,
                      void        * d_out, 
                      const void  * d_in, 
                      size_t        numElements);

Les données à traiter se situent sur le GPU, à l'emplacement pointé par d_in ; la sortie se situe aussi sur le GPU, à l'emplacement pointé par d_out. La fonction utilisera les options prévues dans le plan planHandle et scannera numElements éléments.

CUDPPResult cudppSegmentedScan(CUDPPHandle          planHandle,
                               void               * d_out, 
                               const void         * d_idata,
                               const unsigned int * d_iflags,
                               size_t               numElements);

La fonction traitera les segments précisés dans d_iflags de d_idata.

CUDPPResult cudppMultiScan(CUDPPHandle planHandle,
                           void        *d_out, 
                           const void  *d_in, 
                           size_t      numElements,
                           size_t      numRows);

La fonction effectuera numRows scans sur les données avant d'écrire le résultat.

CUDPPResult cudppCompact(CUDPPHandle          planHandle,
                         void               * d_out, 
                         size_t             * d_numValidElements,
                         const void         * d_in, 
                         const unsigned int * d_isValid,
                         size_t               numElements);

Cette fonction compacte un tableau, avec les éléments d_isValid définis comme valides, les autres étant à éjecter. d_numValidElements sera le nombre d'éléments effectivement gardés.

Par exemple :

d_in        = [ a b c d e f ]
deviceValid = [ 1 0 1 1 0 1 ]
d_out       = [ a c d f ]

CUDPPResult cudppSort(CUDPPHandle planHandle,
                      void        *d_out, 
                      const void  *d_in,
                      size_t      numElements);

Trie le tableau.

Ces algorithmes appliquent une même opération à tous les éléments d'une entrée et les stockent à un autre endroit.

 

Cacher/Afficher le codeSélectionnez

#include <thrust/device_vector.h>
#include <thrust/transform.h>
#include <thrust/sequence.h>
#include <thrust/copy.h>
#include <thrust/fill.h>
#include <thrust/replace.h>
#include <thrust/functional.h>
#include <iostream>

int main(void)
{
    thrust::device_vector<int> X(10);
    thrust::device_vector<int> Y(10);
    thrust::device_vector<int> Z(10);

    // remplit X par des valeurs de 1 à 10
    thrust::sequence(X.begin(), X.end());

    // Y = -X
    thrust::transform(X.begin(), X.end(), Y.begin(), thrust::negate<int>());

    // remplit Z avec des 2
    thrust::fill(Z.begin(), Z.end(), 2);

    // Y = X mod 2
    thrust::transform(X.begin(), X.end(), Z.begin(), Y.begin(), thrust::modulus<int>());

    // remplace tous les 1 de Y par  10
    thrust::replace(Y.begin(), Y.end(), 1, 10);

    // copie le contenu de Y dans le flux de sortie cout
    thrust::copy(Y.begin(), Y.end(), std::ostream_iterator<int>(std::cout, "\n"));

    return 0;
}

thrust::modulus et thrust::negate sont implémentés comme des foncteurs, dans <thrust/functional.h>.

Nous allons poursuivre avec une implémentation personnelle d'une fonction de BLAS : SAXPY, qui effectue les opérations de la forme .

 

Cacher/Afficher le codeSélectionnez

include <thrust/transform.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/functional.h> 
#include <iostream>
#include <iterator>
#include <algorithm>

struct saxpy_functor
{
    const float a;
    saxpy_functor(float _a) : a(_a) {}
    __host__ __device__ float operator() (const float& x, const float& y) const
    {
        return a * x + y;
    }
};

void saxpy_fast(float A, thrust::device_vector<float>& X, thrust::device_vector<float>& Y)
{
   thrust::transform(X.begin(), X.end(), Y.begin(), Y.begin(), saxpy_functor(A));
}

void saxpy_slow(float A, thrust::device_vector<float>& X, thrust::device_vector<float>& Y)
{
    thrust::device_vector<float> temp(X.size());
    thrust::fill(temp.begin(), temp.end(), A);
    thrust::transform(X.begin(), X.end(), temp.begin(), temp.begin(), thrust::multiplies<float>());
    thrust::transform(temp.begin(), temp.end(), Y.begin(), Y.begin(), thrust::plus<float>());
}

Vous pouvez voir deux implémentations : une rapide et une lente. La première effectue lectures et écritures. La seconde effectue lectures et écritures.

Vu que SAXPY est limité par la mémoire et non par la puissance de calcul, la première est plus rapide que la seconde.

La technique utilisée s'appelle la réduction des kernels, car, pour obtenir de meilleures performances, dans ce cas, vous devez effectuer le moins possible de lectures et d'écritures en mémoire.

Cette manière de faire (avec thrust::transform) limite les arguments : un ou deux en entrée. Si vous avez besoin de plus d'arguments, voyez thrust::counting_iterator et thrust::for_each tels qu'utilisés dans ce nouvel exemple.

 

Cacher/Afficher le codeSélectionnez

#include <thrust/for_each.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/generate.h>
#include <thrust/iterator/iterator_traits.h>
#include <thrust/device_ptr.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <iostream>
#include <cstdlib>

struct complicated_functor
{
    const float * A;
    const float * B;
    const float * C;
          float * D;

    complicated_functor(const float * _A, const float * _B, const float * _C, float * _D):
        A(_A), B(_B), C(_C), D(_D) {}

    template <typename IndexType> __host__ __device__
    void operator()(const IndexType & i)
    {
        D[i] = A[i] + B[i] * C[i];
    }
};


int main(void)
{
    thrust::device_vector<float> A(5);
    thrust::device_vector<float> B(5);
    thrust::device_vector<float> C(5);
    thrust::device_vector<float> D(5);

    A[0] = 3;  B[0] = 6;  C[0] = 2; 
    A[1] = 4;  B[1] = 7;  C[1] = 5; 
    A[2] = 0;  B[2] = 2;  C[2] = 7; 
    A[3] = 8;  B[3] = 1;  C[3] = 4; 
    A[4] = 2;  B[4] = 8;  C[4] = 3; 

    thrust::experimental::counting_iterator<int, thrust::random_access_device_iterator_tag> begin(0);

    complicated_functor op(thrust::raw_pointer_cast(&A[0]), 
                           thrust::raw_pointer_cast(&B[0]), 
                           thrust::raw_pointer_cast(&C[0]),
                           thrust::raw_pointer_cast(&D[0]));

    thrust::for_each(begin, begin + 5, op);

    for(int i = 0; i < 5; i++)
        std::cout << A[i] << " + " << B[i] << " * " << C[i] << " = " << D[i] << std::endl;
}

Un algorithme de réduction utilise une opération binaire pour réduire une entrée à une valeur. Voici la manière d'obtenir la somme d'un tableau avec Thrust.

int sum = thrust::reduce(D.begin(), D.end(), (int) 0, thrust::plus<int>());

Cette manière de procéder convient généralement, mais pas toujours, par exemple pour compter. Thrust, comme la STL, fournit quelques autres fonctions, pour votre facilité.

#include <thrust/count.h>
#include <thrust/device_vector.h>

thrust::device_vector<int> vec(5,0);
vec[1] = 1;
vec[3] = 1;
vec[4] = 1;

int result = thrust::count(vec.begin(), vec.end(), 1);
// result == 3

On peut utiliser ces réductions pour calculer la norme d'un vecteur.

 

Cacher/Afficher le codeSélectionnez

#include <thrust/transform_reduce.h>
#include <thrust/functional.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <cmath>

template <typename T>
struct square
{
    __host__ __device__
        T operator()(const T& x) const { 
            return x * x;
        }
};

int main(void)
{
    float x[4] = {1.0, 2.0, 3.0, 4.0};

    thrust::device_vector<float> d_x(x, x + 4);

    square<float>        unary_op;
    thrust::plus<float> binary_op;
    float init = 0;

    float norm = std::sqrt( thrust::transform_reduce(d_x.begin(), d_x.end(), unary_op, init, binary_op) );

    std::cout << norm << std::endl;

    return 0;
}

En fusionnant la multiplication avec l'addition, nous avons un kernel très optimisé, autant qu'un kernel écrit à la main et patiemment optimisé.

Thrust, une fois de plus, se calque sur la STL pour fournir ses algorithmes de tri.

 

Cacher/Afficher le codeSélectionnez

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/generate.h>
#include <thrust/sort.h>
#include <cstdlib>

int main(void)
{
    thrust::host_vector<int> h_vec(20);
    thrust::generate(h_vec.begin(), h_vec.end(), rand);

    thrust::device_vector<int> d_vec = h_vec;

    thrust::sort(d_vec.begin(), d_vec.end());

    return 0;
}

De plus, Thrust propose le tri par clé.

const int N = 6;
int    keys[N] = {  1,   4,   2,   8,   5,   7};
char values[N] = {'a', 'b', 'c', 'd', 'e', 'f'};
thrust::sort_by_key(keys, keys + N, values);
// keys   == {  1,   2,   4,   5,   7,   8}
// values == {'a', 'c', 'b', 'e', 'f', 'd'}

Et selon un critère défini par l'utilisateur.

const int N = 6;
int A[N] = {1, 4, 2, 8, 5, 7};
thrust::stable_sort(A, A + N, thrust::greater<int>());
// A == {8, 7, 5, 4, 2, 1}