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

#include "mpi.h"

#include "mpi_bnc.h"

#define DIM 512

void get_dproblem(DMatrix a, DVector b, DVector ans, long dim)
{
	long int i, j, k;
	double tmp;

	/* Frank Matrix */
	for(i = 0; i < dim; i++)
	{
		for(j = 0; j < dim; j++)
		{
			if(i < j)
				set_dmatrix_ij(a, i, j, (double)(dim - j));
			else
				set_dmatrix_ij(a, i, j, (double)(dim - i));
		}
	}

	/* Answer */
	for(i = 0; i < dim; i++)
		set_dvector_i(ans, i, (double)i);

	/* Make constant vector */
	mul_dmatrix_dvec(b, a, ans);
}

int main(int argc, char *argv[])
{
	int myid, numprocs;
	int  namelen;
	char processor_name[MPI_MAX_PROCESSOR_NAME];

	long int d_ddim[MPI_GMP_MAXPROCS], local_dim;
	DMatrix da, my_da[MPI_GMP_MAXPROCS];
	DVector db, dx, dans, my_db, my_dx, my_dans;
	double start, ftime, dtime, startwtime[2], endwtime[2];

	long int itimes_f, itimes_d, itimes_dm;
	long int i, j;

	MPI_Init(&argc, &argv);
	MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
	MPI_Comm_rank(MPI_COMM_WORLD,&myid);
	MPI_Get_processor_name(processor_name,&namelen);

	fprintf(stdout,"Process %d of %d on %s\n",
		myid, numprocs, processor_name);

	/* divide problem */
	local_dim = _mpi_divide_dim(d_ddim, DIM, numprocs);
	if(myid == 0)
	{
		/* initialize */
		da = init_dmatrix(DIM, DIM);
		db = init_dvector(DIM);
		dx = init_dvector(DIM);
		dans = init_dvector(DIM);

		/* get problem */
		get_dproblem(da, db, dans, DIM);
	
//		print_dmatrix(da);
	}


	my_db = _mpi_init_dvector(d_ddim, DIM, MPI_COMM_WORLD);
	my_dx = _mpi_init_dvector(d_ddim, DIM, MPI_COMM_WORLD);
	_mpi_init_dmatrix(my_da, d_ddim, DIM, MPI_COMM_WORLD);

	_mpi_divide_dvector(my_db, d_ddim, db, MPI_COMM_WORLD);
	_mpi_divide_dmatrix(my_da, d_ddim, da, MPI_COMM_WORLD);

	if(myid == 0) startwtime[0] = MPI_Wtime();
	itimes_dm = _mpi_DCG(my_dx, my_da, my_db, 1.0e-13, 1.0e-99, DIM * 5, DIM, MPI_COMM_WORLD);
	if(myid == 0) endwtime[0] = MPI_Wtime() - startwtime[0];
//	for(i = 0; i < local_dim; i++)
//		printf("%5ld %25.17e\n", i, get_dvector_i(my_dx, i));
	_mpi_collect_dvector(dx, d_ddim, my_dx, MPI_COMM_WORLD);
	if(myid == 0) print_dvector(dx);

	if(myid == 0)
	{

		/* run DCG */
		start = get_secv();
		itimes_d = DCG(dx, da, db, 1.0e-13, 1.0e-99, DIM * 5);
		dtime = get_secv() - start;

		/* print */
		for(i = 0; i < DIM; i++)
			printf("%5ld %25.17e %25.17e\n", i, get_dvector_i(dx, i), get_dvector_i(dans, i));
		
		/* end */
		free_dmatrix(da);
		free_dvector(db);
		free_dvector(dx);
		free_dvector(dans);
	}

	MPI_Finalize();

	if(myid == 0){
	/* print itimes */
	printf("Iterative Times\n");
	printf("double(MPI)    : %ld(%f)\n", itimes_dm, endwtime[0]);
	printf("double         : %ld(%f)\n", itimes_d, dtime);
	}
}