API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Apr 23rd, 2018
57
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 12.78 KB | None | 0 0
raw download clone embed print report
  1. ////////////////////////////////////////////////////////////////////////////////
  2. //  mvec.cpp
  3.  
  4. ////////////////////////////////////////////////////////////////////////////////
  5. //  mvec.cpp
  6.  
  7. #include "stdafx.h"
  8. #include<iostream>
  9. #include <emmintrin.h>
  10. #include <intrin.h>
  11. #include "inc_head.h"
  12.  
  13. using namespace std;
  14.  
  15. void matvec_orgin(double* a, double* x, double* y, int n)
  16. {
  17.     int i, j;
  18.  
  19.     for (i = 0; i < n; i++) {
  20.         y[i] = 0.0;
  21.         for (j = 0; j < n; j++) {
  22.             y[i] += a[i + n * j] * x[j];
  23.         }
  24.     }
  25. }
  26.  
  27.  
  28. void matvec_opt_1(double* a, double* x, double* y, int n)
  29. /*=========================================================================
  30. Usuniecie skokow w danych
  31. ===========================================================================*/
  32. {
  33.     int i, j;
  34.     register double r;
  35.     memset(y, 0, n*sizeof(double));
  36.  
  37.     for (j = 0; j < n; j++) {
  38.         r = x[j];
  39.  
  40.         for (i = 0; i < n; i++) {
  41.             y[i] += a[i + n * j] * r ;
  42.         }
  43.     }
  44. }
  45.  
  46.  
  47.  
  48.  
  49.  
  50. void matvec_opt_2(double* a, double* x, double* y, int n)
  51. /*============================================================================
  52. Rozwijanie petli
  53. =============================================================================*/
  54. {
  55.  
  56.     int i, j;
  57.     int rest = n % 6;
  58.     register double r;
  59.     memset((double*)y, 0, sizeof(double)*n);
  60.  
  61.     for (j = 0; j < n; j++) {
  62.  
  63.         for (i = 0; i < rest; i++) {
  64.             y[i] += a[i + n * j];
  65.         }
  66.  
  67.         for (i = rest; i < n; i+=6) {
  68.             y[i] += a[i + n * j] ;
  69.             y[i+1] += a[i+1 + n * j] ;
  70.             y[i+2] += a[i+2 + n * j] ;
  71.             y[i+3] += a[i+3 + n * j] ;
  72.             y[i+4] += a[i+4 + n * j] ;
  73.             y[i+5] += a[i+5 + n * j] ;
  74.         }
  75.     }
  76.  
  77.    
  78.  
  79. }
  80.  
  81. void matvec_opt_3(double* a, double* x, double* y, int n, int lb)
  82. {
  83.     int i, j;
  84.  
  85.     memset((void *)y, 0, n * sizeof(double));
  86.  
  87.     __m128d ry0, ry1, ry2, ry3, ra0, ra1, ra2, ra3, rx0;
  88.     double *ptr_a, *ptr_x, *ptr_y;
  89.     const int mr = 8;
  90.     if (mr != lb)
  91.     {
  92.         cout << "lb != mr\n";
  93.         system("pause");
  94.         exit(1);
  95.     }
  96.  
  97.     ptr_a = a;
  98.  
  99.     for (i = 0; i < n; i += mr)
  100.     {
  101.         ry0 = _mm_setzero_pd();
  102.         ry1 = _mm_setzero_pd();
  103.         ry2 = _mm_setzero_pd();
  104.         ry3 = _mm_setzero_pd();
  105.  
  106.         //ptr_a = &a[n*i];
  107.         ptr_y = &y[i];
  108.         ptr_x = x;
  109.  
  110.         for (j = 0; j < n; j++)
  111.         {
  112.             //_mm_prefetch((const char *)(ptr_a+8), _MM_HINT_T0);
  113.             _mm_prefetch((const char *)(ptr_a + 8), _MM_HINT_NTA);
  114.             if (j % mr == 0)
  115.                 _mm_prefetch((const char *)(ptr_x + 8), _MM_HINT_T0);
  116.             rx0 = _mm_load1_pd(ptr_x);
  117.             ra0 = _mm_load_pd(ptr_a);
  118.             ra1 = _mm_load_pd(ptr_a + 2);
  119.             ra2 = _mm_load_pd(ptr_a + 4);
  120.             ra3 = _mm_load_pd(ptr_a + 6);
  121.  
  122.             ptr_a += mr;
  123.             ptr_x++;
  124.  
  125.             ra0 = _mm_mul_pd(ra0, rx0);
  126.             ra1 = _mm_mul_pd(ra1, rx0);
  127.             ra2 = _mm_mul_pd(ra2, rx0);
  128.             ra3 = _mm_mul_pd(ra3, rx0);
  129.  
  130.             ry0 = _mm_add_pd(ry0, ra0);
  131.             ry1 = _mm_add_pd(ry1, ra1);
  132.             ry2 = _mm_add_pd(ry2, ra2);
  133.             ry3 = _mm_add_pd(ry3, ra3);
  134.         }
  135.  
  136.         _mm_store_pd(ptr_y, ry0);
  137.         _mm_store_pd(ptr_y + 2, ry1);
  138.         _mm_store_pd(ptr_y + 4, ry2);
  139.         _mm_store_pd(ptr_y + 6, ry3);
  140.     }
  141. }
  142.  
  143. void matvec_opt_4(double* a, double* x, double* y, int n, int lb)
  144. {
  145.     int i, j;
  146.  
  147.     memset((void *)y, 0, n * sizeof(double));
  148.  
  149.     __m128d ry0, ry1, ry2, ry3, ra0, ra1, ra2, ra3, rx0;
  150.     double *ptr_a, *ptr_x, *ptr_y;
  151.     const int mr = 8;
  152.     const int nr = 4;
  153.     if (mr != lb || nr != lb / 2)
  154.     {
  155.         cout << "lb != mr || lb != nr\n";
  156.         system("pause");
  157.         exit(1);
  158.     }
  159.  
  160.     ptr_a = a;
  161.  
  162.     for (i = 0; i < n; i += mr)
  163.     {
  164.         ry0 = _mm_setzero_pd();
  165.         ry1 = _mm_setzero_pd();
  166.         ry2 = _mm_setzero_pd();
  167.         ry3 = _mm_setzero_pd();
  168.  
  169.         //ptr_a = &a[n*i];
  170.         ptr_y = &y[i];
  171.         ptr_x = x;
  172.  
  173.         for (j = 0; j < n; j += nr)
  174.         {
  175.             _mm_prefetch((const char *)(ptr_a + mr * nr), _MM_HINT_NTA);
  176.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 1)), _MM_HINT_NTA);
  177.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 2)), _MM_HINT_NTA);
  178.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 3)), _MM_HINT_NTA);
  179.  
  180.             if (j%nr == 0)
  181.                 _mm_prefetch((const char *)(ptr_x + nr), _MM_HINT_T0);
  182.             //--------------------------0
  183.             rx0 = _mm_load1_pd(ptr_x);
  184.             ra0 = _mm_load_pd(ptr_a);
  185.             ra1 = _mm_load_pd(ptr_a + 2);
  186.             ra2 = _mm_load_pd(ptr_a + 4);
  187.             ra3 = _mm_load_pd(ptr_a + 6);
  188.  
  189.             ra0 = _mm_mul_pd(ra0, rx0);
  190.             ra1 = _mm_mul_pd(ra1, rx0);
  191.             ra2 = _mm_mul_pd(ra2, rx0);
  192.             ra3 = _mm_mul_pd(ra3, rx0);
  193.  
  194.             ry0 = _mm_add_pd(ry0, ra0);
  195.             ry1 = _mm_add_pd(ry1, ra1);
  196.             ry2 = _mm_add_pd(ry2, ra2);
  197.             ry3 = _mm_add_pd(ry3, ra3);
  198.  
  199.             //--------------------------1
  200.             rx0 = _mm_load1_pd(ptr_x + 1);
  201.             ra0 = _mm_load_pd(ptr_a + 8);
  202.             ra1 = _mm_load_pd(ptr_a + 10);
  203.             ra2 = _mm_load_pd(ptr_a + 12);
  204.             ra3 = _mm_load_pd(ptr_a + 14);
  205.  
  206.             ra0 = _mm_mul_pd(ra0, rx0);
  207.             ra1 = _mm_mul_pd(ra1, rx0);
  208.             ra2 = _mm_mul_pd(ra2, rx0);
  209.             ra3 = _mm_mul_pd(ra3, rx0);
  210.  
  211.             ry0 = _mm_add_pd(ry0, ra0);
  212.             ry1 = _mm_add_pd(ry1, ra1);
  213.             ry2 = _mm_add_pd(ry2, ra2);
  214.             ry3 = _mm_add_pd(ry3, ra3);
  215.  
  216.             //--------------------------2
  217.             rx0 = _mm_load1_pd(ptr_x + 2);
  218.             ra0 = _mm_load_pd(ptr_a + 16);
  219.             ra1 = _mm_load_pd(ptr_a + 18);
  220.             ra2 = _mm_load_pd(ptr_a + 20);
  221.             ra3 = _mm_load_pd(ptr_a + 22);
  222.  
  223.             ra0 = _mm_mul_pd(ra0, rx0);
  224.             ra1 = _mm_mul_pd(ra1, rx0);
  225.             ra2 = _mm_mul_pd(ra2, rx0);
  226.             ra3 = _mm_mul_pd(ra3, rx0);
  227.  
  228.             ry0 = _mm_add_pd(ry0, ra0);
  229.             ry1 = _mm_add_pd(ry1, ra1);
  230.             ry2 = _mm_add_pd(ry2, ra2);
  231.             ry3 = _mm_add_pd(ry3, ra3);
  232.  
  233.             //--------------------------3
  234.             rx0 = _mm_load1_pd(ptr_x + 3);
  235.             ra0 = _mm_load_pd(ptr_a + 24);
  236.             ra1 = _mm_load_pd(ptr_a + 26);
  237.             ra2 = _mm_load_pd(ptr_a + 28);
  238.             ra3 = _mm_load_pd(ptr_a + 30);
  239.  
  240.             ptr_a += nr * mr;
  241.             ptr_x += nr;
  242.  
  243.             ra0 = _mm_mul_pd(ra0, rx0);
  244.             ra1 = _mm_mul_pd(ra1, rx0);
  245.             ra2 = _mm_mul_pd(ra2, rx0);
  246.             ra3 = _mm_mul_pd(ra3, rx0);
  247.  
  248.             ry0 = _mm_add_pd(ry0, ra0);
  249.             ry1 = _mm_add_pd(ry1, ra1);
  250.             ry2 = _mm_add_pd(ry2, ra2);
  251.             ry3 = _mm_add_pd(ry3, ra3);
  252.         }
  253.  
  254.         _mm_store_pd(ptr_y, ry0);
  255.         _mm_store_pd(ptr_y + 2, ry1);
  256.         _mm_store_pd(ptr_y + 4, ry2);
  257.         _mm_store_pd(ptr_y + 6, ry3);
  258.     }
  259. }
  260.  
  261. #ifdef YES_AVX
  262. void matvec_opt_5(double* a, double* x, double* y, int n, int lb)
  263. {
  264.     int i, j;
  265.  
  266.     memset((void *)y, 0, n * sizeof(double));
  267.  
  268.     __m256d ry0, ry1, ry2, ry3, ra0, ra1, ra2, ra3, rx0, rx1;
  269.     double *ptr_a, *ptr_x, *ptr_y;
  270.     const int mr = 8, nr = 8;
  271.     if (mr != lb)
  272.     {
  273.         cout << "lb != mr\n";
  274.         system("pause");
  275.         exit(1);
  276.     }
  277.  
  278.     ptr_a = a;
  279.  
  280.     for (i = 0; i < n; i += mr)
  281.     {
  282.         ry0 = _mm256_setzero_pd();
  283.         ry1 = ry2 = ry3 = ry0;
  284.  
  285.         //ptr_a = &a[n*i];
  286.         ptr_y = &y[i];
  287.         ptr_x = x;
  288.  
  289.         for (j = 0; j < n; j += nr)
  290.         {
  291.             _mm_prefetch((const char *)(ptr_a + mr * nr), _MM_HINT_NTA);
  292.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 1)), _MM_HINT_NTA);
  293.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 2)), _MM_HINT_NTA);
  294.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 3)), _MM_HINT_NTA);
  295.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 4)), _MM_HINT_NTA);
  296.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 5)), _MM_HINT_NTA);
  297.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 6)), _MM_HINT_NTA);
  298.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 7)), _MM_HINT_NTA);
  299.  
  300.             _mm_prefetch((const char *)(ptr_x + nr), _MM_HINT_T0);
  301.  
  302.             //---------------------------------------------------0
  303.             rx0 = _mm256_broadcast_sd(ptr_x);
  304.             rx1 = _mm256_broadcast_sd(ptr_x + 1);
  305.             ra0 = _mm256_load_pd(ptr_a);
  306.             ra1 = _mm256_load_pd(ptr_a + 4);
  307.             ra2 = _mm256_load_pd(ptr_a + mr);
  308.             ra3 = _mm256_load_pd(ptr_a + mr + 4);
  309.  
  310.             ra0 = _mm256_mul_pd(ra0, rx0);
  311.             ra1 = _mm256_mul_pd(ra1, rx0);
  312.             ra2 = _mm256_mul_pd(ra2, rx1);
  313.             ra3 = _mm256_mul_pd(ra3, rx1);
  314.  
  315.             ry0 = _mm256_add_pd(ry0, ra0);
  316.             ry1 = _mm256_add_pd(ry1, ra1);
  317.             ry2 = _mm256_add_pd(ry2, ra2);
  318.             ry3 = _mm256_add_pd(ry3, ra3);
  319.  
  320.             //---------------------------------------------------1
  321.  
  322.             rx0 = _mm256_broadcast_sd(ptr_x + 2);
  323.             rx1 = _mm256_broadcast_sd(ptr_x + 3);
  324.             ra0 = _mm256_load_pd(ptr_a + 2 * mr);
  325.             ra1 = _mm256_load_pd(ptr_a + 2 * mr + 4);
  326.             ra2 = _mm256_load_pd(ptr_a + 3 * mr);
  327.             ra3 = _mm256_load_pd(ptr_a + 3 * mr + 4);
  328.  
  329.             ra0 = _mm256_mul_pd(ra0, rx0);
  330.             ra1 = _mm256_mul_pd(ra1, rx0);
  331.             ra2 = _mm256_mul_pd(ra2, rx1);
  332.             ra3 = _mm256_mul_pd(ra3, rx1);
  333.  
  334.             ry0 = _mm256_add_pd(ry0, ra0);
  335.             ry1 = _mm256_add_pd(ry1, ra1);
  336.             ry2 = _mm256_add_pd(ry2, ra2);
  337.             ry3 = _mm256_add_pd(ry3, ra3);
  338.  
  339.             //---------------------------------------------------2
  340.             rx0 = _mm256_broadcast_sd(ptr_x + 4);
  341.             rx1 = _mm256_broadcast_sd(ptr_x + 5);
  342.             ra0 = _mm256_load_pd(ptr_a + 4 * mr);
  343.             ra1 = _mm256_load_pd(ptr_a + 4 * mr + 4);
  344.             ra2 = _mm256_load_pd(ptr_a + 5 * mr);
  345.             ra3 = _mm256_load_pd(ptr_a + 5 * mr + 4);
  346.  
  347.             ra0 = _mm256_mul_pd(ra0, rx0);
  348.             ra1 = _mm256_mul_pd(ra1, rx0);
  349.             ra2 = _mm256_mul_pd(ra2, rx1);
  350.             ra3 = _mm256_mul_pd(ra3, rx1);
  351.  
  352.             ry0 = _mm256_add_pd(ry0, ra0);
  353.             ry1 = _mm256_add_pd(ry1, ra1);
  354.             ry2 = _mm256_add_pd(ry2, ra2);
  355.             ry3 = _mm256_add_pd(ry3, ra3);
  356.  
  357.             //---------------------------------------------------3
  358.             rx0 = _mm256_broadcast_sd(ptr_x + 6);
  359.             rx1 = _mm256_broadcast_sd(ptr_x + 7);
  360.             ra0 = _mm256_load_pd(ptr_a + 6 * mr);
  361.             ra1 = _mm256_load_pd(ptr_a + 6 * mr + 4);
  362.             ra2 = _mm256_load_pd(ptr_a + 7 * mr);
  363.             ra3 = _mm256_load_pd(ptr_a + 7 * mr + 4);
  364.  
  365.             ptr_a += mr * nr;
  366.             ptr_x += nr;
  367.  
  368.             ra0 = _mm256_mul_pd(ra0, rx0);
  369.             ra1 = _mm256_mul_pd(ra1, rx0);
  370.             ra2 = _mm256_mul_pd(ra2, rx1);
  371.             ra3 = _mm256_mul_pd(ra3, rx1);
  372.  
  373.             ry0 = _mm256_add_pd(ry0, ra0);
  374.             ry1 = _mm256_add_pd(ry1, ra1);
  375.             ry2 = _mm256_add_pd(ry2, ra2);
  376.             ry3 = _mm256_add_pd(ry3, ra3);
  377.         }
  378.  
  379.         ry0 = _mm256_add_pd(ry0, ry2);
  380.         ry1 = _mm256_add_pd(ry1, ry3);
  381.         _mm256_store_pd(ptr_y, ry0);
  382.         _mm256_store_pd(ptr_y + 4, ry1);
  383.     }
  384. }
  385. #endif
  386.  
  387.  
  388. #ifdef YES_FMA
  389. void matvec_opt_6(double* a, double* x, double* y, int n, int lb)
  390. {
  391.     int i, j;
  392.  
  393.     memset((void *)y, 0, n * sizeof(double));
  394.  
  395.     __m256d ry0, ry1, ry2, ry3, ra0, ra1, ra2, ra3, rx0, rx1;
  396.     double *ptr_a, *ptr_x, *ptr_y;
  397.     const int mr = 8, nr = 8;
  398.     if (mr != lb)
  399.     {
  400.         cout << "lb != mr\n";
  401.         system("pause");
  402.         exit(1);
  403.     }
  404.  
  405.     ptr_a = a;
  406.  
  407.     for (i = 0; i < n; i += mr)
  408.     {
  409.         ry0 = _mm256_setzero_pd();
  410.         ry1 = ry2 = ry3 = ry0;
  411.  
  412.         ptr_y = &y[i];
  413.         ptr_x = x;
  414.  
  415.         for (j = 0; j < n; j += nr)
  416.         {
  417.             _mm_prefetch((const char *)(ptr_a + mr * nr), _MM_HINT_NTA);
  418.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 1)), _MM_HINT_NTA);
  419.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 2)), _MM_HINT_NTA);
  420.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 3)), _MM_HINT_NTA);
  421.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 4)), _MM_HINT_NTA);
  422.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 5)), _MM_HINT_NTA);
  423.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 6)), _MM_HINT_NTA);
  424.             _mm_prefetch((const char *)(ptr_a + mr * (nr + 7)), _MM_HINT_NTA);
  425.  
  426.             _mm_prefetch((const char *)(ptr_x + nr), _MM_HINT_T0);
  427.  
  428.             //---------------------------------------------------0
  429.             rx0 = _mm256_broadcast_sd(ptr_x);
  430.             rx1 = _mm256_broadcast_sd(ptr_x + 1);
  431.             ra0 = _mm256_load_pd(ptr_a);
  432.             ra1 = _mm256_load_pd(ptr_a + 4);
  433.             ra2 = _mm256_load_pd(ptr_a + mr);
  434.             ra3 = _mm256_load_pd(ptr_a + mr + 4);
  435.  
  436.             ry0 = _mm256_fmadd_pd(ra0, rx0, ry0);
  437.             ry1 = _mm256_fmadd_pd(ra1, rx0, ry1);
  438.             ry2 = _mm256_fmadd_pd(ra2, rx1, ry2);
  439.             ry3 = _mm256_fmadd_pd(ra3, rx1, ry3);
  440.  
  441.             //---------------------------------------------------1
  442.             rx0 = _mm256_broadcast_sd(ptr_x + 2);
  443.             rx1 = _mm256_broadcast_sd(ptr_x + 3);
  444.             ra0 = _mm256_load_pd(ptr_a + 2 * mr);
  445.             ra1 = _mm256_load_pd(ptr_a + 2 * mr + 4);
  446.             ra2 = _mm256_load_pd(ptr_a + 3 * mr);
  447.             ra3 = _mm256_load_pd(ptr_a + 3 * mr + 4);
  448.  
  449.             ry0 = _mm256_fmadd_pd(ra0, rx0, ry0);
  450.             ry1 = _mm256_fmadd_pd(ra1, rx0, ry1);
  451.             ry2 = _mm256_fmadd_pd(ra2, rx1, ry2);
  452.             ry3 = _mm256_fmadd_pd(ra3, rx1, ry3);
  453.  
  454.             //---------------------------------------------------2
  455.             rx0 = _mm256_broadcast_sd(ptr_x + 4);
  456.             rx1 = _mm256_broadcast_sd(ptr_x + 5);
  457.             ra0 = _mm256_load_pd(ptr_a + 4 * mr);
  458.             ra1 = _mm256_load_pd(ptr_a + 4 * mr + 4);
  459.             ra2 = _mm256_load_pd(ptr_a + 5 * mr);
  460.             ra3 = _mm256_load_pd(ptr_a + 5 * mr + 4);
  461.  
  462.             ry0 = _mm256_fmadd_pd(ra0, rx0, ry0);
  463.             ry1 = _mm256_fmadd_pd(ra1, rx0, ry1);
  464.             ry2 = _mm256_fmadd_pd(ra2, rx1, ry2);
  465.             ry3 = _mm256_fmadd_pd(ra3, rx1, ry3);
  466.  
  467.             //---------------------------------------------------2
  468.             rx0 = _mm256_broadcast_sd(ptr_x + 6);
  469.             rx1 = _mm256_broadcast_sd(ptr_x + 7);
  470.             ra0 = _mm256_load_pd(ptr_a + 6 * mr);
  471.             ra1 = _mm256_load_pd(ptr_a + 6 * mr + 4);
  472.             ra2 = _mm256_load_pd(ptr_a + 7 * mr);
  473.             ra3 = _mm256_load_pd(ptr_a + 7 * mr + 4);
  474.  
  475.             ptr_a += mr * nr;
  476.             ptr_x += nr;
  477.  
  478.             ry0 = _mm256_fmadd_pd(ra0, rx0, ry0);
  479.             ry1 = _mm256_fmadd_pd(ra1, rx0, ry1);
  480.             ry2 = _mm256_fmadd_pd(ra2, rx1, ry2);
  481.             ry3 = _mm256_fmadd_pd(ra3, rx1, ry3);
  482.         }
  483.         ry0 = _mm256_add_pd(ry0, ry2);
  484.         ry1 = _mm256_add_pd(ry1, ry3);
  485.  
  486.         _mm256_store_pd(ptr_y, ry0);
  487.         _mm256_store_pd(ptr_y + 4, ry1);
  488.     }
  489. }
  490. #endif
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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!