並列計算コース復習ノート

1
2
3

1.ある直列アプリケーションで、20%の割合が直列実行される必要があります。現在、3倍の性能向上を実現する必要があります。この目標を達成するためには、何個のCPUが必要ですか？5倍の加速比を実現するためには、何個のCPUが必要ですか？
2.5台のコンピュータで実行されている並列プログラムで、10%の並列部分があります。1台のコンピュータでの直列実行に対して、加速比はどれくらいですか？加速比を2倍にしたい場合、何個のCPUが必要ですか？
3.並列化できない部分が5%のアプリケーションを並列プログラムに変更します。現在、市場には2種類の並列コンピュータがあります：コンピュータXは4個のCPUを持ち、各CPUは1時間以内にこのアプリケーションを実行できます。コンピュータYは16個のCPUを持ち、各CPUは2時間以内にこのアプリケーションを実行できます。実行時間を最小化する必要がある場合、どのコンピュータを購入すべきですか？

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

void * MTFlipHM(void * tid){
        struct Pixel pix; //temp swap pixel
        int row, col;
        int id = *((int *) tid);
        int start = id * ip.Vpixels/NumThreads;
        int end = start + ip.Vpixels/NumThreads;
        unsigned char buffer[16384];
        for (row = st ; row < ed; row++)
        {
            memcpy(buffer, TheImage[row], ip.Hbytes);
            col = 0;
            while (col < ip.Hpixels * 3 /2){
            pix.B = buffer[col];
            pix.G = buffer[col+1];
            pix.R = buffer[col+2];
            buffer[col]   = buffer[ip.Hpixels*3 - col -3];
            buffer[col+1] = buffer[ip.Hpixels*3 - col -2];
            buffer[col+2] = buffer[ip.Hpixels*3 - col -1];
            buffer[ip.Hpixels*3 - col -3] = pix.B;
            buffer[ip.Hpixels*3 - col -2] = pix.G;
            buffer[ip.Hpixels*3 - col -1] = pix.R;
            col += 3;
            }
        }
        memcpy(TheImage[row],buffer,ip.Hbytes);
        pthread_exit(NULL);
}

void * MTFlipVM(void * tid){
        struct Pixel pix; //temp swap pixel
        int row, col;
        int id = *((int *) tid);
        int start = id * ip.Vpixels/NumThreads;
        int end = start + ip.Vpixels/NumThreads;
        unsigned char buffer1[16384], buffer2[16384];
        for (row = start; row < end; row ++)
        {
            memcpy(buffer1,TheImage[row],ip.Hbytes);
            int  mirrow = ip.Vpixels - 1 - row;
            memcpy(buffer2,TheImage[mirrow],ip.Hbytes);
            // 再错位拷贝即完成交换
            memcpy(TheImage[row],buffer2,ip.Hbytes);
            memcpy(TheImage[mirrow],buffer1,ip.Hbytes);
        }
}

pthread_create(&ThHandle[i], &ThAttr, MTFlipFunc, (void *)&ThParam[i]);
for(i=0; i<NumThreads; i++)
                pthread_join(ThHandle[i], NULL);

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

const int a=1, b=2, c=3;
__global__ void add(double *x, double *y, double* z){
	const int n = blockIdx.x * blockDim.x + threadIdx.x;
	if (n<N)  z[n] = x[n] + y[n];
}
int main(){
	const int N = 1e9;
	const int M = sizeof(double) * N;
	double *h_x = (double*) malloc(M);
	double *h_y = (double*) malloc(M);
	double *h_z = (double*) malloc(M);
	for (int i=0; i<N; i++)
	{
		h_x[i] = a;
		h_y[i] = b;
	}
	double *d_x, *d_y, *d_z;
	cudaMalloc((void**) &d_x, M);
	cudaMalloc((void**) &d_y, M);
	cudaMalloc((void**) &d_z, M);

	const int block_size =  128;
	int grid_size = (N+block_size-1) / block_size;
	add<<<grid_size, block_size>>>(d_x, d_y, d_z);
	cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost);
	free(h_x),free(h_y),free(h_z);
	cudaFree(d_x),cudaFree(d_y),cudaFree(d_z);
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

