Рішення системи лінійних алгебраїчних рівнянь

#define count 4096 

bool WriteFile(float* a, int n, char* s);   // Функція зчитування даних з файлу в масив.

bool ReadFile(float* a, int n, char* s);   // Запис масиву в файл. 

int _tmain(int argc, _TCHAR* argv[])

{

      //  Виділяємо пам'ять на данні.

      float* matrix  = (float*)malloc(count * count * sizeof(float));

      float* vector = (float*)malloc(count * sizeof(float)); 

      //  Зчитуємо данні з файлу

      ReadFile(matrix, count * count, fileIn);

      ReadFile(vector, count, fileIn); 
 

      int n;  float d, y; 

     float* line;   //вказівник на рядок, що передається для вилучення невідомого 

//прямий хід

      for (int i = 0; i < count - 1; i++)

      {

            line = &matrix[i * count];

            y = vector[i];

//  Цикл проходження всіх рядків, що змінюються

            for (int j = i + 1; j < count; j++)

            {

                  n = j * count;

                  d = matrix[n + i];

//  вилучаємо i-ий невідомий j-го рядка

                  for (int k = i; k < count; k++)

                  {

                        matrix[n + k] = matrix[n + k] / d * line[i] - line[k];

                  }

                  vector[j] = vector[j] / d * line[i] - y;

            }

      }

//  виділення пам’яті на результат

      float* result = (float*)malloc(count * sizeof(float));

      float x; 

//  Зворотній хід

      for (int i = count - 1; i >= 0; i--)

      {

            x = vector[i];

            n = i * count;

//  Підрахунок i-го невідомого

            for (int j = i + 1; j < count; j++)

            {

                  x -= matrix[n + j];

            }

            result[i] = x /= matrix[n + i];

//  Запис знайденого невідомого до матриці, шляхом множення на коефіцієнти,

//  з якими він стоїть

            for (int j = 0; j < count; j++)

            {

                  matrix[j * count + i] *= x;

            }

      } 

//  Запис результату

      WriteFile(result, count, fileOut);

      printf("complited...\n"); 

//  Звільнення пам’яті

      free(matrix);

      free(vector);

      free(result);

      return 0;

} 
 

 

Паралельний метод

     Наведений метод майже нічим не відрізняється  від послідовного. По суті, ми просто розбили масив даних і кожен  процес окремо працює зі своїм блоком. Головна задача (і проблема) це розділити  дані та забезпечити розсилання змінних, необхідних для розрахунків кожному  процесу.

#include "stdafx.h"

#include "mpi.h"

#include "stdio.h"

#include <stdlib.h>

#define count 1024

#define fileIn "D:/temp/tmp/ListNoZero.txt"

#define fileOut "D:/temp/tmp/result.txt" 
 

bool WriteFile(float* a, int n, char* s);   // Функція зчитування даних з файлу в масив.

bool ReadFile(float* a, int n, char* s);   // Запис масиву в файл. 

int main(int argc, char* argv[])

