Logo GenDocs.ru

Поиск по сайту:  

Загрузка...

Параллельное программировани MPI. Решение СЛАУ методом сопряженных градиентов - файл Цель работы.doc


Параллельное программировани MPI. Решение СЛАУ методом сопряженных градиентов
скачать (27.9 kb.)

Доступные файлы (4):

congrad.c
mdatcg.inp
vdatcg.inp
Цель работы.doc77kb.21.01.2010 00:38скачать

содержание
Загрузка...

Цель работы.doc

Реклама MarketGid:
Загрузка...
Цель работы: получить базовые навыки программирования для параллельных вычислительных систем. Разработать программу решения системы линейных уравнений с помощью средств MPI методом сопряженных градиентов.

MPI - это стандарт на программный инструментарий для обеспечения связи между ветвями параллельного приложения.

MPI расшифровывается как "Message passing interface" ("Взаимодействие через передачу сообщений"). Несколько путает дело тот факт, что этот термин уже применяется по отношению к аппаратной архитектуре ЭВМ. Программный инструментарий MPI реализован в том числе и для ЭВМ с такой архитектурой.

MPI предоставляет программисту единый механизм взаимодействия ветвей внутри параллельного приложения независимо от машинной архитектуры (однопроцессорные / многопроцессорные с общей/раздельной памятью), взаимного расположения ветвей (на одном процессоре / на разных) и API операционной системы.

( API = "applications programmers interface" = "интерфейс разработчика приложений" )

Программа, использующая MPI, легче отлаживается (сужается простор для совершения стереотипных ошибок параллельного программирования) и быстрее переносится на другие платформы (в идеале, простой перекомпиляцией).

В настоящее время разными коллективами разработчиков написано несколько программных пакетов, удовлетворяющих спецификации MPI, в частности: MPICH, LAM, HPVM и так далее. Они выступают базовыми при переносе MPI на новые архитектуры ЭВМ. Здесь в пособии рассматриваются разновидности MPICH. Это сделано по двум причинам:

MPICH написан авторами спецификации, и наиболее распространен.

Таким образом, далее по тексту термин MPI используется не только для обозначения изложенных в спецификации сведений, но и, в определенной мере, для описания характеристик конкретной базовой реализации.

Минимально в состав MPI входят: библиотека программирования (заголовочные и библиотечные файлы для языков Си, Си++ и Фортран) и загрузчик приложений.

^ Метод сопряженных градиентов

Может быть применен для решения симметричной положительно определенной системы линейных уравнений

Матрица А является симметричной, если она совпадает со своей транспонированной

матрицей, т.е. А=АТ. Матрица А называется положительно определенной, если для любого вектора x справедливо: xT A x >0 .

Известно, что после выполнения n

Итераций алгоритма ( n есть порядок решаемой системы линейных уравнений ), очередное приближение xn совпадает с точным решением.

Последовательный алгоритм

Если матрица A симметричная и положительно определенная, то функция (8.6)



имеет единственный минимум, который достигается в точке x*, совпадающей с решением системы линейных уравнений (8.2). Метод сопряженных градиентов является одним из широкого класса итерационных алгоритмов, которые позволяют найти решение (8.2) путем минимизации функции q(x).

Итерация метода сопряженных градиентов состоит в вычислении очередного приближения к точному решению в соответствии с правилом: (8.7)



Тем самым, новое значение приближения xk вычисляется с учетом приближения, построенного на предыдущем шаге xk-1, скалярного шага sk и вектора направления dk.

Перед выполнением первой итерации векторы x0 и d0 полагаются равными нулю, а для вектора g0 устанавливается значение, равное - b. Далее каждая итерация для вычисления очередного значения приближения xk включает выполнение четырех шагов:

Шаг 1: Вычисление градиента: (8.8)



Шаг 2: Вычисление вектора направления:

где (gT,g) есть скалярное произведение векторов.

Шаг 3: Вычисление величины смещения по выбранному направлению: (8.10)



Шаг 4: Вычисление нового приближения: (8.11)



Как можно заметить, данные выражения включают две операции умножения матрицы на вектор, четыре операции скалярного произведения и пять операций над векторами. Как результат, общее количество операций, выполняемых на одной итерации, составляет

t1=2n2+13n.

Как уже отмечалось ранее, для нахождения точного решения системы линейных уравнений с положительно определенной симметричной матрицей необходимо выполнить n итераций. Таким образом, для нахождения решения системы необходимо выполнить (8.12)

и, тем самым, сложность алгоритма имеет порядок O(n3).

Поясним выполнение метода сопряженных градиентов на примере решения системы линейных уравнений вида:

3x0 -x1 = 3,

-x0 +3x1 = 7.

Эта система уравнений второго порядка обладает симметричной положительно определенной матрицей, для нахождения точного решения этой системы достаточно провести всего две итерации метода.

