<
#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include <mpi.h>
double *matrix, *vector, *result, *_matrix, *_vector, *_result, m=0.0;
//_matrix, _vector, _result - куски матрицы, вектора b и вектора x на определенном процессе
int size=8, p=1, status, str_count; //str_count-число строк исходной матрицы, хранящихся на конкретном процессе
//status - статус решения задачи
//stringcount - количество строк на процессе
double start, finish, deltime; // переменные для засечения времени
int *IndexVedStr; //массив номеров ведущих строк на каждой итерации // i-ая итерация - j-ая строка ведущая
int *IndexVedIter; /* массив номеров итераций, на которых соответствующие строки системы, расположенные на процессе, выбраны ведущими */
int ProcCount, ProcRank, *ui1,//-массив размером с количество процессов. каждый элемент
//хранит индекс первого элемента из столбца свободных членов, который принадлежит процессу
*range; //массив с количествами пересылаемых элеметов из столбца свободных членов на каждый процесс
//------------------------------------------------------------------------------
void StartProcInitialization()
{
int show;
printf("\nInput matrixix 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();
matrix = new double [size*size];
vector = new double [size];
result = new double [size];
for(int i=0; i<size*size; i++)
{
matrix[i]=rand();
if(show==1) printf("\nmatrix[%i]=%f", i, matrix[i]);
};
for(int i=0; i<size; i++)
{
vector[i]=rand();
if(show==1) printf("\nb[%i]=%f", i, vector[i]);
result[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);
_matrix = new double [str_count*size];//выделяем память под строки матрицы
_vector = new double [str_count]; //память под элементы столбца свободных членов
_result = 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(matrix, sendrange, firstsendindex, MPI_DOUBLE, _matrix, 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;
};
//Рассылка вектора
/*
-что рассылать
-массив с количествами пересылаемых элеметов на каждый процесс
-массив с номерами элеметов, с которых надо начинать рассылку из vector
-куда принять
-сколько элементов на каждый процесс
*/
MPI_Scatterv(vector, range, ui1, MPI_DOUBLE, _vector, 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 = _matrix[i*size+IterationNumber] / globmaxSystemstring[IterationNumber];
for (int j=IterationNumber; j<size; j++)
{
_matrix[i*size + j] -= globmaxSystemstring[j]*delen;
};
_vector[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(_matrix[i+size*j])))
{
MaxValue = fabs(_matrix[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] = _matrix[VedIndex*size + j];
pglobmaxstring[size] = _vector[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 IterResult; // значение Xi, найденное на итерации
double val;
// Основной цикл
for (int i=size-1; i>=0; i--)
{
Frp(IndexVedStr[i], IterProcRank, IndexVed);
// Определили ранг процесса, содержащего текущую ведущую строку, и номер этой строки на процессе
// Вычисляем значение неизвестной
if (ProcRank == IterProcRank)
{
if (_matrix[IndexVed*size+i] == 0)
{
if (_vector[IndexVed] == 0)
IterResult = m;
else
{
status = 0;
MPI_Barrier(MPI_COMM_WORLD);
MPI_Send(&status, 1, MPI_INT, 0, MPI_ANY_TAG, MPI_COMM_WORLD);
break;
};
}
else
IterResult = _vector[IndexVed]/_matrix[IndexVed*size+i];
_result[IndexVed] = IterResult; //нашли значение переменной
};
MPI_Bcast(&IterResult, 1, MPI_DOUBLE, IterProcRank, MPI_COMM_WORLD);
for (int j=0; j<str_count; j++) //подстановка найденной переменной
if (IndexVedIter[j] < i)
{
val = _matrix[size*j + i] * IterResult;
_vector[j] -= val;
};
};
}
//------------------------------------------------------------------------------
void Gauss()
{
GaussElimination();
GaussReverse();
//сбор данных, передача от всех одному (нулевому процессу)
MPI_Gatherv(_result, range[ProcRank], MPI_DOUBLE, result, 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);
if (ProcRank == 0) setvbuf(stdout, 0, _IONBF, 0);
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();
deltime = (finish-start);
if (ProcRank == 0)
{
if (status == 1) // решение найдено
{
printf("\nRoots: \n");
//for (int i = 0; i<size; i++)
//printf(" \nresult[%i] = %f ", i, result[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("\nresult[%d] = %f ", i, result[IndexVedStr[i]]);
};
printf("\n\nTime= %f\n", deltime);
};
if (ProcRank == 0)
{
printf("\nFinished!");
_getch();
};
MPI_Finalize();
delete [] matrix;
delete [] vector;
delete [] result;
}
