#include "stdafx.h"
#include  <  stdio.h  > 
#include  <  stdlib.h  > 
#include  <  conio.h  > 
#include  <  math.h  > 
#include  <  mpi.h  > 

float* pMatrix;  		// Матрица линейной системы
float* pVector;  		// Вектор правых частей линейной системы
int* pProcNum;			// Номера строк, обрабатываемых данным процессом
float* pProcRows; 		// Строки матрицы A
float* pProcVector; 	// Блок вектора b
int size;       		// Размер матрицы и векторов
int RowNum;     		// Количество строк матрицы 
double start, finish, duration;
int ProcNum;
int ProcRank;

void ProcessInitialization()
{
	if (ProcRank == 0) 	
	{
		int show = 0, fill;
		printf("\nEnter size of matrix: ");
		scanf("%i", &size);
		printf("Do you want to fill matrix manually(0) or automatically(1)?: ");
		scanf("%i", &fill);		
		if(fill)
		{
			printf("Do you want to see matrix(1) or not(0)?: ");
			scanf("%i", &show);
		}

		pMatrix = new float[size*size];
		pVector = new float[size];

		//автоматическое заполнение матрицы
		if(fill==1)
		{
			//int d = 1000;
			srand(1);
			for(int i=0; i  <  size*size; i++)
				pMatrix[i]=(float)(rand() % 1000);
			for(int i=0; i  <  size; i++)
				pVector[i]=(float)(rand() % 1000);
			if(show)
			{
				for(int i=0; i < size*size; i++)				 
					printf("\npMatrix[%i]=%f", i, pMatrix[i]);
				printf("\n");
				for(int i=0; i < size; i++)				
					printf("\npVector[%i]=%f", i, pVector[i]);
				printf("\n");
			}
		}
		//ручное заполнение матрицы
		else
		{
			printf("\nEnter values of pMatrix:\n");
			for(int i=0; i < size*size; i++)
			{
				printf("pMatrix[%i] = ", i);
				scanf("%f",&pMatrix[i]);
			}
			printf("\nEnter values of pVector:\n");
			for(int i=0; i < size; i++)
			{
				printf("pVector[%i] = ", i);
				scanf("%f",&pVector[i]);
			}
		}
	}

	//рассылка размера матрицы pMatrix на все процессы 0-м процессом
	MPI_Bcast(&size,1,MPI_INT,0,MPI_COMM_WORLD);	
	
	int StringNum;
	RowNum = 0;
	
	//определение числа строк, обрабатываемых данным процессом
	for(int j=0;j  <  =size/ProcNum;j++)
	{
		StringNum = ProcRank + j*ProcNum;
		if(StringNum  >  =size)
			break;
		RowNum++;		
	}
	pProcRows = new float[RowNum*size];
	pProcVector = new float[RowNum];
	pProcNum = new int[RowNum];	
}

void DataDistribution()
{
	//заполнение массива номеров строк pProcNum
	int StringNum;
	MPI_Status status;
	for(int j=0;j < RowNum;j++)
		pProcNum[j]=ProcRank + j*ProcNum;	
	
	//рассылка 0-м процессом остальным процессам значений матрицы pMatrix
	if(ProcRank==0)
	{
		int rank;
		for(int i=0;i < size*size;i++)
		{
			StringNum = i / size;
			rank = StringNum % ProcNum;
			if(rank!=0)
				MPI_Send(&pMatrix[i],1,MPI_FLOAT,rank,0,MPI_COMM_WORLD);
		}		
		for(int j=0;j < RowNum;j++)
			for(int k=0;k < size;k++)
				pProcRows[j*size+k]=pMatrix[pProcNum[j]*size+k];
	}
	else
		for(int i=0;i < RowNum*size;i++)
			MPI_Recv(&pProcRows[i],1,MPI_FLOAT,0,0,MPI_COMM_WORLD,&status);
	
	//рассылка 0-м процессом остальным процессам значений вектора pVector
	if(ProcRank==0)
	{
		int rank;
		for(int i=0;i < size;i++)
		{
			rank = i % ProcNum;
			if(rank!=0)
				MPI_Send(&pVector[i],1,MPI_FLOAT,rank,1,MPI_COMM_WORLD);
		}
		for(int j=0;j < RowNum;j++)
			pProcVector[j]=pVector[pProcNum[j]];
	}
	else
		for(int i=0;i < RowNum;i++)
			MPI_Recv(&pProcVector[i],1,MPI_FLOAT,0,1,MPI_COMM_WORLD,&status);
}