На первой итерации было получено значение градиента g1=(-3,-7), значение вектора направления d1=(3, 7), значение величины смещения s1=0,439. Соответственно, очередное приближение к точному решению системы x1=(1,318, 3,076).

На второй итерации было получено значение градиента g2=(-2,121, 0,909), значение вектора направления d2=(2,397, -0,266), а величина смещения – s2=0,284. Очередное приближение совпадает с точным решением системы x2=(2, 3).

^ Организация параллельных вычислений

При разработке параллельного варианта метода сопряженных градиентов для решения систем линейных уравнений в первую очередь следует учесть, что выполнение итераций метода осуществляется последовательно и, тем самым, наиболее целесообразный подход состоит в распараллеливании вычислений, реализуемых в ходе выполнения итераций.

Анализ соотношений (8.8) – (8.11) показывает, что основные вычисления, выполняемые в соответствии с методом, состоят в умножении матрицы A на векторы x и d, и, как результат, при организации параллельных вычислений может быть полностью использован материал, изложенный в лекции 6.

Дополнительные вычисления в (8.8) – (8.11), имеющие меньший порядок сложности, представляют собой различные операции обработки векторов (скалярное произведение, сложение и вычитание, умножение на скаляр). Организация таких вычислений, конечно же, должна быть согласована с выбранным параллельным способом выполнения операции умножения матрицы на вектор. Общие же рекомендации могут состоять в следующем: при малом размере векторов можно применить дублирование векторов между процессорами, при большом порядке решаемой системы линейных уравнений более целесообразно осуществлять блочное разделение векторов.

Программа:
#include<math.h>

#include<stdio.h>

#include<stdlib.h>

#include<mpi.h>
#define EPSILON 1.0E-20

#define MAX_ITERATIONS 10000
/******************************************************************************/

void

GetPreconditionMatrix(double **Bloc_Precond_Matrix, int NoofRows_Bloc,

int NoofCols)

{

/*... Preconditional Martix is identity matrix .......*/

int Bloc_MatrixSize;

int irow, icol, index;

double *Precond_Matrix;
Bloc_MatrixSize = NoofRows_Bloc*NoofCols;
Precond_Matrix = (double *) malloc(Bloc_MatrixSize * sizeof(double));
index = 0;

for(irow=0; irow<NoofRows_Bloc; irow++){

for(icol=0; icol<NoofCols; icol++){

Precond_Matrix[index++] = 1.0;

}

}

*Bloc_Precond_Matrix = Precond_Matrix;

}

/******************************************************************************/

double

ComputeVectorDotProduct(double *Vector1, double *Vector2, int VectorSize)

{

int index;

double Product;
Product = 0.0;

for(index=0; index<VectorSize; index++)

Product += Vector1[index]*Vector2[index];
return(Product);

}

/******************************************************************************/

void

CalculateResidueVector(double *Bloc_Residue_Vector, double *Bloc_Matrix_A, double *Input_B, double *Vector_X, int NoofRows_Bloc, int VectorSize, int MyRank)

{

/*... Computes residue = AX - b .......*/

int irow, index, GlobalVectorIndex;

double value;
GlobalVectorIndex = MyRank * NoofRows_Bloc;

for(irow=0; irow<NoofRows_Bloc; irow++){

index = irow * VectorSize;

value = ComputeVectorDotProduct(&Bloc_Matrix_A[index], Vector_X,

VectorSize);

Bloc_Residue_Vector[irow] = value - Input_B[GlobalVectorIndex++];

}

}

/******************************************************************************/

void

SolvePrecondMatrix(double *Bloc_Precond_Matrix, double *HVector, double *Bloc_Residue_Vector, int Bloc_VectorSize)

{

/*...HVector = Bloc_Precond_Matrix inverse * Bloc_Residue_Vector.......*/

int index;
for(index=0; index<Bloc_VectorSize; index++){

HVector[index] = Bloc_Residue_Vector[index]/1.0;

}

}

/******************************************************************************/

main(int argc, char *argv[])

