#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include <mpi.h>
double *matr, *vect, *res, *_matr, *_vect, *_res, m=0.0;
//_matr, _vect, _res - куски матрицы, вектора b и вектора x на определенном процессе
int size=8, p=1, status, str_count; //str_count-число строк исходной матрицы, хранящихся на конкретном процессе
//status - статус решения задачи
//stringcount - количество строк на процессе
double start, finish, delta; // переменные для засечения времени
int *IndexVedStr; //массив номеров ведущих строк на каждой итерации // i-ая итерация - j-ая строка ведущая
int *IndexVedIter; /* массив номеров итераций, на которых соответствующие строки системы, расположенные на процессе, выбраны ведущими */
int ProcCount, ProcRank, *ui1,//-массив размером с количество процессов. каждый элемент
//хранит индекс первого элемента из столбца свободных членов, который принадлежит процессу
*range; //массив с количествами пересылаемых элеметов из столбца свободных членов на каждый процесс
//------------------------------------------------------------------------------
void StartProcInitialization()
{
int show;
printf("\nInput matrix rank: ");
scanf("%i", &size);
//printf("\nInput srand argument: ");
//scanf("%i", &p);
printf("\nDo you want to see values?(1/0): ");
scanf("%i", &show);
//size=8;
srand(p);
//randomize();
matr = new double [size*size];
vect = new double [size];
res = new double [size];
for(int i=0; i<size*size; i++)
{
matr[i]=rand();
if(show==1) printf("\nmatr[%i]=%f", i, matr[i]);
};
for(int i=0; i<size; i++)
{
vect[i]=rand();
if(show==1) printf("\nb[%i]=%f", i, vect[i]);
res[i]=0;
};
}
//------------------------------------------------------------------------------
void ProcInitialization()
{
int nstrings; //Число строк, ещё не распределённых по процессам
if (ProcRank == 0)
{
printf("\nStarted!");
StartProcInitialization();
}
MPI_Bcast(&size, 1, MPI_INT, 0, MPI_COMM_WORLD);
//Определение размера части данных, расположенных на конкретном процессе
nstrings = size;
for (int i=0; i<ProcRank; i++)
nstrings = nstrings-nstrings/(ProcCount-i);
str_count = nstrings/(ProcCount-ProcRank);
_matr = new double [str_count*size];//выделяем память под строки матрицы
_vect = new double [str_count]; //память под элементы столбца свободных членов
_res = new double [str_count]; //память под элементы вектора результата
IndexVedStr = new int [size]; //массив индексов ведущих строк системы на каждой итерации
IndexVedIter = new int [str_count]; /* итерация, на которой соответствующая строка системы, расположенная на процессе, выбрана ведущей */
ui1 = new int [ProcCount];
range = new int [ProcCount];
for (int i=0; i<str_count; i++)
IndexVedIter[i] = -1;
}
//------------------------------------------------------------------------------
// Распределение исходных данных между процессами
void Distribution()
{
int *firstsendindex; //Индекс первого элемента матрицы, передаваемого процессу
int *sendrange; //Число элементов матрицы, передаваемых процессу
int CountProcessstrings;
int nstrings;
firstsendindex = new int [ProcCount];
sendrange = new int [ProcCount];
nstrings = size;
for (int i=0; i<ProcCount; i++) //Определяем, сколько элементов матрицы будет передано каждому процессу
{
int Previoussize = (i==0) ? 0 : sendrange[i-1];
int PreviousIndex = (i==0) ? 0 : firstsendindex[i-1];
int Portion = (i==0) ? 0 : CountProcessstrings; //число строк, отданных предыдущему процессу
nstrings -= Portion;
CountProcessstrings = nstrings/(ProcCount-i);
sendrange[i] = CountProcessstrings*size;
firstsendindex[i] = PreviousIndex+Previoussize;
};
MPI_Scatterv(matr, sendrange, firstsendindex, MPI_DOUBLE, _matr, sendrange[ProcRank], MPI_DOUBLE, 0, MPI_COMM_WORLD);
//Разослали
nstrings = size;
for (int i=0; i<ProcCount; i++)
{
int Previoussize = (i==0) ? 0 : range[i-1];
int PreviousIndex = (i==0) ? 0 : ui1[i-1];
int Portion = (i==0) ? 0 : range[i-1];
nstrings -= Portion;
range[i] = nstrings/(ProcCount-i);
ui1[i] = PreviousIndex+Previoussize;
};
//Рассылка вектора
/*
-что рассылать
-массив с количествами пересылаемых элеметов на каждый процесс
-массив с номерами элеметов, с которых надо начинать рассылку из vect
-куда принять
-сколько элементов на каждый процесс
*/
MPI_Scatterv(vect, range, ui1, MPI_DOUBLE, _vect, range[ProcRank], MPI_DOUBLE, 0, MPI_COMM_WORLD);
status = 1;
MPI_Bcast(&status, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&m, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
delete [] sendrange;
delete [] firstsendindex;
}
//------------------------------------------------------------------------------
void StolbElimination(int IterationNumber, double *globmaxSystemstring)
{
double delen;
for (int i=0; i<str_count; i++)//для каждой строки в процессе
{
if (IndexVedIter[i] == -1)
{
delen = _matr[i*size+IterationNumber] / globmaxSystemstring[IterationNumber];
for (int j=IterationNumber; j<size; j++)
{
_matr[i*size + j] -= globmaxSystemstring[j]*delen;
};
_vect[i] -= globmaxSystemstring[size]*delen;
};
};
}
//------------------------------------------------------------------------------
void GaussElimination()
{
int VedIndex; // индекс ведущей строки на конкретном процессе
struct { double MaxValue; int ProcRank; } Localmax, globmax; //максимальный элемент+номер процесса, у которого он
double* pglobmaxstring = new double [size+1]; //т.е. строка матрицы+значение вектора
for (int i=0; i<size; i++) //главный цикл решения
{
// Вычисление ведущей строки
double MaxValue = 0;
int tmp_index = -1;
for (int j=0; j<str_count; j++)
{
tmp_index = j;
if ((IndexVedIter[j] == -1)&&(MaxValue < fabs(_matr[i+size*j])))
{
MaxValue = fabs(_matr[i+size*j]);
VedIndex = j;
};
};
Localmax.MaxValue = MaxValue;
Localmax.ProcRank = ProcRank;
//каждый процесс рассылает свой локально максимальный элемент по всем столцам, все процесы принимают уже глобально максимальный элемент
MPI_Allreduce(&Localmax, &globmax, 1, MPI_DOUBLE_INT, MPI_MAXLOC, MPI_COMM_WORLD);
//Вычисление ведущей строки всей системы
if ( ProcRank == globmax.ProcRank )
{
if (globmax.MaxValue == 0)
{
status = 2;
MPI_Barrier(MPI_COMM_WORLD);
MPI_Send(&status, 1, MPI_INT, 0, MPI_ANY_TAG, MPI_COMM_WORLD);
IndexVedIter[tmp_index] = i;
IndexVedStr[i] = ui1[ProcRank] + VedIndex;
continue;
}
else
{
IndexVedIter[VedIndex] = i; // Номер итерации, на которой строка с локальным номером является ведущей для всей системы
IndexVedStr[i] = ui1[ProcRank] + VedIndex; //Вычисленный номер ведущей строки системы
};
};
MPI_Bcast(&IndexVedStr[i], 1, MPI_INT, globmax.ProcRank, MPI_COMM_WORLD);
if (ProcRank == globmax.ProcRank)
{
for (int j=0; j<size; j++)
pglobmaxstring[j] = _matr[VedIndex*size + j];
pglobmaxstring[size] = _vect[VedIndex];
};
//Рассылка ведущей строки всем процессам
MPI_Bcast(pglobmaxstring, size+1, MPI_DOUBLE, globmax.ProcRank, MPI_COMM_WORLD);
//Исключение неизвестных в столбце с номером i
StolbElimination(i, pglobmaxstring);
};
}
//------------------------------------------------------------------------------
void Frp(int stringIndex, // номер строки, которая была ведущей на определеной итерации
int &iterationprocrank, // процесс, на котором эта строка
int &IterationItervedindex) // локальный номер этой строки (в рамках одного процесса)
{
for (int i=0; i<ProcCount-1; i++) //Определяем ранг процесса, содержащего данную строку
{
if ((ui1[i]<=stringIndex)&&(stringIndex<ui1[i+1]))
iterationprocrank = i;
}
if (stringIndex >= ui1[ProcCount-1])
iterationprocrank = ProcCount-1;
IterationItervedindex = stringIndex - ui1[iterationprocrank];
}
//------------------------------------------------------------------------------
void GaussReverse()
{
int IterProcRank; // Ранг процесса, хранящего текущую ведущую строку
int IndexVed; // локальный на своем процессе номер текущей ведущей
double IterRes; // значение Xi, найденное на итерации
double val;
// Основной цикл
for (int i=size-1; i>=0; i--)
{
Frp(IndexVedStr[i], IterProcRank, IndexVed);
// Определили ранг процесса, содержащего текущую ведущую строку, и номер этой строки на процессе
// Вычисляем значение неизвестной
if (ProcRank == IterProcRank)
{
if (_matr[IndexVed*size+i] == 0)
{
if (_vect[IndexVed] == 0)
IterRes = m;
else
{
status = 0;
MPI_Barrier(MPI_COMM_WORLD);
MPI_Send(&status, 1, MPI_INT, 0, MPI_ANY_TAG, MPI_COMM_WORLD);
break;
};
}
else
IterRes = _vect[IndexVed]/_matr[IndexVed*size+i];
_res[IndexVed] = IterRes; //нашли значение переменной
};
MPI_Bcast(&IterRes, 1, MPI_DOUBLE, IterProcRank, MPI_COMM_WORLD);
for (int j=0; j<str_count; j++) //подстановка найденной переменной
if (IndexVedIter[j] < i)
{
val = _matr[size*j + i] * IterRes;
_vect[j] -= val;
};
};
}
//------------------------------------------------------------------------------
void Gauss()
{
GaussElimination();
GaussReverse();
//сбор данных, передача от всех одному (нулевому процессу)
MPI_Gatherv(_res, range[ProcRank], MPI_DOUBLE, res, range, ui1, MPI_DOUBLE, 0, 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, &ProcCount);
MPI_Barrier(MPI_COMM_WORLD);
ProcInitialization();
Distribution();
MPI_Barrier(MPI_COMM_WORLD);
start = MPI_Wtime();
Gauss();
MPI_Barrier(MPI_COMM_WORLD);
finish = MPI_Wtime();
delta = (finish-start);
if (ProcRank == 0)
{
if (status == 1) // решение найдено
{
printf("\nRoots: \n");
for (int i = 0; i<size; i++)
printf(" \nres[%i] = %f ", i, res[IndexVedStr[i]]);
};
if (status == 0) printf("\nNo roots.\n"); // решений нет
if (status == 2) // решений бесконечно много
{
printf("\nContinum. Solution: \n\n");
for (int i = 0; i<size; i++)
printf("\nres[%d] = %f ", i, res[IndexVedStr[i]]);
};
printf("\n\nTime= %f\n", delta);
};
if (ProcRank == 0)
{
printf("\nFinished!");
_getch();
};
MPI_Finalize();
delete [] matr;
delete [] vect;
delete [] res;
}