void GaussForward() 
{	
	float maxProc;
	int StringNumMaxProc;
	float max;
	int StringNumMax;	//номер строки матрицы pMatrix, содержащей максимальный элемент на данном процессе
	float* swap;
	float swap_vect;
	int Proc_StringNumMax;
	int Proc_i;
	int rank;
	float coef;
	MPI_Status status;
	for(int i=0;i < size-1;i++)
	{
		//определение максимального элемента в столбце i (ниже главной диагонали) в строках данного процесса
		StringNumMax = size-1;		
		if(ProcRank == 0)
			max = fabs(pMatrix[i+(size-1)*size]);
		MPI_Bcast(&max,1,MPI_FLOAT,0,MPI_COMM_WORLD);
		for(int j=RowNum-1;j > =0;j--)
		{
			if(pProcNum[j] < i)
				break;
			if(max < fabs(pProcRows[i+j*size]))
			{
				max = fabs(pProcRows[i+j*size]);
				StringNumMax = pProcNum[j];
			}
		}
		//отправка максимального значения и номера строки 0-му процессу
		if(ProcRank != 0)
		{
			MPI_Send(&max,1,MPI_FLOAT,0,2,MPI_COMM_WORLD);
			MPI_Send(&StringNumMax,1,MPI_INT,0,3,MPI_COMM_WORLD);
		}
		//получение максимального значения и номера строки 0-ым процессом
		if(ProcRank == 0)
		{			
			for(int rank=1;rank < ProcNum;rank++)
			{
				MPI_Recv(&maxProc,1,MPI_FLOAT,rank,2,MPI_COMM_WORLD,&status);
				MPI_Recv(&StringNumMaxProc,1,MPI_INT,rank,3,MPI_COMM_WORLD,&status);
				if(max < maxProc)
					{
						max = maxProc;
						StringNumMax = StringNumMaxProc;
					}
			}			
		}
		MPI_Bcast(&StringNumMax,1,MPI_INT,0,MPI_COMM_WORLD);
		
		//определение номеров процессов, которым принадлежат строки, которые необходимо поменять местами
		Proc_StringNumMax = StringNumMax % ProcNum;
		Proc_i = i % ProcNum;
		if((StringNumMax != i) && (ProcRank == Proc_StringNumMax || ProcRank == Proc_i))
		{	
			//меняем местами строки StringNumMax и i
			if(Proc_StringNumMax == Proc_i)
			{
				swap = new float[size];
				int k_i = -1;
				int k_StringNumMax = -1;
				for(int j=0;j < RowNum;j++)
				{
					if(pProcNum[j] == i)
						k_i = j;
					if(pProcNum[j] == StringNumMax)
						k_StringNumMax = j;
				}
				for(int j=0;j < size;j++)
				{
					swap[j] = pProcRows[j+k_i*size];
					pProcRows[j+k_i*size] = pProcRows[j+k_StringNumMax*size];
					pProcRows[j+k_StringNumMax*size] = swap[j];					
				}
				swap_vect = pProcVector[k_i];
				pProcVector[k_i] = pProcVector[k_StringNumMax];
				pProcVector[k_StringNumMax] = swap_vect;
			}
			else
			{
				if(ProcRank == Proc_StringNumMax)
				{
					//передача своей строки процессу Proc_i
					int k = -1;
					for(int j=0;j < RowNum;j++)
						if(pProcNum[j] == StringNumMax)
						{
							k = j;
							break;
						}
					for(int j=0;j < size;j++)
							MPI_Send(&pProcRows[j+k*size],1,MPI_FLOAT,Proc_i,4,MPI_COMM_WORLD);
					MPI_Send(&pProcVector[k],1,MPI_FLOAT,Proc_i,7,MPI_COMM_WORLD);
				}
				if(ProcRank == Proc_i)
				{
					//запись строки процесса Proc_StringNumMax в массив swap
					swap = new float[size];
					for(int j=0;j < size;j++)
						MPI_Recv(&swap[j],1,MPI_FLOAT,Proc_StringNumMax,4,MPI_COMM_WORLD,&status);
					MPI_Recv(&swap_vect,1,MPI_FLOAT,Proc_StringNumMax,7,MPI_COMM_WORLD,&status);
					//передача своей строки процессу Proc_StringNumMax
					int k = -1;
					for(int j=0;j < RowNum;j++)
						if(pProcNum[j] == i)
						{
							k = j;
							break;
						}
					for(int j=0;j < size;j++)
						MPI_Send(&pProcRows[j+k*size],1,MPI_FLOAT,Proc_StringNumMax,5,MPI_COMM_WORLD);
					MPI_Send(&pProcVector[k],1,MPI_FLOAT,Proc_StringNumMax,8,MPI_COMM_WORLD);
				}
				if(ProcRank == Proc_StringNumMax)
				{
					//запись строки процесса Proc_i
					int k = -1;
					for(int j=0;j < RowNum;j++)
						if(pProcNum[j] == StringNumMax)
						{
							k = j;
							break;
						}
					for(int j=0;j < size;j++)
						MPI_Recv(&pProcRows[j+k*size],1,MPI_FLOAT,Proc_i,5,MPI_COMM_WORLD,&status);
					MPI_Recv(&pProcVector[k],1,MPI_FLOAT,Proc_i,8,MPI_COMM_WORLD,&status);
				}
				if(ProcRank == Proc_i)
				{
					int k = -1;
					for(int j=0;j < RowNum;j++)
						if(pProcNum[j] == i)
						{
							k = j;
							break;
						}
					for(int j=0;j < size;j++)
						pProcRows[j+k*size] = swap[j];
					pProcVector[k] = swap_vect;
				}
			}
		}	
		//
		//MPI_Barrier(MPI_COMM_WORLD);
		//рассылка строки i всем процессам, которым принадлежат строки c номером, большим i
		if(ProcRank == Proc_i)
		{
			int k = -1;
			for(int j=0;j < RowNum;j++)
				if(pProcNum[j] == i)
				{
					k = j;
					break;
				}
			for(int string=i+1;string < size;string++)
			{	
				rank = string % ProcNum;
				if(rank != Proc_i)
				{
					for(int j=0;j < size;j++)
						MPI_Send(&pProcRows[j+k*size],1,MPI_FLOAT,rank,6,MPI_COMM_WORLD);
					MPI_Send(&pProcVector[k],1,MPI_FLOAT,rank,9,MPI_COMM_WORLD);
				}
				else
					for(int j=0;j < RowNum;j++)
						if(pProcNum[j] == string)
							if(pProcRows[i+j*size] != 0)
							{
								coef = pProcRows[i+j*size]/pProcRows[i+k*size];
								pProcRows[i+j*size] = 0;
								pProcVector[j] = pProcVector[j] - coef*pProcVector[k];
								for(int m=i+1;m < size;m++)
									pProcRows[m+j*size] = pProcRows[m+j*size] - coef*pProcRows[m+k*size];						
							}
			}
		}
		if(ProcRank != Proc_i)
		{
			for(int j=0;j < RowNum;j++)
				if(pProcNum[j]  >  i)
				{
					swap = new float[size];
					for(int k=0;k < size;k++)
						MPI_Recv(&swap[k],1,MPI_FLOAT,Proc_i,6,MPI_COMM_WORLD,&status);
					MPI_Recv(&swap_vect,1,MPI_FLOAT,Proc_i,9,MPI_COMM_WORLD,&status);
					if(pProcRows[i+j*size] != 0)
					{
						coef = pProcRows[i+j*size]/swap[i];
						pProcRows[i+j*size] = 0;
						pProcVector[j] = pProcVector[j] - coef*swap_vect;
						for(int k=i+1;k < size;k++)
							pProcRows[k+j*size] = pProcRows[k+j*size] - coef*swap[k];
					}
				}
		}		
	}	
}

