Точечный продукт с использованием Cblas медленно - PullRequest
Новая
       7

Точечный продукт с использованием Cblas медленно

2 голосов
/ 27 февраля 2012

Я хочу рассчитать произведение A ^ T * A (A - матрица 2000x1000).Также я хочу решить только верхнюю треугольную матрицу.Во внутреннем цикле я должен решить скалярное произведение двух векторов.

Теперь вот проблема.Использование cblas ddot () не быстрее, чем вычисление скалярного произведения с помощью цикла.Как это возможно?(с использованием процессора Intel Core i7 M720 @ 2,67 ГГц, 1,92 ГБ ОЗУ)

Ответы [ 3 ]

4 голосов
/ 30 ноября 2013

Проблема вызвана в основном размером матрицы, а не ddot. Ваши матрицы настолько велики, что не помещаются в кеш-память. Решение состоит в том, чтобы переставить три вложенных цикла таким образом, чтобы как можно больше можно было выполнить с помощью строки в кеше, что сокращает обновления кеша. Реализация модели следует как для подхода ddot, так и для подхода daxpy. На моем компьютере потребление времени было около 15: 1. Другими словами: никогда, никогда, никогда не программируйте умножение матриц по схеме «столбец времен строки», которую мы изучали в школе. 

    /*
      Matrix product of A^T * A by two methods.
      1) "Row times column" as we learned in school.
      2) With rearranged loops such that need for cash refreshes is reduced
         (this can be improved even more).

      Compile:  gcc -o aT_a aT_a.c -lgslcblas -lblas -lm
    */

    #include <stdio.h>
    #include <stdlib.h>
    #include <time.h>
    #include <cblas.h>

    #define ROWS 2000
    #define COLS 1000

    static double a[ROWS][COLS];
    static double c[COLS][COLS];

    static void dot() {
      int i, j;
      double *ai, *bj;
      ai = a[0];
      for (i=0; i<COLS; i++) {
        bj = a[0];
        for (j=0; j<COLS; j++) {
          c[i][j] = cblas_ddot(ROWS,ai,COLS,bj,COLS);
          bj += 1;
        }
        ai += 1;
      }
    }

    static void axpy() {
      int i, j;
      double *ci, *bj, aij;
      for (i=0; i<COLS; i++) {
        ci = c[i];
        for (j=0; j<COLS; j++) ci[j] = 0.;
        for (j=0; j<ROWS; j++) {
          aij = a[j][i];
          bj = a[j];
          cblas_daxpy(COLS,aij,bj,1,ci,1);
        }
      }
    }

    int main(int argc, char** argv) {
      clock_t t0, t1;
      int i, j;

      for (i=0; i<ROWS; ++i) 
        for (j=0; j<COLS; ++j)
          a[i][j] = i+j;

      t0 = clock();
      dot();
      t0 = clock();
      printf("Time for DOT : %f sec.\n",(double)t0/CLOCKS_PER_SEC);
      axpy();
      t1 = clock();
      printf("Time for AXPY: %f sec.\n",(double)(t1-t0)/CLOCKS_PER_SEC);

      return 0;
    }
1 голос
/ 27 февраля 2012

Точечный продукт CBLAS - это просто вычисление в слегка развернутом цикле. Netlib Fortran - это просто: 

     DO I = MP1,N,5
      DTEMP = DTEMP + DX(I)*DY(I) + DX(I+1)*DY(I+1) +
 $            DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4)
     END DO

есть. просто петля развернута до шага 5.

Если вам необходимо использовать точечное изделие в стиле ddot, вы можете получить повышение производительности, переписав цикл для использования встроенных функций SSE2: 

#include <emmintrin.h>

double ddotsse2(const double *x, const double *y, const int n)
{
    double result[2];
    int n2 = 2 * (n/2);
    __m128d dtemp;

    if ( (n % 2) == 0) {
        dtemp = _mm_setzero_pd(); 
    }  else {
        dtemp = _mm_set_sd(x[n] * y[n]);
    }

    for(int i=0; i<n2; i+=2) {
        __m128d x1 = _mm_loadr_pd(x+i);
        __m128d y1 = _mm_loadr_pd(y+i);
        __m128d xy = _mm_mul_pd(x1, y1);
        dtemp = _mm_add_pd(dtemp, xy);
    }

    _mm_store_pd(&result[0],dtemp);

    return result[0] + result[1];
}

(не проверено, никогда не компилировалось, покупатель остерегается).

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

0 голосов
/ 18 июня 2018

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

void matMulDotProduct(int n, float *A, float* B, int a_size, int b_size, int a_row, int a_col, int b_row, int b_col, float *C) {
    int i, j, k;
    MKL_INT incx, incy;

    incx = 1;
    incy = b_size;

    //copy out multiplying matrix from larger matrix
    float *temp = (float*) malloc(n * n * sizeof(float));
    for (i = 0; i < n; ++i) {
        cblas_scopy(n, &B[(b_row * b_size) + b_col + i], incy, &temp[i * n], 1);
    }

    //transpose
    mkl_simatcopy('R', 'T', n, n, 1.0, temp, 1, 1);

    for (i = 0; i < n; i+= BLOCK_SIZE) {
        for (j = 0; j < n; j++) {
            for (k = 0; k < BLOCK_SIZE; ++k) {
                C[((i + k) * n) + j] = cblas_sdot(n, &A[(a_row + i + k) * a_size + a_col], incx, &temp[n * j], 1);
            }
        }
    }
    free(temp);
}

На моей машине этот код примерно на 1 порядок быстрее, чем код с 3 циклами (но также на 1 порядок медленнее, чем вызов cblas_? Gemm) для операций с плавающей запятой одинарной точности и 2K на 2Kматрицы.(Я использую Intel MKL).

...