并行算法之前缀和

#include<stdio.h>
#include<stdlib.h>
#include<fcntl.h>
#include<pthread.h>
#include<math.h>
#include<time.h>

int threadnum = 1025;
int n = 0;//n is less than threadnum
int A[1025], B[1025][1025], C[1025][1025];

struct param
{
	int i;
	int j;
	int value;
}param;

int input()
{
	int i = 0, j = 0;
	//printf("input the element number of array A: ");
	scanf("%d", &n);
	printf("size = %d\n", n);
	if(n == 0) return 0;
	A[0] = 0;
	//printf("init the array of A:\n");
	for(i=0; i<n; i++)
		scanf("%d", &A[i]);
	return 1;
}

void output()
{
	int i, j;
	printf("matrix B:\n");
	for(i=0; i<=(int)(log(n)/log(2)); i++)
	{
		for(j=0; j<(int)(n/pow(2, i)); j++)
			printf("%3d ", B[i][j]);
		printf("\n");
	}
	printf("matrix C\n");
	for(i=0; i<=(int)(log(n)/log(2)); i++)
	{
		for(j=0; j<(int)(n/pow(2, i)); j++)
			printf("%3d ", C[i][j]);
		printf("\n");
	}
}

void *initMatrixB(void *arg)
{	
	struct param *data = (struct param *)arg;
	int j = data->j;
	B[0][j] = data->value;
	//printf("B[0][%d] = %d ", j, B[0][j]);
	free(data);
	data = NULL;
	pthread_exit(0);
}


void *caculateB(void *arg)
{
	struct param *data = (struct param *)arg;
	int h = data->i;
	int j = data->j;
	B[h][j] = B[h-1][2*j] + B[h-1][2*j+1];
	//printf("B[%d][%d]=%d ", h, j, B[h][j]);
	free(data);
	data = NULL;
	pthread_exit(0);
}

void *caculateC(void *arg)
{
	struct param *data = (struct param *)arg;
	int h = data->i;
	int j = data->j;
	if(j == 0)
	{
		C[h][0] = B[h][0];
		//printf("C[%d][%d]=%d ", h, j, C[h][j]);
	}
	else if(j%2 != 0)
	{
		C[h][j] = C[h+1][j/2];
		//printf("C[%d][%d]=%d ", h, j, C[h][j]);
	}
	else
	{
		C[h][j] = C[h+1][(j-1)/2] + B[h][j];
		//printf("C[%d][%d]=%d ", h, j, C[h][j]);
	}
	free(data);
	data = NULL;
	pthread_exit(0);
}

int main(int argc, char **argv)
{
	clock_t start, end;
	freopen("data.dat", "r", stdin);
	while(input())
	{
		int i, j, h;
		start = clock();
		pthread_t threadids[threadnum];
		//printf("init matrix B\n");
		for(i=0; i<n; i++)
		{
			struct param *data = (struct param *)malloc(sizeof(struct param));
			data->i = 0;
			data->j = i;
			data->value = A[i];
			pthread_create(&threadids[i], NULL, initMatrixB, (void *)data);
		}
		for(i=0; i<n; i++)
		{
			pthread_join(threadids[i], NULL);
		}
		
		//printf("caculate matrix B\n");
		for(h=1; h<=(int)(log(n*1.0)/log(2*1.0)); h++)
		{
			for(j=0; j<(int)(n/pow(2, h)); j++)
			{
				struct param *data = (struct param *)malloc(sizeof(struct param));
				data->i = h;
				data->j = j;
				data->value = 0;
				pthread_create(&threadids[j], NULL, caculateB, (void *)data);
			}
			for(i=0; i<(int)(n/pow(2, h)); i++)
			{
				pthread_join(threadids[i], NULL);
			}
		}
		
		//printf("caculate matrix C\n");
		for(h=(int)(log(n*1.0)/log(2*1.0)); h>=0; h--)
		{
			for(j=0; j<(int)(n/pow(2*1.0, h*1.0)); j++)
			{
				struct param *data = (struct param *)malloc(sizeof(struct param));
				data->i = h;
				data->j = j;
				data->value = 0;
				pthread_create(&threadids[j], NULL, caculateC, (void *)data);
			}
			for(i=0; i<(int)(n/pow(2*1.0, h*1.0)); i++)
			{
				pthread_join(threadids[i], NULL);
			}
		}
		end = clock();
		//output();
		printf("the non-Recursive prefix sum algorithm runs %f seconds\n", (double)(end-start)/CLOCKS_PER_SEC);
	}
	return 0;
}