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

int main(int argc, char** argv) {
	FILE *f;
	int dim;
	int res;
	int nprocs, rank;
	double* vector;
	int width;
	int nc;
	int i, j;
	int start, size;
	int proc, p2;
	double tmp;
	double* result;
	double* presult;
	double data;

	double time_start;
	double time_end;
	
	MPI_Status status;

	/**
	 * Init
	 */
	res = MPI_Init(&argc, &argv);
	if (res != MPI_SUCCESS) {
		printf("Error.\n");
		MPI_Abort(MPI_COMM_WORLD, res);
		return 2;
	}

	MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
	
	if (nprocs < 2) {
		if (rank == 0) {
			printf("Requires at least 2 procs.\n");
		}
		MPI_Abort(MPI_COMM_WORLD, res);
		return 1;
	}

	time_start = MPI_Wtime();

	/*
	 * Read data
	 */
	f = fopen(argv[1], "r");
	fscanf(f, "%d", &dim);

	if (rank == 0) {
		printf("Matrix dimension: %d\n", dim);
	}

	vector = (double*)malloc(sizeof(double) * dim);
	for (i=0; i<dim; i++) {
		fscanf(f, "%lf", &vector[i]);
	}

	if (rank == 0) {
		/*
		 * Distribute
		 */
		nc = nprocs - 1;
		width = (dim + nc/2 - 1) / nc;
		start = 0;

		if (width < 1) {
			width = 1;
		}

		proc = 1;
		size = width;
		while (1) {
			if (proc == nprocs-1) {
				size = dim - start;
			}
			MPI_Send(&start, 1, MPI_INT, proc, 1, MPI_COMM_WORLD);
			MPI_Send(&size, 1, MPI_INT, proc, 1, MPI_COMM_WORLD);

			start += size;
			proc++;
			p2 = proc;

			if (start >= dim) {
				break;
			}
		}

		start = -1;
		while (proc < nprocs) {
			MPI_Send(&start, 1, MPI_INT, proc, 1, MPI_COMM_WORLD);
			proc++;
		}

	}
	else {
		MPI_Recv(&start, 1, MPI_INT, 0, 1, MPI_COMM_WORLD, &status);
		if (start > -1) {
			MPI_Recv(&size, 1, MPI_INT, 0, 1, MPI_COMM_WORLD, &status);
			printf("[%d] start=%d, size=%d\n", rank, start, size);

			/*
			 * read data
			 */

			// skip first 'start' lines
			for (i=0; i<start; i++) {
				for (j=0; j<dim; j++) {
					fscanf(f, "%lf", &tmp);
				}
			}

			// read
			result = (double*)malloc(sizeof(double) * (size));
			for (i=0; i<size; i++) {
				result[i] = 0;
				for (j=0; j<dim; j++) {
					fscanf(f, "%lf", &data);
					result[i] += data * vector[j];
				}
			}

			/*
			 * send result
			 */
			MPI_Send(result, size, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD);
		}
		else {
			printf("[%d] skip\n", rank);
		}
	}

	fflush(stdout);
	MPI_Barrier(MPI_COMM_WORLD);

	/*
	 * retrieve result
	 */
	if (rank == 0) {
		result = (double*)malloc(sizeof(double) * dim);
		size = width;
		start = 0;
		presult = result;
		for (i=1; i<p2; i++) {
			if (i == p2-1) {
				size = dim - start;
			}
			MPI_Recv(presult, size, MPI_DOUBLE, i, 1, MPI_COMM_WORLD, &status);

			presult += size;
			start += size;
		}
	}

	time_end = MPI_Wtime();

	/*
	 * print result
	 */
	if (rank == 0) {
		printf("Result:\n  ");
		for (i=0; i<dim; i++) {
			printf("%5.3f ", result[i]);
		}
		printf("\n");

		printf("Execution time: %.4f\n", time_end - time_start);
	}

	MPI_Finalize();
}