#include < stdio.h >
#include < stdlib.h >
#include < time.h >
#include < math.h >
#include < mpi.h >
#include < locale >
#include < iostream >

using namespace std;

int ProcNum = 0;     // Количество используемых процессов 
int ProcRank = 0;    // Ранг текущего процесса
int GridSize;        // Размер виртуальной решетки
int Coord[2];    	 // Координаты текущего процесса в решетке
MPI_Comm Grid;   	 // Коммуникатор-решетка
MPI_Comm Col;     	 // Коммуникатор столбца
MPI_Comm Row;     	 // Коммуникатор строки

void CreateCommunicators()//Создаем коммуникаторы
{
	int dims[2];
	int periods[2];//=1 для каждого измерения, являющегося периодическим
	int subdims[2];//=1 для каждого измерения в подрешетке
	periods[0]=0;
	periods[1]=0;
	dims[0]=GridSize;
	dims[1]=GridSize;
	
	MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, 2, &Grid); //создаем решетку
	//определение координат процесса в решетке
	MPI_Cart_coords(Grid, ProcRank, 2, Coord);
	//создание коммуникаторов для столбцов процессорной решетки
	subdims[0]=1;
	subdims[1]=0;
	MPI_Cart_sub(Grid, subdims, &Col);
	//создание комуникаторов для строк процессорной решетки
	subdims[0]=0;
	subdims[1]=1;
	MPI_Cart_sub(Grid, subdims, &Row);
}
void MemoryAllocation(double* &A, double* &B, double* &C, double* &BlockA, double* &BlockB, double* &BlockC, int &Size, int &BlockSize)
{	
	if (ProcRank == 0) 
	{
		do 
		{
			wcout<< L"Введите размер для квадратных матриц: ";
			wcin>> Size;
			if (Size%GridSize != 0) 
			{
				wcout<< L"Размер матриц должен быть кратен числу процессов"<i))
		{
			ATransit(BlockA, BlockSize);
		}
		if((Coord[1]!=0) && (Coord[1]>i))
		{
			BTransit(BlockB, BlockSize);
		}
	}
}

void MatrixMult(double* A, double* B, double* &C, int Size)
{
	for (int i = 0; i < Size; i++)  
		for (int j = 0; j < Size; j++)
			for (int k = 0; k < Size; k++)
				C[i*Size+j] += A[i*Size+k]*B[k*Size+j];
}


void Calculation(double* BlockA, double* BlockB, double* BlockC, int BlockSize)
{
	 for (int i = 0; i < GridSize; i++)
	 {
		
		MatrixMult(BlockA, BlockB, BlockC, BlockSize);
		ATransit(BlockA, BlockSize);
		BTransit(BlockB,BlockSize);
		
	}
}
void Result(double* C, double* BlockC, int Size, int BlockSize)
{
	double * Buff = new double [Size*BlockSize];
	for (int i = 0; i < BlockSize; i++) 
	{
		 MPI_Gather( &BlockC[i*BlockSize], BlockSize, MPI_DOUBLE, &Buff[i*Size], BlockSize, MPI_DOUBLE, 0, Row);//сборка строк матрицы С из блоков
	}
	if (Coord[1] == 0) 
	{
		 MPI_Gather(Buff, BlockSize*Size, MPI_DOUBLE, C, BlockSize*Size, MPI_DOUBLE, 0, Col);//сборка матрицы С из строк
	}
	delete [] Buff;

}

// Вывод матрицы
void Print (double* Matrix, int Size) 
{
 for (int i = 0; i < Size; i++) 
 {
    for (int j = 0; j < Size; j++)
		if ( j == Size-1 ) 
			cout << Matrix[i*Size+j] << "\n" ;
		else cout << Matrix[i*Size+j] << " " ;
 }
}
void Clear(double* A, double* B, double* C, double* BlockA, double* BlockB, double* BlockC)
{
	if (ProcRank == 0) 
	{
		delete [] A; 
		delete [] B;
		delete [] C; 
	}	
	delete [] BlockA;
	delete [] BlockB;
	delete [] BlockC;
}

void main(int argc, char* argv[]) 
{
  double* A;  // Матрица А
  double* B;  // Матрица В
  double* C;  // Результат
  int Size;          // Размер матриц
  int BlockSize;     // Размер блоков
  double* BlockA;   // блок матрицы А
  double* BlockB;   // блок матрицы В
  double* BlockC;   // блок результирующей матрицы
  clock_t Start, med, End;
  double Finish;

  locale loc("rus_rus.866");
  locale::global(loc);
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &ProcNum);
  MPI_Comm_rank(MPI_COMM_WORLD, &ProcRank);
  GridSize = sqrt((double)ProcNum);
  if (ProcNum != GridSize*GridSize) 
  {
    if (ProcRank == 0) 
	{
      wcout << L"Ошибка! Число процессов должно быть полным квадратом" << endl;
    }
  }
  else 
  {
	  if (ProcRank == 0)
	  {
		wcout << L"Алгорит Кэннона умножения матриц при блочном разделении данных" << endl;
	  }
	  CreateCommunicators();	// Создаем коммуникаторы 
          MemoryAllocation( A, B, C, BlockA, BlockB, BlockC, Size, BlockSize );	// выделение памяти и инициализация элементов массивов
	  if (ProcRank == 0) Start=clock();// замер времени
	  DataDistribution(A, B, BlockB, BlockA, Size, BlockSize);//разделение данных
	  MPI_Barrier(Grid);// синхронизация
          Calculation(BlockA, BlockB, BlockC, BlockSize);// вычисления
	  Result(C, BlockC, Size, BlockSize);//сбор результата
	  MPI_Barrier(Grid);// синхронизация
	  if (ProcRank == 0)
	  {
		  End=clock();
		  Finish=(double)(End-Start)/ CLOCKS_PER_SEC;
	  }
      if (ProcRank == 0) 
	  {
			if(Size < 8)
			{
				wcout << L"Матрица А:" << endl;
				Print(A, Size);
				wcout << L"Матрица В:" << endl;
				Print(B, Size);
				wcout << L"Матрица С:" << endl;
				Print(C, Size);
			}
			wcout << L"Время с рассылкой блоков: " << Finish << endl;
	  }
      Clear (A, B, C, BlockA, BlockB, BlockC); // освобождение памяти
  }
  MPI_Finalize();
}