void GaussBack()
{
	int Proc_i;
	int rank;
	float res;
	MPI_Status status;
	for(int i=size-1;i > 0;i--)
	{
		Proc_i = i % ProcNum;
		if(ProcRank == Proc_i)
		{
			int k;
			for(int j=0;j < RowNum;j++)
				if(pProcNum[j] == i)
				{
					k = j;
					break;
				}
			pProcVector[k] = pProcVector[k]/pProcRows[i+k*size];
			pProcRows[i+k*size] = 1;
			for(int j=i-1;j > =0;j--)
			{
				rank = j % ProcNum;
				if(rank != Proc_i)
					MPI_Send(&pProcVector[k],1,MPI_FLOAT,rank,10,MPI_COMM_WORLD);
				else
				{
					int k_j;
					for(int m=0;m < RowNum;m++)
						if(pProcNum[m] == j)
						{
							k_j = m;
							break;
						}
					pProcVector[k_j] = pProcVector[k_j] - pProcRows[i+k_j*size]*pProcVector[k];
					pProcRows[i+k_j*size] = 0;
				}
			}			 
		}
		else
		{
			for(int j=RowNum-1;j > =0;j--)
				if(pProcNum[j] < i)
				{
					MPI_Recv(&res,1,MPI_FLOAT,Proc_i,10,MPI_COMM_WORLD,&status);
					pProcVector[j] = pProcVector[j] - pProcRows[i+j*size]*res;
					pProcRows[i+j*size] = 0;
				}			
		}		
		//MPI_Barrier(MPI_COMM_WORLD);		
	}
	if(ProcRank == 0)
	{
		pProcVector[0] = pProcVector[0]/pProcRows[0];
		pProcRows[0] = 1;
	}	
}

