Survol du Advanced Massive Parallelism (AMP)

#include <iostream>
#include <vector>
#include <chrono>
#include <string>
#include <amp.h>
using namespace concurrency;
using namespace std;
using namespace std::chrono;
void test_amp();

template <class F>
   void tester(F fct, const string &msg, ostream &out) {
      out << msg << "\n\t" << flush;
      auto avant = system_clock::now();
      fct();
      auto apres = system_clock::now();
      out << "Temps ecoule: " << duration_cast<milliseconds>(apres-avant).count() << " ms." << endl;
   }

void calculs_v0(vector<int>& vA, vector<int>& vB, vector<int>& vC, int M, int N) {
   for (int i = 0; i < M; i++)
      for (int j = 0; j < N; j++)
         vC[i * N + j] = vA[i * N + j] + vB[i * N + j];
}

void calculs_v1(vector<int>& vA, vector<int>& vB, vector<int>& vC, int M, int N) {
   array_view<int> a(M*N, vA), b(M*N, vB);
   array_view<int> c(M*N, vC);
   for (int i = 0; i < M; i++)
      for (int j = 0; j < N; j++)
         c(i * N + j) = a(i * N + j) + b(i * N + j);
}

void calculs_v2(vector<int>& vA, vector<int>& vB, vector<int>& vC, int M, int N) {
   concurrency::extent<1> e(M*N);
   array_view<int, 1> a(e, vA), b(e, vB);
   array_view<int, 1> c(e, vC);
   for (int i = 0; i < M; i++)
      for (int j = 0; j < N; j++)
         c(i * N + j) = a(i * N + j) + b(i * N + j);
}

void calculs_v3(vector<int>& vA, vector<int>& vB, vector<int>& vC, int M, int N) {
   concurrency::extent<2> e(M, N);
   array_view<int, 2> a(e, vA), b(e, vB);
   array_view<int, 2> c(e, vC);

   for (int i = 0; i < e[0]; i++)
      for (int j = 0; j < e[1]; j++)
         c(i, j) = a(i, j) + b(i, j);
}

void calculs_v4(vector<int>& vA, vector<int>& vB, vector<int>& vC, int M, int N) {
   concurrency::extent<2> e(M, N);
   array_view<int, 2> a(e, vA), b(e, vB);
   array_view<int, 2> c(e, vC);
   index<2> idx(0, 0);
   for (idx[0] = 0; idx[0] < e[0]; idx[0]++)
      for (idx[1] = 0; idx[1] < e[1]; idx[1]++)
         c[idx] = a[idx] + b[idx];
         //c(idx[0], idx[1]) = a(idx[0], idx[1]) + b(idx[0], idx[1]);
}

void calculs_v5(vector<int>& vA, vector<int>& vB, vector<int>& vC, int M, int N) {
   concurrency::extent<2> e(M, N);
   array_view<int, 2> a(e, vA), b(e, vB);
   array_view<int, 2> c(e, vC);
   c.discard_data(); // ICI: éviter de copier inutilement des données
   //
   // Le code dans la λ sera exécuté sur un accélérateur par défaut (il est possible
   // de choisir un accélérateur par programmation si on le souhaite, par exemple dans un
   // cas où on saurait que celui par défaut n'est pas idéal pour nos types de données)
   //
   // Par conséquent, cette λ est sollicitée par plusieurs threads concurrents. Le
   // paramètre «index» dans le «extent» guide le découpage
   //
   // Quelques assertions qui tiendront ici, pour tout idx passé à la lambda:
   //
   // - e.rank == idx.rank
   // - e.containx(idx) == true
   // - le nombre d'appels à la lambda sera e.size()
   //
   // Dans ce cas bien précis, e.rank vaut 2 et e.size() vaut 3 * 2, donc 6. Les six
   // threads qui appelleront la lambda utiliseront les idx suivants:
   //
   //    { 0,0 } { 0,1 } { 1,0 } { 1,1 } { 2,0 } { 2,1 }
   //
   // Notez que la capture d'une lambda assujettie à restrict(amp) ne peut se faire que
   // par valeur (une capture par référence serait illégale)
   //
   // Plus de détails: http://blogs.msdn.com/b/nativeconcurrency/archive/2011/12/19/restrict-amp-restrictions-part-0-of-n-introduction.aspx
   //
   parallel_for_each(e, [=](index<2> idx) restrict (amp) {
      c[idx] = a[idx] + b[idx];
   });
   // Le synchronize() ci-dessous complète le calcul. Cette étape est optionnelle (l'extent la
   // fait dans son destructeur), sauf si on veut traiter d'éventuelles exceptions
   c.synchronize();
}

void MatMulSeq(vector<int>& vC, const vector<int>& vA, const vector<int>& vB, int M, int N, int W) {
   for (int row = 0; row < M; row++)
      for (int col = 0; col < N; col++) {
         int sum = 0;
         for(int i = 0; i < W; i++)
            sum += vA[row * W + i] * vB[i * N + col];
         vC[row * N + col] = sum;
      }
}

void MatMulAmp(vector<int>& vC, const vector<int>& vA, const vector<int>& vB, int M, int N, int W) {
   array_view<const int, 2> a(M, W, vA), b(W, N, vB);
   array_view<int, 2> c(M, N, vC);
   c.discard_data();
   parallel_for_each(c.extent, [=](index<2> idx) restrict(amp) {
      int row = idx[0]; int col = idx[1];
      int sum = 0;
      for(int i = 0; i < b.extent[0]; i++)
         sum += a(row, i) * b(i, col);
      c[idx] = sum;
   });
   c.synchronize();
}

int main() {
   enum { M = 1024, N = 1024, W = 1024 };

   vector<int> vA(M * N);
   vector<int> vB(M * N);

   int i = 0;
   generate(vA.begin(), vA.end(), [&i] {return i++;});
   generate(vB.begin(), vB.end(), [&i] {return i--;});

   vector<int> vC(M * N);

   tester([&](){ calculs_v0(vA, vB, vC, M, N); }, "Calculs " + to_string(M) + "x" + to_string(N) + ", v0", cout);
   tester([&](){ calculs_v1(vA, vB, vC, M, N); }, "Calculs " + to_string(M) + "x" + to_string(N) + ", v1", cout);
   tester([&](){ calculs_v2(vA, vB, vC, M, N); }, "Calculs " + to_string(M) + "x" + to_string(N) + ", v2", cout);
   tester([&](){ calculs_v3(vA, vB, vC, M, N); }, "Calculs " + to_string(M) + "x" + to_string(N) + ", v3", cout);
   tester([&](){ calculs_v4(vA, vB, vC, M, N); }, "Calculs " + to_string(M) + "x" + to_string(N) + ", v4", cout);
   tester([&](){ calculs_v5(vA, vB, vC, M, N); }, "Calculs " + to_string(M) + "x" + to_string(N) + ", v5", cout);

   tester([&](){ MatMulSeq(vA, vB, vC, M, N, W); }, "Multiplication de matrices " + to_string(M) + "x" + to_string(N) + " et " + to_string(N) + "x" + to_string(W) + ", sequentiel", cout);
   tester([&](){ MatMulAmp(vA, vB, vC, M, N, W); }, "Multiplication de matrices " + to_string(M) + "x" + to_string(N) + " et " + to_string(N) + "x" + to_string(W) + ", AMP", cout);
}