Lähdetään etsimään motivaatiota kurssin suorittamiseen seuraavasta ohjelmanpätkästä:
b
ja laskee kullekkin b
:n alkiolle käänteisluvun käyttäen neljää Newtonin iteraatiota. Lopputulokset tallennetaan a
-taulukkoon.void process(double *a, const double *b, int n) {
#pragma omp parallel for
for(int i = 0; i < n; i++) {
const double md = -b[i];
double y = (48.0/17.0) + (32.0/17.0) * md;
y = y + y * (1.0 + md * y);
y = y + y * (1.0 + md * y);
y = y + y * (1.0 + md * y);
y = y + y * (1.0 + md * y);
a[i] = y;
}
}
process
-aliohjelmaa kutsutaan seuraavasti:#define N (32*1024*1024)
// .......................
double *a = new double[N];
double *b = new double[N];
for(int i = 0; i < N; i++)
b[i] = 0.5 + 0.5*i/N;
process(a, b, N);
// .......................
delete [] a;
delete [] b;
Time: 0.0456872 s
Flops: 13.9543 GFlops
Max error: 2.22045e-16
Tarkastellaan seuraavaksi CUDA:lla toteutettua versiota samasta numeerisesta algoritmistä. Ohjelmakoodin yksityiskohtia ei tarvitse vielä tässä vaiheessa ymmärtää sillä niihin palataan myöhemmin:
threadId
ja säikeiden määrä tallennetaan threadCount
-muuttujaan__global__ void process(double *a, const double *b, int n) {
int threadId = blockIdx.x * blockDim.x + threadIdx.x;
int threadCount = gridDim.x * blockDim.x;
for(int i = threadId; i < n; i += threadCount) {
const double md = -b[i];
double y = (48.0/17.0) + (32.0/17.0) * md;
y = y + y * (1.0 + md * y);
y = y + y * (1.0 + md * y);
y = y + y * (1.0 + md * y);
y = y + y * (1.0 + md * y);
a[i] = y;
}
}
double *a = new double[N];
double *b = new double[N];
...
double *d_a, *d_b;
cudaMalloc((void **)&d_a, N*sizeof(double));
cudaMalloc((void **)&d_b, N*sizeof(double));
cudaMemcpy(d_b, b, N*sizeof(double), cudaMemcpyHostToDevice);
...
// Käynnistetään noin miljoona säiettä GPU:lla
process<<<1024, 1024>>>(d_a, d_b, N);
...
cudaMemcpy(a, d_a, N*sizeof(double), cudaMemcpyDeviceToHost);
Time: 0.00306797 s
Flops: 207.803 GFlops
Max err: 2.22045e-16
D\I | Single Instruction | Multiple Instruction |
---|---|---|
Single Data | SISD | MISD |
Multiple Data | SIMD | MIMD |
#pragma omp
-direktiividouble y = x + 7.0;
for(int i = 0; i < N; i++)
y[i] = 4.0 * x[i] + 7.0;
for(int i = 0; i < N; i++) {
x[i] = x[i] + 5.0;
for(int j = 0; j < M; j++)
y[i] = A[i][j] * x[i] + 7.8;
}
\[ v = \frac{W}{T} \]
\[ Q = \rho vA = \frac{\rho AW}{T} = \frac{\rho A}{L} \]
#define N 10000000
...
void axpy(double *y, const double *x, double a, int n) {
#pragma omp parallel for
for(int i = 0; i < n; i++)
y[i] = a * x[i] + y[i];
}
...
axpy(x, y, a, N);
$ for((i=1;i<=12;i++)); do OMP_NUM_THREADS=$i ./strong; done
Suoritusaika 1 säikeellä: 0.0114419
Suoritusaika 2 säikeellä: 0.00692391
Suoritusaika 3 säikeellä: 0.00632405
Suoritusaika 4 säikeellä: 0.006181
...
#define N 10000000
...
void axpy(double *y, const double *x, double a, int n) {
#pragma omp parallel for
for(int i = 0; i < n; i++)
y[i] = a * x[i] + y[i];
}
...
int n = threadCount*N; // Säikeen tekemä työ pysyy vakiona
...
axpy(x, y, a, n);
$ for((i=1;i<=12;i++)); do OMP_NUM_THREADS=$i ./weak; done
Suoritusaika 1 säikeellä: 0.0115819
Suoritusaika 2 säikeellä: 0.013674
Suoritusaika 3 säikeellä: 0.0196369
Suoritusaika 4 säikeellä: 0.024848
...