__global__ void Vflip(uch *ImgDst, uch *ImgSrc, ui Hpixels, ui Vpixels){
ui ThrPerBlk = blockDim.x;
ui MYbid = blockIdx.x;
ui MYtid = threadIdx.x;
ui MYgtid = ThrPerBlk * MYbid + MYtid;
ui BlkPerRow = (Hpixels + ThrPerBlk - 1) / ThrPerBlk; // ceil
ui RowBytes = (Hpixels * 3 + 3) & (~3);
ui MYrow = MYbid / BlkPerRow;
ui MYcol = MYgtid - MYrow*BlkPerRow*ThrPerBlk;
if (MYcol >= Hpixels) return;// col out of range
ui MYmirrorrow = Vpixels - 1 - MYrow;
ui MYsrcOffset = MYrow * RowBytes;
ui MYdstOffset = MYmirrorrow * RowBytes;
ui MYsrcIndex = MYsrcOffset + 3 * MYcol;
ui MYdstIndex = MYdstOffset + 3 * MYcol;
// swap pixels RGB @MYcol , @MYmirrorcol
ImgDst[MYdstIndex] = ImgSrc[MYsrcIndex];
ImgDst[MYdstIndex + 1] = ImgSrc[MYsrcIndex + 1];
ImgDst[MYdstIndex + 2] = ImgSrc[MYsrcIndex + 2];}

__global__ void Hflip(uch *ImgDst, uch *ImgSrc, ui Hpixels){
ui ThrPerBlk = blockDim.x;
ui MYbid = blockIdx.x;
ui MYtid = threadIdx.x;
ui MYgtid = ThrPerBlk * MYbid + MYtid;
ui BlkPerRow = (Hpixels + ThrPerBlk -1 ) / ThrPerBlk; // ceil
ui RowBytes = (Hpixels * 3 + 3) & (~3);
ui MYrow = MYbid / BlkPerRow;
ui MYcol = MYgtid - MYrow*BlkPerRow*ThrPerBlk;
if (MYcol >= Hpixels) return;// col out of range
ui MYmirrorcol = Hpixels - 1 - MYcol;
ui MYoffset = MYrow * RowBytes;
ui MYsrcIndex = MYoffset + 3 * MYcol;
ui MYdstIndex = MYoffset + 3 * MYmirrorcol;
// swap pixels RGB @MYcol , @MYmirrorcol
ImgDst[MYdstIndex] = ImgSrc[MYsrcIndex];
ImgDst[MYdstIndex + 1] = ImgSrc[MYsrcIndex + 1];
ImgDst[MYdstIndex + 2] = ImgSrc[MYsrcIndex + 2];}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

