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

#include "mpi.h"

#define USE_GMP
#define USE_MPFR
#include "mpi_bnc.h"

#define DIM 5

int main(int argc, char *argv[])
{
	long int i, j, local_dim;
	MPFMatrix c, x, b;
	MPFMatrix local_c[10], local_x[10], local_b[10];

	long int d_dim[10];
	int myrank, num_procs;

	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
	MPI_Comm_size(MPI_COMM_WORLD, &num_procs);

	_mpi_set_bnc_default_prec_decimal(50, MPI_COMM_WORLD);
	commit_mpf(&(MPI_MPF), ceil(50/log10(2.0)), MPI_COMM_WORLD);

	local_dim = _mpi_divide_dim(d_dim, DIM, num_procs);
	if(myrank == 0)
	{
		x = init_mpfmatrix(num_procs * local_dim, num_procs * local_dim);
		b = init_mpfmatrix(num_procs * local_dim, num_procs * local_dim);
		c = init_mpfmatrix(num_procs * local_dim, num_procs * local_dim);
	
		for(i = 0; i < DIM; i++)
		{
			for(j = 0; j < DIM; j++)
			{
				set_mpfmatrix_ij_d(x, i, j, (double)(i * DIM + j + 1));
				set_mpfmatrix_ij_d(b, i, j, (double)(DIM * DIM - (i * DIM + j)));
			}
		}
	}

	_mpi_init_mpfmatrix(local_c, d_dim, DIM, MPI_COMM_WORLD);
	_mpi_init_mpfmatrix(local_x, d_dim, DIM, MPI_COMM_WORLD);
	_mpi_init_mpfmatrix(local_b, d_dim, DIM, MPI_COMM_WORLD);

	_mpi_divide_mpfmatrix(local_c, d_dim, c, MPI_COMM_WORLD);
	_mpi_divide_mpfmatrix(local_x, d_dim, x, MPI_COMM_WORLD);
	_mpi_divide_mpfmatrix(local_b, d_dim, b, MPI_COMM_WORLD);

	for(i = 0; i < num_procs; i++)
		add_mpfmatrix(local_c[i], local_x[i], local_b[i]);

	_mpi_collect_mpfmatrix(c, d_dim, local_c, MPI_COMM_WORLD);

	_mpi_free_mpfmatrix(local_c, MPI_COMM_WORLD);
	_mpi_free_mpfmatrix(local_x, MPI_COMM_WORLD);
	_mpi_free_mpfmatrix(local_b, MPI_COMM_WORLD);
	
	if(myrank == 0)
	{
		print_mpfmatrix(c);

		free_mpfmatrix(x);
		free_mpfmatrix(b);
		free_mpfmatrix(c);
	}

	free_mpf(&(MPI_MPF));
	MPI_Finalize();

	return EXIT_SUCCESS;
}