#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys/timeb.h>
#include <mpi.h>
#include <string.h>
#include <unistd.h>
#include <mkl.h>
#include <mkl_blas.h>
#include <mkl_cblas.h>
#include <mkl_scalapack.h>
#include <mkl_pblas.h>
#include <mkl_blacs.h>
#include <mkl_lapacke.h>
#include <math.h>
#include <ctype.h>
#include "matrices.h"
#include "common.h"

void GenerateRandomMatrix(real* M, size_t n, size_t m)
{
	register size_t i;
	for (i = 0; i < n * m; i++)
		M[i] = rand() / ((real)RAND_MAX);
}


int main(int argc, const char** argv){



	size_t N = 1000;
	size_t M = 1000;
	srand(0);

	if (argc > 1)
	{
		if (argc < 3)
		{
			fprintf(stderr, "Not enough input arguments.\n");
			printf("Usage: AtA [N M]\n");
			exit(1);
		}

		N = atol(argv[1]);
		M = atol(argv[2]);

		if (argc > 3)
			srand(atoi(argv[3]));
	}


	MPI_Init(NULL, NULL);
	int id, size;

	MPI_Comm_rank(MPI_COMM_WORLD, &id);
	MPI_Comm_size(MPI_COMM_WORLD, &size);

	FILE *deb;
	char s[16];
	sprintf(s, "log%d.txt", id);
	deb = fopen(s, "w");


	double t1, t2, diff;

	printf("id %d\n", id);
	printf("Generating random %lu-by-%lu matrix... ", N, M);
	t1 = MPI_Wtime();
	real* A = (real*)malloc(N * M * sizeof(real));
	if (A == NULL)
	{
		fprintf(stderr, "Cannot allocate A->");
		exit(1);
	}
	if (id == 0){
		GenerateRandomMatrix(A, N, M);
	} else{
		free(A);
		A = NULL;
	}
	real *C = (real*)calloc(M * M, sizeof(real));

	// long long bhandle = sys2blacs_handle(MPI_COMM_WORLD);
 //    long long icontxt = bhandle;
	long long icontxt;
	long long what = 0;
    long long npcol;// 	qui npcol = 1;
    long long nprow;//	qui nprow = size;
    long long myrow, mycol;
    int dims[] = {0, 0};

    MPI_Dims_create(size, 2, dims);
    nprow = dims[0];
    npcol = dims[1];
//    fprintf(deb, "grid is %d x %d\n", nprow, npcol);
    long long negone = -1;
    long long one = 1;
    long long zero = 0;
    long long n = (long long)N;
    long long m = (long long)M;
    MKL_INT descA[9], descC[9], descSubA[9], descSubC[9];
    


    blacs_get(&negone, &what, &icontxt);
    blacs_gridinit( &icontxt, "R", &nprow, &npcol );
    blacs_gridinfo( &icontxt, &nprow, &npcol, &myrow, &mycol );
    MKL_INT np, mp, nb = 1, mb = 1;

    np = numroc_( &n, &nb, &myrow, &zero, &nprow );
    mp = numroc_( &m, &mb, &mycol, &zero, &npcol );

//    printf("I am proc %d, the grid is %dx%d\n my coordinates are (%d,%d)\nnb = %d, mb= %d, np = %d, mp = %d\n", id, nprow, npcol, myrow, mycol, nb, mb, np, mp);

    long long lldA = n;
    long long lldC = m;
    long long lldSubA = np;
    long long lldSubC = mp;
    MKL_INT info;
    real dzero = 0;
    real done = 1;
    real *subA = (real*)calloc(np * mp, sizeof(real));
    real *subC = (real*)malloc(mp * mp * sizeof(real));


	descinit_(descA, &n, &m, &n, &m, &zero, &zero, &icontxt, &lldA, &info);    
	descinit_(descC, &m, &m, &m, &m, &zero, &zero, &icontxt, &lldC, &info); 

        descinit_(descSubA, &n, &m, &nb, &mb, &zero, &zero, &icontxt, &lldSubA, &info);    
	descinit_(descSubC, &m, &m, &mb, &mb, &zero, &zero, &icontxt, &lldSubC, &info); 

	pdgeadd_("N", &n, &m, &done, A, &one, &one, descA, &dzero, subA, &one, &one, descSubA);
	t1 = MPI_Wtime();
	pdsyrk_("U", "N", &m, &n, &done, subA, &one, &one, descSubA, &dzero, subC, &one, &one, descSubC);

	//fprintf(deb,"C before pdgeadd\n");
	//fprint_mat(C, m, m, m, m, deb);

	pdgeadd_("N", &m, &m, &done, subC, &one, &one, descSubC, &dzero, C, &one, &one, descC);
	//fprintf(deb,"\nC after pdgeadd\n");
	//fprint_mat(C, m, m, m, m, deb);
	
	//MPI_Bcast(C, m * m, MPI_RREAL, 0, MPI_COMM_WORLD);
	//fprintf(deb,"\nC after MPI_Bcast\n");
	//fprint_mat(C, m, m, m, m, deb);
	t2 = MPI_Wtime();
	diff = t2 - t1;
	printf("%1.2f seconds id %d.\n", diff, id);
	// if (id == 0){
	// 	fprint_mat(C, m, m, m, m, deb);
	// 	fprintf(deb, "\n");
	// 	memset(C, 0, m * m * sizeof(real));
	// 	cblas_dsyrk(CblasRowMajor, CblasLower, CblasTrans,
	// 		   m, n, done, A, m, 0, C, m);
	// 	fprint_mat(C, m, m, m, m, deb);
	//}

    blacs_freebuff(&icontxt, &zero);
    blacs_gridexit(&icontxt);
    blacs_exit(&zero);
    fclose(deb);

    return 0;
}