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

#include "mpi.h"

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

#define DIM 512

void get_mpfproblem(MPFMatrix a, MPFVector b, MPFVector ans, long dim)
{
	long int i, j, k;
	mpf_t tmp, sqr2;

	mpf_init2(tmp, prec_mpfvector(ans));
	mpf_init2(sqr2, prec_mpfvector(ans));

	mpf_set_ui(sqr2, 2UL); mpf_sqrt(sqr2, sqr2);

	/* Frank Matrix */
	for(i = 0; i < dim; i++)
	{
		for(j = 0; j < dim; j++)
		{
			if(i < j)
			{
				mpf_set_si(tmp, dim - j);
				set_mpfmatrix_ij(a, i, j, tmp);
			}
			else
			{
				mpf_set_si(tmp, dim - i);
				set_mpfmatrix_ij(a, i, j, tmp);
			}
			mpf_mul(tmp, tmp, sqr2);
			set_mpfmatrix_ij(a, i, j, tmp);
		}
	}

	/* Answer */
	for(i = 0; i < dim; i++)
	{
		mpf_set_si(tmp, i);
		set_mpfvector_i(ans, i, tmp);
	}

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

	mpf_clear(tmp);
	mpf_clear(sqr2);
}

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;
	double start, ftime, dtime, startwtime[2], endwtime[2];

	MPFMatrix mpfa, my_mpfa[MPI_GMP_MAXPROCS];
	MPFVector mpfb, mpfx, mpfans;
	MPFVector my_mpfb, my_mpfx, my_mpfans;
	mpf_t reps, aeps;
	long int itimes_mpf, itimes_mpfm;
	double mpftime[3];

	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);

#define MPF_PREC 128

	_mpi_set_bnc_default_prec(MPF_PREC, MPI_COMM_WORLD);
	commit_mpf(&(MPI_MPF), MPF_PREC, MPI_COMM_WORLD);
	create_mpf_op(&(MPI_MPF_SUM), _mpi_mpf_add, MPI_COMM_WORLD);

	/* initialize */
	mpf_init(reps);
	mpf_init(aeps);

	/* divide problem */
	local_dim = _mpi_divide_dim(d_ddim, DIM, numprocs);
	if(myid == 0)
	{
		mpfa  = init_mpfmatrix(DIM, DIM);
		mpfb  = init_mpfvector(DIM);
		mpfx  = init_mpfvector(DIM);
		mpfans  = init_mpfvector(DIM);

		/* get problem */
		get_mpfproblem(mpfa, mpfb, mpfans, DIM);

		// print_mpfmatrix(mpfa);
	}

	/* run MPFFCG */
	mpf_set_d(reps, 1.0e-20);
	mpf_set_d(aeps, 1.0e-50);

	my_mpfb = _mpi_init_mpfvector(d_ddim, DIM, MPI_COMM_WORLD);
	my_mpfx = _mpi_init_mpfvector(d_ddim, DIM, MPI_COMM_WORLD);
	_mpi_init_mpfmatrix(my_mpfa, d_ddim, DIM, MPI_COMM_WORLD);

	_mpi_divide_mpfvector(my_mpfb, d_ddim, mpfb, MPI_COMM_WORLD);
	_mpi_divide_mpfmatrix(my_mpfa, d_ddim, mpfa, MPI_COMM_WORLD);

	if(myid == 0) startwtime[1] = MPI_Wtime();
	itimes_mpfm = _mpi_MPFCG(my_mpfx, my_mpfa, my_mpfb, reps, aeps, DIM * 5, DIM, MPI_COMM_WORLD);
	if(myid == 0) endwtime[1] = MPI_Wtime() - startwtime[1];

/*	for(i = 0; i < local_dim; i++)
	{
		printf("%5ld ", i);
		mpf_out_str(stdout, 10, 0, get_mpfvector_i(my_mpfx, i));
		printf("\n");
	}
*/
	_mpi_collect_mpfvector(mpfx, d_ddim, my_mpfx, MPI_COMM_WORLD);
//	if(myid == 0) print_mpfvector(mpfx);
	if(myid == 0)
	{
		i = 0; printf("%5d, ", i); mpf_out_str(stdout, 10, 0, gmpfvi(mpfx, i)); printf("\n");
		i = DIM/2 - 1; printf("%5d, ", i); mpf_out_str(stdout, 10, 0, gmpfvi(mpfx, i)); printf("\n");
		i = DIM - 1; printf("%5d, ", i); mpf_out_str(stdout, 10, 0, gmpfvi(mpfx, i)); printf("\n");
	}

	if(myid == 0)
	{
		free_mpfmatrix(mpfa);
		free_mpfvector(mpfb);
		free_mpfvector(mpfx);
		free_mpfvector(mpfans);
	}

	/* end */
	mpf_clear(reps); mpf_clear(aeps);

	free_mpf(&(MPI_MPF));
	free_mpf_op(&(MPI_MPF_SUM));

end:
	MPI_Finalize();

	if(myid == 0){
		/* print itimes */
		printf("Iterative Times\n");
		printf("mpf_t(MPI, %d)  : %ld(%f)\n", MPF_PREC, itimes_mpfm, endwtime[1]);
		printf("1 iter(millisec): %f milli-sec\n", 1000 * endwtime[1] / (double)itimes_mpfm);
	}

	return EXIT_SUCCESS;
}