2014-10-26 2 views
1

Я пытаюсь векторизовать умножение матрицы, используя блокировку и векторные свойства. Мне кажется, что добавочная часть в векторном умножении не может быть векторизована. Не могли бы вы посмотреть, могу ли я улучшить свой код, чтобы векторизовать дальше?Векторизация добавочной части матричного умножения с использованием встроенных функций?

double dd[4], bb[4]; 
    __m256d op_a, op_b, op_d; 
    for(i = 0; i < num_blocks; i++){ 
     for(j = 0; j < num_blocks; j++){ 
      for(k = 0; k < num_blocks; k++){ 
       for(ii = 0; ii < block_size ; ii++){ 
        for(kk = 0; kk < block_size; kk++){ 
         for(jj = 0; jj < block_size ; jj+=4){ 

          aoffset=n*(i*block_size+ii)+j*block_size +jj ; 
          boffset=n*(j*block_size+jj)+k*block_size +kk; 
          coffset=n*(i*block_size+ii)+ k*block_size + kk; 

          bb[0]=b[n*(j*block_size+jj)+k*block_size +kk]; 
          bb[1]=b[n*(j*block_size+jj+1)+k*block_size +kk]; 
          bb[2]=b[n*(j*block_size+jj+2)+k*block_size +kk]; 
          bb[3]=b[n*(j*block_size+jj+3)+k*block_size +kk]; 

          op_a = _mm256_loadu_pd (a+aoffset); 
          op_b= _mm256_loadu_pd (bb); 
          op_d = _mm256_mul_pd(op_a, op_b); 
          _mm256_storeu_pd (dd, op_d); 
          c[coffset]+=(dd[0]+dd[1]+dd[2]+dd[3]); 

         } 
        } 
       } 
      } 
     } 
    } 

Спасибо.

+1

У вашего ядра несколько неэффективности. Для этого вы должны переупорядочить внутренний блок, чтобы вы могли читать его горизонтально, а не вертикально, на 'b'. Вы также должны хранить непосредственно в 'c', а не использовать' dd'. Но есть больше проблем. Я постараюсь ответить в течение следующих дней, если кто-то еще этого не сделает. Вы должны измерить эффективность вашего метода, чтобы у вас была ссылка. –

+1

Не можете ли вы использовать существующую реализацию? Матричное умножение - это хорошо протекающий дерн, и есть несколько существующих библиотек, которые ускоряют его. Что вы пытаетесь достичь за пределами предшествующего уровня техники? –

+1

@JohnZwinck Я пытаюсь научиться делать Matix для умножения матриц с помощью векторизации и векторных свойств. Затем один раз, я управляю своим пониманием по векторизации в одиночном, я попытаюсь узнать о MPI и OpenMP. Это не профессиональная задача, это мое самообучение. –

ответ

1

Вы можете использовать эту версию матричного умножения (с [I, J] = а [I, K] * б [K, J]) алгоритм (скалярную версия):

for(int i = 0; i < i_size; ++i) 
{ 
    for(int j = 0; j < j_size; ++j) 
     c[i][j] = 0; 

    for(int k = 0; k < k_size; ++k) 
    { 
     double aa = a[i][k]; 
     for(int j = 0; j < j_size; ++j) 
      c[i][j] += aa*b[k][j]; 
    } 
} 

И векторизованную версия :

for(int i = 0; i < i_size; ++i) 
{ 
    for(int j = 0; j < j_size; j += 4) 
     _mm256_store_pd(c[i] + j, _mm256_setzero_pd()); 

    for(int k = 0; k < k_size; ++k) 
    { 
     __m256d aa = _mm256_set1_pd(a[i][k]); 
     for(int j = 0; j < j_size; j += 4) 
     { 
      _mm256_store_pd(c[i] + j, _mm256_add_pd(_mm256_load_pd(c[i] + j), _mm256_mul_pd(aa, _mm256_load_pd(b[k] + j)))); 
     } 
    } 
} 
+1

У этого есть ужасные характеристики кеширования по сравнению с кодом в вопросе, и он также не считывает правильные значения из 'b' (или' a'). Один - это вектор столбца, один - вектор строки, и вы относитесь к ним одинаково. –

0

«Горизонтальные добавить» является более недавним изданием для набора инструкций SSE, так что вы не можете использовать ускоренную версию, если совместимость со многими различными процессорами вашей цели.

Однако, вы определенно можете векторизовать дополнения. Обратите внимание, что внутренний цикл влияет только на один coffset. Вы должны переместить вычисление coffset наружу (компилятор сделает это автоматически, но код будет более читабельным, если вы это сделаете), а также используйте четыре аккумулятора в самом внутреннем цикле, выполняя горизонтальное добавление только один раз за coffset. Это улучшение, даже если используется векторное горизонтальное добавление, а для скалярного горизонтального добавления оно довольно большое.

Что-то вроде:

for(kk = 0; kk < block_size; kk++){ 
    op_e = _mm256_setzero_pd(); 

    for(jj = 0; jj < block_size ; jj+=4){ 
     aoffset=n*(i*block_size+ii)+j*block_size +jj ; 
     boffset=n*(j*block_size+jj)+k*block_size +kk; 

     bb[0]=b[n*(j*block_size+jj)+k*block_size +kk]; 
     bb[1]=b[n*(j*block_size+jj+1)+k*block_size +kk]; 
     bb[2]=b[n*(j*block_size+jj+2)+k*block_size +kk]; 
     bb[3]=b[n*(j*block_size+jj+3)+k*block_size +kk]; 

     op_a = _mm256_loadu_pd (a+aoffset); 
     op_b= _mm256_loadu_pd (bb); 
     op_d = _mm256_mul_pd(op_a, op_b); 
     op_e = _mm256_add_pd(op_e, op_d); 
    } 
    _mm256_storeu_pd(dd, op_e); 
    coffset = n*(i*block_size+ii)+ k*block_size + kk; 
    c[coffset] = (dd[0]+dd[1]+dd[2]+dd[3]); 
} 

Вы также можете ускорить этот процесс, делая транспонирование на b заранее, вместо того, чтобы собрать вектор внутри самого внутреннего цикла.