SHARE
TWEET

Untitled

a guest Nov 14th, 2019 114 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /* 10.  Вычисление обратной матрицы для заданной матрицы*/
  2.  
  3.  
  4. #include <cuda.h>
  5. #include <cuda_runtime.h>
  6. #include "device_launch_parameters.h"
  7. #include <time.h>
  8. #include <iostream>
  9. #include <fstream>
  10. #include <stdlib.h>
  11. #include <stdio.h>
  12. using namespace std;
  13.  
  14. double *createMass1(double *mass, int x, int y, int N) //создание подматрицы путем вычеркивания столбца и строки
  15. {
  16.      double *mass1 = new double[N*N];
  17.     int mi = 0;
  18.     for (int i = 0; i < N*N; i++, mi++)
  19.     {
  20.         while (mi / (N + 1) == x || mi % (N + 1) == y)
  21.             mi++;
  22.         mass1[i] = mass[mi];
  23.     }
  24.     return mass1;
  25. }
  26.  
  27. __device__ double* createMass( double *mass, int x, int y, int N)
  28. {
  29.     double *mass1 = new double[N*N];
  30.     int mi = 0;
  31.     for (int i = 0; i < N*N; i++, mi++)
  32.     {
  33.         while (mi / (N + 1) == x || mi % (N + 1) == y)
  34.             mi++;
  35.         mass1[i] = mass[mi];
  36.     }
  37.     return mass1;
  38. }
  39.  
  40. double determinant1( double *mass, int N)
  41. {
  42.      double determ = 1;
  43.     for (int i = 0; i < N*N; i += (N + 1)) //идем по диагонали
  44.     {
  45.         bool d = false;
  46.         for (int j = i; j < N*N; j += i%N)
  47.         {
  48.              double first = mass[j]; // Выбираем первый слева столбец матрицы
  49.             determ *= first;
  50.             do
  51.             {
  52.                 mass[j] /= first; //Все элементы первой строки делят на верхний элемент выбранного столбца
  53.                 if (d)
  54.                     mass[j] -= mass[i + (j%N - i%N)]; //вычеркивается первая строка и первый столбец
  55.                 j++;
  56.  
  57.             } while (j % N);
  58.             d = true;
  59.         }
  60.     }
  61.     return determ;
  62. }
  63.  
  64. __device__ double determinant(double *mass, int N)
  65. {
  66.      double determ = 1;
  67.     for (int i = 0; i < N*N; i += (N + 1)) //идем по диагонали
  68.     {
  69.         bool d = false;
  70.         for (int j = i; j < N*N; j += i%N) // от текущего элемента главной диагонали до последней строки
  71.         {
  72.              double first = mass[j];// Выбираем первый слева столбец матрицы
  73.             determ *= first;
  74.             do
  75.             {
  76.                 mass[j] /= first; //Все элементы первой строки делят на верхний элемент выбранного столбца
  77.                 if (d)
  78.                     mass[j] -= mass[i + (j%N - i%N)]; //вычеркивается первая строка и первый столбец
  79.                 j++;
  80.  
  81.             } while (j % N);
  82.             d = true;
  83.         }
  84.     }
  85.     return determ;
  86. }
  87.  
  88. __global__ void division(int i, double m,  double *mass)
  89. {
  90.     mass[threadIdx.x + i] /= m;  
  91.  
  92. }
  93.  
  94. __global__ void subtraction(int i, int j,  double *mass)
  95. {
  96.     mass[threadIdx.x + i] -= mass[threadIdx.x + j];
  97. }
  98.  
  99.  double deter(int N,  double *mass)
  100. {
  101.      double *mass1 = 0;
  102.     // Выделение памяти GPU
  103.     cudaMalloc((void**)&mass1, N * N * sizeof(double));
  104.     // Скопировать входные данные в GPU для обработки
  105.     cudaMemcpy(mass1, mass, N * N * sizeof( double), cudaMemcpyHostToDevice);
  106.      double determ = 1;
  107.     for (int i = 0; i < N*N; i += (N + 1))
  108.     {
  109.         bool d = false;
  110.         int j, k;
  111.         k = j = i;
  112.         do
  113.         {
  114.             determ *= mass[k]; // Выбираем первый слева столбец матрицы
  115.             division <<<1, (N - i%N) >>>((i + N*(j - i)), mass[k], mass1);
  116.             if (d)
  117.                 subtraction <<<1, (N - i%N) >>>((i + N*(j - i)), i, mass1);
  118.             d = true;
  119.             j++;
  120.             k += N;
  121.  
  122.         } while (j % N);
  123.         // Переписать результаты работы GPU в память CPU:
  124.         cudaMemcpy(mass, mass1, N * N * sizeof( double), cudaMemcpyDeviceToHost);
  125.     }
  126.     return determ;
  127. }
  128.  
  129. void start1( double *mass, double *temp, int N,  double *d)
  130. {
  131.      double det = determinant1(d, N);//вычисление детерминанта исходной матрицы
  132.     if (det == 0)
  133.     {
  134.         cout << "Obratnoy matrici ne suchestvuet" << endl;
  135.         return;
  136.     }
  137.     for (int p = 0; p < N*N; p++)
  138.     {
  139.         int i = p;
  140.         int j = i / N;
  141.         i = i % N;
  142.         //создание подматрицы путем вычеркивания i-ого и j-ого столбца для нахождения минора
  143.          double* A = createMass1(mass, i, j, N - 1);
  144.         //нахождение алгебраического дополнения и деление его на определитель
  145.         temp[p] = (determinant1(A, N - 1) * ((2 + i + j) % 2 ? -1 : 1)) / det;
  146.         delete[]A;
  147.     }
  148. }
  149.  
  150. __global__ void transp(double *mass, double det, int N)
  151. {
  152.     int k = blockIdx.x * blockDim.x + threadIdx.x;
  153.     mass[k] *= ((2 + (k % N) + (k / N)) % 2 ? -1 : 1);
  154.     mass[k] /= det;
  155. }
  156.  
  157. __global__ void comp( double *mass, double *temp, int N, double det)
  158. {
  159.     int i = threadIdx.x%N;
  160.     int j = threadIdx.x / N;
  161.     double* A = createMass(mass, i, j, N - 1);
  162.     temp[threadIdx.x] = determinant(A, N - 1);
  163.     delete[]A;
  164. }
  165.  
  166.  
  167. void start( double *mass, double *temp, int N, double *d)
  168. {
  169.     size_t free_before, free_after, total;
  170.  
  171.     cudaMemGetInfo(&free_before,&total);
  172.     cout << ("%d KB free_before of total %d KB\n",free_before/1024,total/1024);
  173.  
  174.     if (N >= 32)
  175.     {
  176.         cout << "Matrica vichodit za granici pamyati GPU adaptera" << endl;
  177.         return;
  178.     }
  179.  
  180.     double det = deter(N, d);
  181.     if (det == 0)
  182.     {
  183.         cout << "Obratnoy matrici ne suchestvuet" << endl;
  184.         return;
  185.     }
  186.  
  187.     double *temp1 = 0, *mass1 = 0;
  188.     // Выделение памяти GPU
  189.     cudaMalloc((void**)&temp1, N * N * sizeof(double));
  190.     cudaMalloc((void**)&mass1, N * N * sizeof(double));
  191.  
  192.     // Скопировать входные данные в GPU для обработки
  193.     cudaMemcpy(mass1, mass, N * N * sizeof(double), cudaMemcpyHostToDevice);
  194.    
  195.     // Инициализация ядра в GPU
  196.     comp<<< 1, N*N >>>(mass1, temp1, N, det);
  197.     transp<<< 1, N*N >>>(temp1, det, N);
  198.  
  199.     // Переписать результаты работы GPU в память CPU:
  200.     cudaMemcpy(temp, temp1, N * N * sizeof(double), cudaMemcpyDeviceToHost);
  201.  
  202.     cudaMemGetInfo(&free_after,&total);
  203.     cout << ("%d KB free_after of total %d KB\n",free_after/1024,total/1024);
  204.  
  205.     // Выделено в памяти 2 массива
  206.     int count_of_threads = (free_after - free_before) /  (2 * N * N * sizeof(double));  
  207.  
  208.     cout << "\n" << count_of_threads << endl;
  209.  
  210.     cudaFree(temp1);
  211.     cudaFree(mass1);
  212. }
  213.  
  214. bool compare(double *mass1, double *mass2, int N)
  215. {
  216.     for (int i = 0; i < N; i++)
  217.     {
  218.         if (mass1[i] != mass2[i])
  219.             return false;
  220.     }
  221.     return true;
  222. }
  223.  
  224. int main()
  225. {
  226.     setlocale(LC_CTYPE, "rus");
  227.     srand(time(0));
  228.     int N;
  229.     cout << "Vvedite razmer matrici: ";
  230.     cin >> N;
  231.     double *mass = new double[N*N];
  232.     double *temp = new double[N*N];
  233.     double *_mass = new double[N*N];
  234.     double *_temp = new double[N*N];
  235.     double *d = new double[N*N];
  236.     double *_d = new double[N*N];
  237.  
  238.     for (int i = 0; i < N*N; i++)
  239.         mass[i] = _mass[i] = d[i] = _d[i] = ((double)rand()) / 100000;
  240.  
  241.  
  242.  
  243.     double t1 = clock();
  244.     start(mass, temp, N, d);
  245.     t1 = clock() - t1;
  246.     cout << endl << endl << "Parallel: " << t1 / CLOCKS_PER_SEC << " sec." << endl;
  247.  
  248.     double t = clock();
  249.     start1(_mass, _temp, N, _d);
  250.     t = clock() - t;
  251.     cout << endl << "Posledovat: " << t / CLOCKS_PER_SEC << " sec." << endl;
  252.  
  253.     if (t - t1 > 0)
  254.         cout << endl << endl << "Parallel bystree posledovat v " << int(t / t1) << " raz." << endl;
  255.     else
  256.         cout << endl << endl << "Parallel medlennee posledovat v " << int(t1 / t) << " raz." << endl;
  257.  
  258.     if (compare(temp, _temp, N))
  259.         cout << endl << "Resultat identichni!" << endl << endl;
  260.     else
  261.         cout << endl << "Resultat razlichni!" << endl << endl;
  262.  
  263.     delete[]temp;
  264.     delete[]mass;
  265.     delete[]_temp;
  266.     delete[]_mass;
  267. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top