__global__ void transpose(int a[], int b[], int N){
    //分配共享内存
    __shared__ int S[TILE][TILE + 1];
    int bx = blockIdx.x * TILE;
    int by = blockIdx.y * TILE;
    int ix = bx + threadIdx.x;
    int iy = by + threadIdx.y;
    if (ix < N && iy < N)// 读入共享内存
        S[threadIdx.y][threadIdx.x] = a[iy * N + ix];
    __syncthreads();//同步，这是必不可少的
    int ix2 = bx + threadIdx.y;
    int iy2 = by + threadIdx.x;
    if (ix2 <N && iy2 < N)// 写回
        b[ix2 * N + iy2 ] = S[threadIdx.x][threadIdx.y];
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

__shared__ float Mds[WIDTH][TILE_WIDTH];
__shared__ float Nds[TILE_WIDTH][WIDTH];
int bx=blockIdx.x ; int by = blockIdx.y;
int tx=threadIdx.x ; int ty = threadIdx.y;
int Row = by *TILE_WIDTH +ty;
int Col = bx*TILE_WIDTH + tx;
float Pvalue = 0;
for(int m=0; m<WIDTH/TILE_WIDTH; ++m)
    {
    // 每个线程载入M的子矩阵的一个元素
    Mds[ty][tx] = Md[Row*width+(m*TILE_WIDTH+tx)];
    //每个线程载入N的子矩阵的一个元素
    Nds[ty][tx] = Nd[(m*TILE_WIDTH+ty)*width+Col];
	 __syncthreads();
	 for (int k = 0; k < TILE_WIDTH; ++k)
	    Pvalue += Mds[ty][k] * Nds[k][tx];
	 __syncthreads();
	 }
Pd[Row*WIDTH+Col] = Pvalue;//将结果写回P矩阵

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

#define SIZE (100*1024*1024)
//通过工具函数big_random_block()来生成随机的字节流
unsigned char *buffer =(unsigned char*)big_random_block( SIZE );
unsigned int histo[256];
for (int i = 0; i<256; i++)
	histo[i] = 0;
for (int i = 0; i < SIZE; i++)
	histo[buffer[i]]++;
long histoCount = 0;
for (int i = 0; i<256; i++) {
	histoCount += histo[i]; }

__global__ void histo_kernel(unsigned char *buffer, long size, unsigned int *histo){
__shared__ unsigned int temp[256];
temp[threadIdx.x] = 0;
__syncthreads();
int i = threadIdx.x + blockIdx.x * blockDim.x;
int offset = blockDim.x *gridDim.x;
while (i<size){
	atomicAdd(&temp[buffer[i]], 1);
	i += offset;
}
__syncthreads();
atomicAdd(&(histo[threadIdx.x]), temp[threadIdx.x]);
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

__global__ void _sum_gpu(int *input, int count, int *output)
{
    __shared__ int sum_per_block[BLOCK_SIZE];

    int temp = 0;
    for (int idx = threadIdx.x + blockDim.x * blockIdx.x;
         idx < count; idx += gridDim.x * blockDim.x
	)
    {// 跨网格循环，一个线程加多个数据，应对海量数据
        temp += input[idx];
    }

    sum_per_block[threadIdx.x] = temp;  //the per-thread partial sum is temp!
    __syncthreads();

    //**********shared memory summation stage***********
    for (int length = BLOCK_SIZE / 2; length >= 1; length /= 2)
    {
        int double_kill = -1;
	if (threadIdx.x < length)
	{
	    double_kill = sum_per_block[threadIdx.x] + sum_per_block[threadIdx.x + length];
	}
	__syncthreads();  //why we need two __syncthreads() here, and,
	
	if (threadIdx.x < length)
	{
	    sum_per_block[threadIdx.x] = double_kill;
	}
	__syncthreads();  //....here ?
	
    } //the per-block partial sum is sum_per_block[0]

    if (blockDim.x * blockIdx.x < count) //in case that our users are naughty
    {
        //the final reduction performed by atomicAdd()
        if (threadIdx.x == 0) atomicAdd(output, sum_per_block[0]);
    }
}

1
2
3
4
5

typedef struct {     float value;     int index; } Tuple;  
Tuple *d_tuples; 
int threadsPerBlock = 256; 
int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
initializeTuples<<<blocksPerGrid, threadsPerBlock>>>(d_data, d_tuples, n);

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

int *d_indices;
kReduceKernel<<<blocksPerGrid, threadsPerBlock>>>(d_tuples, d_indices, n, k);
__global__ void kReduceKernel(Tuple *input, int *output, int n, int k) {
    extern __shared__ Tuple shared[];
    int tid = threadIdx.x;
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    shared[tid] = (i < n) ? input[i] : Tuple{0, 0};
    __syncthreads();
    for (int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s)
            shared[tid] = (shared[tid].value > shared[tid + s].value) ? shared[tid] : shared[tid + s];
        __syncthreads();
    }

    if (tid == 0)
        output[blockIdx.x] = shared[0].index;
}

1
2
3
4
5

cudaMemcpy(h_indices, d_indices, size, cudaMemcpyDeviceToHost);  
for (int i = 0; i < k; ++i) {     
    int index = h_indices[i];     
    h_result[i] = h_data[index]; }  
std::sort(h_result, h_result + k);