void Result()
{
	int rank;
	int k;
	for(int j=0;j < size;j++)
	{
		rank = j % ProcNum;
		if(ProcRank == rank)
			{				
				for(int m=0;m < RowNum;m++)
					if(pProcNum[m] == j)
					{
						k = m;
						break;
					}
				printf("\nResult[%i] = %f",j,pProcVector[k]);
			}
		MPI_Barrier(MPI_COMM_WORLD);
	}
}

void main(int argc, char* argv[])
{
	setvbuf(stdout, 0, _IONBF, 0);

	MPI_Init(&argc, &argv);
	MPI_Comm_rank (MPI_COMM_WORLD, &ProcRank);
	MPI_Comm_size (MPI_COMM_WORLD, &ProcNum);

	if (ProcRank == 0)
		printf("Gauss method is starting!");

	ProcessInitialization();
	DataDistribution();
	MPI_Barrier(MPI_COMM_WORLD);
	start = MPI_Wtime();
	GaussForward();
	GaussBack();
	//MPI_Barrier(MPI_COMM_WORLD);
	finish = MPI_Wtime();
	duration = finish - start;
	Result();
	if (ProcRank == 0)
	{
		printf("\n\nGauss method has finished!");
		printf("\n\nTime: %f",duration);
		getchar();
		getchar();
	}

	MPI_Finalize();
}