{

      int p;   //  кількість процесів

     int rc;

     rank;  //  ранг процесу

      rc = MPI_Init(&argc,&argv);

      if (rc != MPI_SUCCESS)

      {

            printf ("Error starting MPI program. Terminating.\n");

            MPI_Abort(MPI_COMM_WORLD, rc);

      }

      MPI_Comm_size(MPI_COMM_WORLD, &p);

      MPI_Comm_rank(MPI_COMM_WORLD, &rank); 

      float* matrix = 0;

      float* bVector = 0; 

      //  Завантаження даних в майстер-процесі

      if (rank == 0)

      {

            matrix = (float*)malloc(count * count * sizeof(float));

            if (matrix == NULL){ MPI_Finalize(); printf ("not enough memory\n"); return -1; }

            bVector = (float*)malloc(count * sizeof(float));

            if (bVector == NULL){ MPI_Finalize(); printf ("not enough memory\n"); return -1; }

            ReadFile(matrix, count * count, fileIn);

            ReadFile(bVector, count, fileIn);

      } 

//  розрахунок параметрів розгалуження

      int *displsA = (int*)malloc(p * sizeof(int));

      int *scountsA = (int*)malloc(p * sizeof(int));

      int *displsB = (int*)malloc(p * sizeof(int));

      int *scountsB = (int*)malloc(p * sizeof(int));

      int h, block = count / p;

      for (h = 0; h < p; h++)

      {

            displsA[h] = h * block * count;

            scountsA[h] = block * count;

            displsB[h] = block * h;

            scountsB[h] = block;

      }

      scountsA[h - 1] += (count % p) * count;

      scountsB[h - 1] += (count % p);

      block = (rank != (p - 1) ? (count / p) : (count / p + (count % p))); 

//  виділення пам’яті на дані в процесах

// слід зауважити, що змінна block містить кількість рядків,

//  котрі було передано в даний процес

      float *procMatrix = (float*)malloc(block * count * sizeof(float));

      if (procMatrix == NULL){ MPI_Finalize(); printf ("not enough memory\n"); return -1; }

      float *procVector = (float*)malloc(block * count * sizeof(float));

      if (procVector == NULL){ MPI_Finalize(); printf ("not enough memory\n"); return -1; } 

//  розсилання даних по процесам

      MPI_Scatterv(matrix, scountsA, displsA, MPI_FLOAT, procMatrix, block * count, MPI_FLOAT, 0, MPI_COMM_WORLD);

      MPI_Scatterv(bVector, scountsB, displsB, MPI_FLOAT, procVector, block, MPI_FLOAT, 0, MPI_COMM_WORLD); 

//  звільнення пам’яті майстер-процесі

      if (rank == 0) {

            free(matrix); free(bVector);

             free(scountsA); free(displsA);

            free(scountsB); free(displsB);

      } 

//  прямий хід

      float* line;

      float* msLine;

//  формуємо масив, котрий буде приймати рядок, за допомоги якого ми будемо

//  вилучати невідомий в кожному блоці даних процесу.

      msLine = line = (float*)malloc(count * sizeof(float));

      if (msLine == NULL) { MPI_Finalize(); printf ("not enough memory\n"); return -1; } 

      int n, index = 0, w = count / p; float a, y = 0;

      for (int i = 0; i < count - 1; i++){

//  визначення з якого процесу і який рядок треба передати всім процесам

//  для вилучення невідомого

            if (i / w == rank){

                  // якщо це процес з якого передають данні

                  msLine = &procMatrix[count * (i % w)];

                  y = procVector[i % w];

                  index = rank;

            }

            else

            {

                  //  перевірка на обробок залишка

                  if (i / w == p)

                  {

                        if (rank == p - 1)

                        {

                              // якщо це процес з якого передають данні

                              msLine = &procMatrix[count * (i % w + w)];

                              y = procVector[i % w + w];

                        }

                        else

                        {

                              msLine = line;

                        }

                        index = p - 1;

                  }

                  else

                  {

                        //  встановлення вказівників на приймач

                        msLine = line;

                        index = i / w;

                  }

            } 

//  передача даних від одного процеса всім

      MPI_Bcast(msLine, count, MPI_FLOAT, index, MPI_COMM_WORLD);

      MPI_Bcast(&y, 1, MPI_FLOAT, index, MPI_COMM_WORLD); 

//  обробка прийнятих даних

      for (int j = 0; j < block; j++){

            if (count / p * rank + j > i) {

                  n = j * count;

                  a = procMatrix[n + i];

                  for (int k = i; k < count; k++){

                        procMatrix[n + k] = procMatrix[n + k] / a * msLine[i] - msLine[k]; }

                        procVector[j] = procVector[j] / a * msLine[i] - y; }

     }

//  закінчення прямого ходу

//  залишаємо масив лише на майстер-процесі, де ми будемо зберігати результат

//  на інших процесах даний масив вже не знадобиться

      if (rank)

      {

            free(line);

      } 
 

//  зворотній хід

      for (int i = count - 1; i >= 0; i--)

      {

            if (i / w == rank)

            {

                  y = procVector[i % w];

                  n = (i % w) * count;

                  for (int j = i + 1; j < count; j++)

                  {

                        y -= procMatrix[n + j];

                  }

                  y /= procMatrix[n + i];

                  index = rank;

            }

            else

            {

                  if (i / w == p)

                  {

                        if (rank == p - 1)

                        {

                              y = procVector[i % w + w];

                              n = (i % w + w) * count;

                              for (int j = i + 1; j < count; j++)

                              {

                                    y -= procMatrix[n + j];

                              }

                              y /= procMatrix[n + i];

                        }

                        index = p - 1;

                  }

                  else

                  {

                        index = i / w;

                  }

            } 

      //  розсилання знайденого невідомого всім процесам

            MPI_Bcast(&y, 1, MPI_FLOAT, index, MPI_COMM_WORLD); 
 

            if (!rank) { //  збереження значення невідомого

                  line[i] = y; } 

//  «запис» невідомого в масив