Быстрый цикл умножения матриц для C = AA'B - PullRequest
Новая
       76

Быстрый цикл умножения матриц для C = AA'B

0 голосов
/ 26 февраля 2019

Задача состоит в том, чтобы сделать операцию C = AA'B , где A, B и C - прямоугольные массивы.Это может быть достигнуто через BLAS с помощью последовательности вызовов gemm с использованием промежуточной матрицы T , например. T = A'B и C = AT '.Мне интересно, существует ли «более умный» способ сделать это, поскольку AA'B можно переформулировать как последовательность продуктов внешнего вектора.Таким образом, простой цикл в Фортране выглядит следующим образом: 

    allocate(xxx(size(A,1)),source=0.0)
    Do i=1,size(B,2)
      !$OMP PARALLEL PRIVATE(tmp) FIRSTPRIVATE(xxx) REDUCTION(+:C)
      xxx=0.0
      !$OMP DO
      Do j=1,size(A,2)
        tmp=0.0D0;
        Do k=1,size(A,1)
          tmp=tmp+A(k,j)*B(k,i)
        End Do
        if(j==1) Then
          Do k=1,size(A,1)
            xxx(k)=tmp*A(k,j)
          end Do
        Else
          Do k=1,size(A,1)
            xxx(k)=xxx(k)+tmp*A(k,j)
          end Do
        End if
      End Do
      !$OMP END DO
      C(:,i)=C(:,i)+xxx
      !$OMP END PARALLEL
    end Do

Это оптимальная схема с точки зрения структуры доступа к памяти?И если да, то как это организовать с помощью OpenMP.
Все мои попытки создания потоков были примерно в 15 раз медленнее, чем две последовательные операции с гемами через MKL.Так что я думаю, что я делаю что-то не так здесь.

...