#include <stdio.h>
#include <stdlib.h>
#include "mpi.h"

#define BLOCK_LOW(id,p,n)  ((id)*(n)/(p))
#define BLOCK_HIGH(id,p,n) (BLOCK_LOW((id)+1,p,n)-1)
#define BLOCK_SIZE(id,p,n) (BLOCK_HIGH(id,p,n)-BLOCK_LOW(id,p,n)+1)

double* allocate1D(int n) 
{
	double *a;
	a = (double *)malloc(n * sizeof(double));
	return a;
}

double** allocate2D(int m, int  n) 
{
	double **a;
	double *a_storage;
	int i;

	a_storage = (double*)malloc(m*n*sizeof(double));
	a = (double**)malloc((unsigned)m*sizeof(double*));

	for(i=0; i<m; i++)
	{
		a[i] = &a_storage[i*n];
	}

	return a;
}

void printArray(double *a, int n) 
{
	FILE *file = fopen("outputCol.dat","w");
	int i;

	for(i = 0; i < n; i++) 
	{
		fprintf(file,"%1.2f\n", a[i]);
	}

	fclose(file);
}

void printMatrix(double **a, int m, int n) 
{
	int i, j;

	for(i=0; i<m; i++) 
	{
		for(j=0; j<n; j++) 
		{ 
			printf("%1.2f\t", a[i][j]);
		}
		printf("\n");
	}
	  
	printf("\n");
}

int main(int argc, char *argv[]) 
{
	int m, n, nv;
	int id, np;
	int i, j, r;
	int *counts, *displs;
	double **local_A;
	double *local_b, *c, *local_c, *row;	
	double elapsed_time;
	MPI_Status status;
	
	MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &id);
    MPI_Comm_size(MPI_COMM_WORLD, &np);
	
	counts = (int *)malloc((unsigned)np * sizeof(int));	
	displs = (int *)malloc((unsigned)np * sizeof(int));	
	
	if(id == 0)	
	{
		scanf("%d%d%d",&m, &n, &nv);
		
		if(n != nv) 
		{
			printf("Broj kolona matrice se ne poklapa sa brojem elemenata vektora\n");
			MPI_Abort(MPI_COMM_WORLD,1);
			MPI_Finalize();
		}
	}		
	
	MPI_Bcast(&m, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&nv, 1, MPI_INT, 0, MPI_COMM_WORLD);
		
	counts[0] = BLOCK_SIZE(0,np,n);		
	displs[0] = 0;
	
	for(i = 1; i < np; i++)	
	{
		displs[i] = displs[i-1] + (counts[i-1]);
		counts[i] = BLOCK_SIZE(i,np,n);
	}	
	
	local_A = allocate2D(m,counts[id]); 	
	local_b = allocate1D(counts[id]); 
	
	if(id == 0)	
	{		
		row = allocate1D(n);
	}
		
	for (i = 0; i < m; i++)	
	{
		if(id == 0)	
		{
			for(j = 0; j < n; j++) 
			{
				scanf("%lf", &row[j]);
			}	
		}		
		MPI_Scatterv (row, counts, displs, MPI_DOUBLE, &(local_A[i][0]), counts[id], MPI_DOUBLE, 0, MPI_COMM_WORLD);
	}
	
	if(id == 0)	
	{			
		for(j = 0; j < n; j++) 
		{
			scanf("%lf", &row[j]);	
		}
	}
	
	MPI_Scatterv (row, counts, displs, MPI_DOUBLE, local_b, counts[id], MPI_DOUBLE, 0, MPI_COMM_WORLD);
				
	MPI_Barrier(MPI_COMM_WORLD);
		
	if(id == 0)	
	{
		free(row);				
		c = allocate1D(m);		
		elapsed_time = -MPI_Wtime();   
	}
		
	local_c = allocate1D(m); 	
	
	for(i = 0; i < m; i++) 
	{
		local_c[i] = 0.0;
		for(j = 0; j < counts[id]; j++) 
		{
			local_c[i] += local_A[i][j] * local_b[j];
		}
	
		MPI_Reduce(local_c, c, m, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); 
	}
	
	if(id == 0)	
	{	
		elapsed_time += MPI_Wtime();
		
		printf("Proteklo vreme : %lf,  procesora %d, m = %d, n = %d\n", elapsed_time, np, m, n);		
		printArray(c, m);
		
		free(c);
	}	
		
	free(*local_A);
	free(local_b);
	free(local_c);
	
	free(counts);
	free(displs);
	
	MPI_Finalize();
	return 0;
}