Facebook
From Sludgy Agouti, 6 Years ago, written in C++.
Embed
Download Paste or View Raw
Hits: 258
  1. // classic.cpp : "Textbook" implementation of matrix multiply
  2.  
  3. // Author:  Paul J. Drongowski
  4. // Address: Boston Design Center
  5. //          Advanced Micro Devices, Inc.
  6. //          Boxborough, MA 01719
  7. // Date:    20 October 2005
  8. //
  9. // Copyright (c) 2005 Advanced Micro Devices, Inc.
  10.  
  11. // Celem tego programu jest prezentacja pomiaru i analizy
  12. //efektywnosci programu za pomocą  CodeAnalyst(tm).
  13. // Implementacja mnożenia macierzy jest realizowana za pomoca typowego
  14. // algorytmu podręcznikowego.
  15. //  Efektywność tego podejscia jest niska poprzez
  16. // nieefektywną  kolejnosć odwołań do elementów macierzy.
  17. #include <stdio.h>
  18. #include <time.h>
  19. #include <windows.h>
  20. #include "omp.h"
  21. #include <string>
  22.  
  23. #define USE_MULTIPLE_THREADS true
  24. #define MAXTHREADS 128
  25. int NumThreads;
  26. double start;
  27.  
  28. static const int ROWS = 2048;     // liczba wierszy macierzy
  29. static const int COLUMNS = 2048;  // lizba kolumn macierzy
  30. int r = 256;
  31.  
  32. float matrix_a[ROWS][COLUMNS];    // lewy operand
  33. float matrix_b[ROWS][COLUMNS];    // prawy operand
  34. float matrix_r[ROWS][COLUMNS];    // wynik
  35.  
  36. float result = 0;
  37. float resultParallel = 0;
  38. float result6 = 0;
  39. float resultParallel6 = 0;
  40.  
  41. float minResult = 100;
  42. float minResultParallel = 100;
  43. float minResult6 = 100;
  44. float minResultParallel6 = 100;
  45.  
  46. float maxResult = 0;
  47. float maxResultParallel = 0;
  48. float maxResult6 = 0;
  49. float maxResultParallel6 = 0;
  50.  
  51. float tmp;
  52.  
  53. FILE *result_file;
  54.  
  55. void initialize_matrices()
  56. {
  57.         // zdefiniowanie zawarosci poczatkowej macierzy
  58. #pragma omp parallel for
  59.         for (int i = 0; i < ROWS; i++) {
  60.                 for (int j = 0; j < COLUMNS; j++) {
  61.                         matrix_a[i][j] = (float)rand() / RAND_MAX;
  62.                         matrix_b[i][j] = (float)rand() / RAND_MAX;
  63.                         matrix_r[i][j] = 0.0;
  64.                 }
  65.         }
  66. }
  67.  
  68. void initialize_matricesZ()
  69. {
  70.         // zdefiniowanie zawarosci poczatkowej macierzy
  71.  
  72. #pragma omp parallel for
  73.         for (int i = 0; i < ROWS; i++) {
  74.                 for (int j = 0; j < COLUMNS; j++) {
  75.                         matrix_r[i][j] = 0.0;
  76.                 }
  77.         }
  78. }
  79.  
  80. void print_result()
  81. {
  82.         // wydruk wyniku
  83.         for (int i = 0; i < ROWS; i++) {
  84.                 for (int j = 0; j < COLUMNS; j++) {
  85.                         fprintf(result_file, "%6.4f ", matrix_r[i][j]);
  86.                 }
  87.                 fprintf(result_file, "\n");
  88.         }
  89. }
  90.  
  91. void multiply_matrices_IKJ()
  92. {
  93.         for (int i = 0; i < ROWS; i++)
  94.                 for (int k = 0; k < COLUMNS; k++)
  95.                         for (int j = 0; j < COLUMNS; j++)
  96.                                 matrix_r[i][j] += matrix_a[i][k] * matrix_b[k][j];
  97. }
  98.  
  99. void multiply_matrices_IKJ_Parallel()
  100. {
  101.  
  102. #pragma omp parallel for
  103.         for (int i = 0; i < ROWS; i++)
  104.                 for (int k = 0; k < COLUMNS; k++)
  105.                         for (int j = 0; j < COLUMNS; j++)
  106.                         {
  107.                                 matrix_r[i][j] += matrix_a[i][k] * matrix_b[k][j];
  108.                         }
  109.  
  110. }
  111.  
  112. void multiply_matrices_JKI()
  113. {
  114.         // mnozenie macierzy
  115.         for (int j = 0; j < COLUMNS; j++)
  116.                 for (int k = 0; k < COLUMNS; k++)
  117.                         for (int i = 0; i < ROWS; i++)
  118.                                 matrix_r[i][j] += matrix_a[i][k] * matrix_b[k][j];
  119.  
  120. }
  121.  
  122. void multiply_matrices_JKI_Parallel()
  123. {
  124.         // mnozenie macierzy
  125. #pragma omp parallel for
  126.         for (int j = 0; j < COLUMNS; j++)
  127.                 for (int k = 0; k < COLUMNS; k++)
  128.                         for (int i = 0; i < ROWS; i++)
  129.                                 matrix_r[i][j] += matrix_a[i][k] * matrix_b[k][j];
  130.  
  131. }
  132.  
  133. void multiply_matrices_IJK_IKJ()
  134. {
  135.  
  136.         for (int i = 0; i < ROWS; i += r)
  137.                 for (int j = 0; j < COLUMNS; j += r)
  138.                         for (int k = 0; k < COLUMNS; k += r)
  139.                                 for (int ii = i; ii < i + r; ii++)  // WYNIK CZĘŚCIOWY
  140.                                         for (int kk = k; kk < k + r; kk++)
  141.                                                 for (int jj = j; jj < j + r; jj++)
  142.                                                         matrix_r[ii][jj] += matrix_a[ii][kk] * matrix_b[kk][jj];
  143. }
  144.  
  145. void multiply_matrices_IJK_IKJ_Parallel()
  146. {
  147.  
  148. #pragma omp parallel for
  149.         for (int i = 0; i < ROWS; i += r)
  150.                 for (int j = 0; j < COLUMNS; j += r)
  151.                         for (int k = 0; k < COLUMNS; k += r)
  152.                                 for (int ii = i; ii < i + r; ii++)  // WYNIK CZĘŚCIOWY
  153.                                         for (int kk = k; kk < k + r; kk++)
  154.                                                 for (int jj = j; jj < j + r; jj++)
  155.                                                         matrix_r[ii][jj] += matrix_a[ii][kk] * matrix_b[kk][jj];
  156. }
  157.  
  158. float print_elapsed_time(std::string info)
  159. {
  160.         double elapsed;
  161.         double resolution;
  162.  
  163.         // wyznaczenie i zapisanie czasu przetwarzania
  164.         elapsed = (double)clock() / CLK_TCK;
  165.         resolution = 1.0 / CLK_TCK;
  166.  
  167.         fprintf(result_file,
  168.                 "Czas wykonania programu %s: %8.4f sec (%6.4f sec rozdzielczosc pomiaru)\n",
  169.                 info.c_str(), elapsed - start, resolution);
  170.         return elapsed - start;
  171. }
  172.  
  173. void enterLine()
  174. {
  175.         fprintf(result_file,
  176.                 "\n"
  177.         );
  178.         fprintf(result_file,
  179.                 "\n"
  180.         );
  181. }
  182.  
  183. void load_IKJ()
  184. {
  185.         for (int i = 0;i < 100;i++)
  186.         {
  187.                 initialize_matricesZ();
  188.                 start = (double)clock() / CLK_TCK;
  189.                 multiply_matrices_IKJ();
  190.                 tmp = print_elapsed_time("Sekwencyjnie IKJ ");
  191.                 result += tmp;
  192.  
  193.                 if (tmp > maxResult)
  194.                 {
  195.                         maxResult = tmp;
  196.                 }
  197.  
  198.                 if (tmp < minResult)
  199.                 {
  200.                         minResult = tmp;
  201.                 }
  202.         }
  203.         enterLine();
  204. }
  205.  
  206. void loadParallel_IKJ()
  207. {
  208.         for (int i = 0; i < 100; i++)
  209.         {
  210.                 initialize_matricesZ();
  211.                 start = (double)clock() / CLK_TCK;
  212.                 multiply_matrices_IKJ_Parallel();
  213.                 tmp = print_elapsed_time("Rownolegle IKJ ");
  214.                 resultParallel += tmp;
  215.  
  216.                 if (tmp > maxResultParallel)
  217.                 {
  218.                         maxResultParallel = tmp;
  219.                 }
  220.  
  221.                 if (tmp < minResultParallel)
  222.                 {
  223.                         minResultParallel = tmp;
  224.                 }
  225.         }
  226.         enterLine();
  227. }
  228.  
  229. void load_IJK_IKJ()
  230. {
  231.         for (int i = 0; i < 100; i++)
  232.         {
  233.                 initialize_matricesZ();
  234.                 start = (double)clock() / CLK_TCK;
  235.                 multiply_matrices_IJK_IKJ();
  236.                 tmp = print_elapsed_time("Sekwencyjnie IJK_IKJ ");
  237.                 result6 += tmp;
  238.  
  239.                 if (tmp > maxResult6)
  240.                 {
  241.                         maxResult6 = tmp;
  242.                 }
  243.  
  244.                 if (tmp < minResult6)
  245.                 {
  246.                         minResult6 = tmp;
  247.                 }
  248.         }
  249.         enterLine();
  250. }
  251.  
  252. void loadParallel_IJK_IKJ()
  253. {
  254.         for (int i = 0; i < 100; i++)
  255.         {
  256.                 initialize_matricesZ();
  257.                 start = (double)clock() / CLK_TCK;
  258.                 multiply_matrices_IJK_IKJ_Parallel();
  259.                 tmp = print_elapsed_time("Rownolegle IJK_IKJ");
  260.                 resultParallel6 += tmp;
  261.  
  262.                 if (tmp > maxResultParallel6)
  263.                 {
  264.                         maxResultParallel6 = tmp;
  265.                 }
  266.  
  267.                 if (tmp < minResultParallel6)
  268.                 {
  269.                         minResultParallel6 = tmp;
  270.                 }
  271.         }
  272.         enterLine();
  273. }
  274.  
  275. void save_IKJ()
  276. {
  277.         float avgResult = result / 100;
  278.         fprintf(result_file,
  279.                 "Średnia sekwencyjnego %8.4f \n", avgResult
  280.         );
  281.         fprintf(result_file,
  282.                 "Max sekwencyjnego %8.4f \n", maxResult
  283.         );
  284.         fprintf(result_file,
  285.                 "Min sekwencyjnego %8.4f \n", minResult
  286.         );
  287.         enterLine();
  288. }
  289.  
  290. void saveParallel_IKJ()
  291. {
  292.         float avgResultParallel = resultParallel / 100;
  293.         fprintf(result_file,
  294.                 "Średnia równoległego %8.4f \n", avgResultParallel
  295.         );
  296.         fprintf(result_file,
  297.                 "Min równoległego %8.4f \n", minResultParallel
  298.         );
  299.         fprintf(result_file,
  300.                 "Max równoległego %8.4f \n", maxResultParallel
  301.         );
  302.         fclose(result_file);
  303. }
  304.  
  305. void saveParallel_IJK_IKJ()
  306. {
  307.         float avgResultParallel = resultParallel6 / 100;
  308.         fprintf(result_file,
  309.                 "Średnia równoległego IJK_IKJ   %8.4f \n", avgResultParallel
  310.         );
  311.         fprintf(result_file,
  312.                 "Min równoległego IJK_IKJ   %8.4f \n", minResultParallel6
  313.         );
  314.         fprintf(result_file,
  315.                 "Max równoległego IJK_IKJ   %8.4f \n", maxResultParallel6
  316.         );
  317.         fclose(result_file);
  318. }
  319.  
  320. void save_IJK_IKJ()
  321. {
  322.         float avgResult = result6 / 100;
  323.         fprintf(result_file,
  324.                 "Średnia sekwencyjnego IJK_IKJ %8.4f \n", avgResult
  325.         );
  326.         fprintf(result_file,
  327.                 "Max sekwencyjnego IJK_IKJ %8.4f \n", maxResult6
  328.         );
  329.         fprintf(result_file,
  330.                 "Min sekwencyjnego IJK_IKJ %8.4f \n", minResult6
  331.         );
  332.         enterLine();
  333. }
  334.  
  335. int main(int argc, char* argv[])
  336. {
  337.         omp_set_num_threads(4);
  338.  
  339.         //       start = (double) clock() / CLK_TCK ;
  340.         if ((result_file = fopen("classic.txt", "w")) == NULL) {
  341.                 fprintf(stderr, "nie mozna otworzyc pliku wyniku \n");
  342.                 perror("classic");
  343.                 return(EXIT_FAILURE);
  344.         }
  345.  
  346.  
  347.         //Determine the number of threads to use
  348.         if (USE_MULTIPLE_THREADS) {
  349.                 SYSTEM_INFO SysInfo;
  350.                 GetSystemInfo(&SysInfo);
  351.                 NumThreads = SysInfo.dwNumberOfProcessors;
  352.                 if (NumThreads > MAXTHREADS)
  353.                         NumThreads = MAXTHREADS;
  354.         }
  355.         else
  356.                 NumThreads = 1;
  357.         fprintf(result_file, "Klasyczny algorytm mnozenia macierzy, liczba watkow %d \n", NumThreads);
  358.         printf("liczba watkow  = %d\n\n", NumThreads);
  359.  
  360.         //load_IKJ();
  361.         //loadParallel_IKJ();
  362.  
  363.         //save_IKJ();
  364.         //saveParallel_IKJ();
  365.  
  366.         //load_IJK_IKJ();
  367.         //loadParallel_IJK_IKJ();
  368.         //
  369.         //
  370.         //save_IJK_IKJ();
  371.         //saveParallel_IJK_IKJ();
  372.  
  373.         initialize_matricesZ();
  374.         start = (double)clock() / CLK_TCK;
  375.         multiply_matrices_IJK_IKJ_Parallel();
  376.         tmp = print_elapsed_time("Sekwencyjnie IJK_IKJ ");
  377.  
  378.         return(0);
  379. }