Решенние систем линейных уравнений методом Якоби

Отчет по лабораторной работе по курсу "Многопроцессорные вычислительные системы и параллельное программирование"

Тема "Решение системы СЛУ итеративным методом Якоби"

Выполнили студенты группы 84-10 Штыркова О.А., Сериков Е.В.

Постановка задачи

Дана система линейных уравнений:
, где

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

где в принятых обозначениях D означает матрицу, у которой на главной диагонали стоят соответствующие элементы матрицы A, а все остальные нули; тогда как матрицы U и L содержат верхнюю и нижнюю треугольные части A, на главной диагонали которых нули, E — единичная матрица. Тогда процедура нахождения решения имеет вид: где k счётчик итерации. Условие окончания итерационного процесса при достижении заданной точности имеет вид: .

Описание схемы параллельной реализации

Динамическая визуализация

Теоретический анализ эффективности

Результаты проведенных экспериментов

Матрица	Число итераций	Последовательные (в 1 процесс)	Параллельные в 2 процесса	Ускорение при 2 процессаx	Параллельные в 4 процесса	Ускорение при 4 процессаx
50x50	366	0,37	0,19	0,18	0,17	0,2
100x100	744	3,1	1,5	1,6	1,4	1,7

Исходный код программы

#include "stdafx.h"
#include "math.h"
#include "windows.h"
#include "stdio.h"
#include "iostrea"
#include "./include/mpi.h"
#pragma comment(lib,"lib/i386/msmpi.lib")

int      N      = 0;
double   EPS    = 1; 
double*  XS     = NULL;
double*  X      = NULL;
double** A      = NULL;
double*  B      = NULL;
void Jacobi (int rank);

void initProcData( );
void calc( int rank );
int getRangeLen( int rank)
{
    int size = 0;
    MPI_Comm_size( MPI_COMM_WORLD, &size);
    return (N/size)+((rank<(N-(N/size)*size))?1:0);
}
int getStartIndex(int rank)
{
    int index = 0;
    for(int i = 0; i norm)
                    norm = fabs(X[h] - XS[h]);
                XS[h] = X[h];
            }

        }while(norm>EPS);
        int flag = 0;
        
        for(int i = 1; i< size;i++)
        {
            MPI_Send(XS, N, MPI_DOUBLE, i, 99, MPI_COMM_WORLD);
            MPI_Send(&flag, 1, MPI_INT, i, 98, MPI_COMM_WORLD);
        }

        for(int i = 0; i< N;i++)
        {
            std::cout<< XS[i] <<",";
        }

    } 
    else                /* code for process one */ 
    {
        int flag = 0;
        do {

            MPI_Send(X+getStartIndex(rank), getRangeLen(rank), MPI_DOUBLE, 0, 99, MPI_COMM_WORLD);
            MPI_Recv(XS, N, MPI_DOUBLE, 0, 99, MPI_COMM_WORLD, &status); 
            MPI_Recv(&flag, 1, MPI_INT, 0, 98, MPI_COMM_WORLD, &status); 
            calc( rank );

        }while(flag);
    } 


    MPI_Finalize();


}
void calc( int rank )
{

    Jacobi (rank);

}
 
void initProcData( )
{
    char buf[256] = {0};
    char tmp[256] = {0};
    FILE *file = fopen("input.txt","r");
    if(!file)
        return;
    fgets(buf,256,file);
    N=atoi(buf);
    ZeroMemory(buf,256);

    fgets(buf,256,file);
    EPS=atof(buf);
    ZeroMemory(buf,256);

    XS = new double[N];
    X = new double[N];
    A = new double*[N];
    for(int i = 0; i< N;i++)
    {
        A[i] = new double[N];
    }
    B = new double[N];

    for(int i = 0; i< N;i++)
    {
        fgets(buf,256,file);
        char*softbuf = buf;
        if(strlen(buf) == 0)
            return;
        int len = strlen(buf);
        for(int k = 0; k< len; k++)
        {
            if( softbuf[k] == ',')
                softbuf[k] = 0;
        }
        for(int j = 0; j< N; j++)
        {
            A[i][j] = atof(softbuf);
            softbuf += (strlen( softbuf )+1);
        }
        B[i] = atof(softbuf);
        ZeroMemory(buf,256);
    }

    fgets(buf,256,file);
    char*softbuf = buf;
    if(strlen(buf) == 0)
        return;
    int len = strlen(buf);
    for(int k = 0; k< len; k++)
    {
        if( softbuf[k] == ',')
            softbuf[k] = 0;
    }
    for(int j = 0; j< N; j++)
    {
        XS[j] = atof(softbuf);
        softbuf += (strlen( softbuf )+1);
    }
    ZeroMemory(buf,256);

    fclose(file);
   return;
}

void Jacobi (int rank)
{

        for (int i = getStartIndex(rank); i< getStartIndex(rank)+getRangeLen(rank); i++) {
            X[i] =- B[i];
            for (int g = 0; g < N; g++) {
                if (i != g)
                    X[i] += A[i][g] * XS[g];
            }
            X[i] /= -A[i][i];
        }
 }

Отчет