{

int NumProcs, MyRank, Root=0;

int NoofRows, NoofCols, VectorSize;

int n_size, NoofRows_Bloc, Bloc_MatrixSize, Bloc_VectorSize;

int Iteration = 0, irow, icol, index, CorrectResult;

double **Matrix_A, *Input_A, *Input_B, *Vector_X, *Bloc_Vector_X;

double *Bloc_Matrix_A, *Bloc_Precond_Matrix, *Buffer;

double *Bloc_Residue_Vector ,*Bloc_HVector, *Bloc_Gradient_Vector;

double *Direction_Vector, *Bloc_Direction_Vector;

double Delta0, Delta1, Bloc_Delta0, Bloc_Delta1;

double Tau, val, temp, Beta;
double *AMatDirect_local, *XVector_local;

double *ResidueVector_local, *DirectionVector_local;

double StartTime, EndTime;

MPI_Status status;

FILE *fp;
/*...Initialising MPI .......*/

MPI_Init(&argc,&argv);

MPI_Comm_size(MPI_COMM_WORLD,&NumProcs);

MPI_Comm_rank(MPI_COMM_WORLD,&MyRank);


/* .......Read the Input file ......*/

if(MyRank == Root) {
if ((fp = fopen ("./data/mdatcg.inp", "r")) == NULL) {

printf("Can't open input matrix file");

exit(-1);

}

fscanf(fp, "%d %d", &NoofRows,&NoofCols);

n_size=NoofRows;
/* ...Allocate memory and read data .....*/

Matrix_A = (double **) malloc(n_size*sizeof(double *));
for(irow = 0; irow < n_size; irow++){

Matrix_A[irow] = (double *) malloc(n_size * sizeof(double));

for(icol = 0; icol < n_size; icol++)

fscanf(fp, "%lf", &Matrix_A[irow][icol]);

}

fclose(fp);
if ((fp = fopen ("./data/vdatcg.inp", "r")) == NULL){

printf("Can't open input vector file");

exit(-1);

}
fscanf(fp, "%d", &VectorSize);

n_size=VectorSize;

Input_B = (double *)malloc(n_size*sizeof(double));

for (irow = 0; irow<n_size; irow++)

fscanf(fp, "%lf",&Input_B[irow]);

fclose(fp);
/* ...Convert Matrix_A into 1-D array Input_A ......*/

Input_A = (double *)malloc(n_size*n_size*sizeof(double));

index = 0;

for(irow=0; irow<n_size; irow++)

for(icol=0; icol<n_size; icol++)

Input_A[index++] = Matrix_A[irow][icol];
}
MPI_Barrier(MPI_COMM_WORLD);

StartTime = MPI_Wtime();
/*...Broadcast Matrix and Vector size and perform input validation tests...*/

MPI_Bcast(&NoofRows, 1, MPI_INT, Root, MPI_COMM_WORLD);

MPI_Bcast(&NoofCols, 1, MPI_INT, Root, MPI_COMM_WORLD);

MPI_Bcast(&VectorSize, 1, MPI_INT, Root, MPI_COMM_WORLD);


if(NoofRows != NoofCols){

MPI_Finalize();

if(MyRank == Root)

printf("Error : Coefficient Matrix Should be square matrix");

exit(-1);

}
if(NoofRows != VectorSize){

MPI_Finalize();

if(MyRank == Root)

printf("Error : Matrix Size should be equal to VectorSize");

exit(-1);

}
if(NoofRows % NumProcs != 0){

MPI_Finalize();

if(MyRank == Root)

printf("Error : Matrix cannot be evenly striped among processes");

exit(-1);

}
/*...Allocate memory for Input_B and BroadCast Input_B.......*/

if(MyRank != Root)

Input_B = (double *) malloc(VectorSize*sizeof(double));

MPI_Bcast(Input_B, VectorSize, MPI_DOUBLE, Root, MPI_COMM_WORLD);
/*...Allocate memory for Block Matrix A and Scatter Input_A .......*/

NoofRows_Bloc = NoofRows / NumProcs;

Bloc_VectorSize = NoofRows_Bloc;

Bloc_MatrixSize = NoofRows_Bloc * NoofCols;

Bloc_Matrix_A = (double *) malloc (Bloc_MatrixSize*sizeof(double));

MPI_Scatter(Input_A, Bloc_MatrixSize, MPI_DOUBLE,

Bloc_Matrix_A, Bloc_MatrixSize, MPI_DOUBLE, Root, MPI_COMM_WORLD);


/*... Allocates memory for solution vector and intialise it to zero.......*/

Vector_X = (double *) malloc(VectorSize*sizeof(double));

for(index=0; index<VectorSize; index++)

Vector_X[index] = 0.0;
/*...Calculate RESIDUE = AX - b .......*/

Bloc_Residue_Vector = (double *) malloc(NoofRows_Bloc*sizeof(double));

CalculateResidueVector(Bloc_Residue_Vector,Bloc_Matrix_A, Input_B, Vector_X, NoofRows_Bloc, VectorSize, MyRank);
/*... Precondtion Matrix is identity matrix ......*/

GetPreconditionMatrix( &Bloc_Precond_Matrix, NoofRows_Bloc, NoofCols);
/*...Bloc_HVector = Bloc_Precond_Matrix inverse * Bloc_Residue_Vector......*/

Bloc_HVector = (double *) malloc(Bloc_VectorSize*sizeof(double));

SolvePrecondMatrix(Bloc_Precond_Matrix, Bloc_HVector, Bloc_Residue_Vector, Bloc_VectorSize);


/*...Initailise Bloc Direction Vector = -(Bloc_HVector).......*/

Bloc_Direction_Vector = (double *) malloc(Bloc_VectorSize*sizeof(double));

for(index=0; index<Bloc_VectorSize; index++)

Bloc_Direction_Vector[index] = 0 - Bloc_HVector[index];
/*...Calculate Delta0 and check for convergence .......*/

Bloc_Delta0 = ComputeVectorDotProduct(Bloc_Residue_Vector,Bloc_HVector, Bloc_VectorSize);

MPI_Allreduce(&Bloc_Delta0, &Delta0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
if(Delta0 < EPSILON){

MPI_Finalize();

exit(0);

}
/*...Allocate memory for Direction Vector.......*/

Direction_Vector = (double *) malloc(VectorSize*sizeof(double));
/*...Allocate temporary buffer to store Bloc_Matrix_A*Direction_Vector...*/

Buffer = (double *) malloc(Bloc_VectorSize*sizeof(double));
/*...Allocate memory for Bloc_Vector_X .......*/

Bloc_Vector_X = (double *)malloc(Bloc_VectorSize*sizeof(double));
Iteration = 0;

do{
Iteration++;

/*

if(MyRank == Root)

printf("Iteration : %d\n",Iteration);

*/
/*...Gather Direction Vector on all processes.......*/

MPI_Allgather(Bloc_Direction_Vector, Bloc_VectorSize, MPI_DOUBLE,

Direction_Vector, Bloc_VectorSize, MPI_DOUBLE, MPI_COMM_WORLD);
/*...Compute Tau = Delta0 / (DirVector Transpose*Matrix_A*DirVector)...*/

for(index=0; index<NoofRows_Bloc; index++){

Buffer[index] = ComputeVectorDotProduct(&Bloc_Matrix_A[index*NoofCols],

Direction_Vector, VectorSize);

}

temp = ComputeVectorDotProduct(Bloc_Direction_Vector, Buffer, Bloc_VectorSize);
MPI_Allreduce(&temp, &val, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

Tau = Delta0 / val;
/*Compute new vector Xnew = Xold + Tau*Direction........................*/

/*Compute BlocResidueVec = BlocResidueVect + Tau*Bloc_MatA*DirVector...*/

for(index = 0; index<Bloc_VectorSize; index++){

Bloc_Vector_X[index] = Vector_X[MyRank*Bloc_VectorSize + index] + Tau*Bloc_Direction_Vector[index];

Bloc_Residue_Vector[index] = Bloc_Residue_Vector[index] + Tau*Buffer[index];

}

/*...Gather New Vector X at all processes......*/

MPI_Allgather(Bloc_Vector_X, Bloc_VectorSize, MPI_DOUBLE, Vector_X,

Bloc_VectorSize, MPI_DOUBLE, MPI_COMM_WORLD);
SolvePrecondMatrix(Bloc_Precond_Matrix, Bloc_HVector, Bloc_Residue_Vector, Bloc_VectorSize);

Bloc_Delta1 = ComputeVectorDotProduct(Bloc_Residue_Vector, Bloc_HVector, Bloc_VectorSize);

MPI_Allreduce(&Bloc_Delta1, &Delta1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
if(Delta1 < EPSILON)

break;
Beta = Delta1 / Delta0;

Delta0 = Delta1;

for(index=0; index<Bloc_VectorSize; index++){

Bloc_Direction_Vector[index] = -Bloc_HVector[index] + Beta*Bloc_Direction_Vector[index];

}

}while(Delta0 > EPSILON && Iteration < MAX_ITERATIONS);
MPI_Barrier(MPI_COMM_WORLD);

EndTime = MPI_Wtime();
if (MyRank == 0) {
printf ("\n");

printf(" ------------------------------------------- \n");

printf("Results on processor %d: \n", MyRank);

printf ("\n");
printf("Matrix Input_A \n");

printf ("\n");

for (irow = 0; irow < n_size; irow++) {

for (icol = 0; icol < n_size; icol++)

printf("%.3lf ", Matrix_A[irow][icol]);

printf("\n");

}

printf ("\n");

printf("Matrix Input_B \n");

printf("\n");

for (irow = 0; irow < n_size; irow++) {

printf("%.3lf\n", Input_B[irow]);

}

printf ("\n");

printf("Solution vector \n");

printf("Number of iterations = %d\n",Iteration);

printf ("\n");

for(irow = 0; irow < n_size; irow++)

printf("%.12lf\n",Vector_X[irow]);

printf(" --------------------------------------------------- \n");

}

MPI_Finalize();

}


Скачать файл (27.9 kb.)

Поиск по сайту:  

© gendocs.ru
При копировании укажите ссылку.
обратиться к администрации