Оптимизирует ли большинство компиляторов MATMUL (TRANSPOSE (A), B)?

Question

В программе на Фортране мне нужно вычислить несколько выражений, таких как M · v , M ^T · v , M ^T · M , M · M ^T и т. Д. Здесь M и v - это двумерные и одномерные массивы небольшого размера (меньше чем100, обычно около 2-10)

Мне было интересно, развернется ли запись MATMUL(TRANSPOSE(M),v) во время компиляции в некоторый код, эффективный как MATMUL(N,v), где N явно хранится как N=TRANSPOSE(M).Меня особенно интересуют компиляторы gnu и ifort с «сильными» флагами оптимизации (например, -O2, -O3 или -Ofast).

Answer 1

Ниже вы найдете несколько времен выполнения различных методов.

system:

• Процессор Intel® Core ™ TM i5-6500T @ 2,50 ГГц
• Размер кэша: 6144 КБ
• ОЗУ: 16 МБ
• GNU Fortran (GCC) 6.3.1 20170216 (Red Hat 6.3.1-3)
• ifort (IFORT) 18.0.5 20180823
• BLAS: для компилятора gnu по умолчанию используется blasверсия

сборник:

[gnu] $ gfortran -O3 x.f90 -lblas
[intel] $ ifort -O3 -mkl x.f90

исполнение:

[gnu] $ ./a.out > matmul.gnu.txt
[intel] $ EXPORT MKL_NUM_THREADS=1; ./a.out > matmul.intel.txt

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

вектор матрицы времени

Сравнивалось шесть различных реализаций:

1. manual: do j=1,n; do k=1,n; w(j) = P(j,k)*v(k); end do; end do
2. matmul: matmul(P,v)
3. blas N: dgemv('N',n,n,1.0D0,P,n,v,1,0,w,1)
4. matmul-transpose: matmul(transpose(P),v)
5. matmul-transpose-tmp: Q=transpose(P); w=matmul(Q,v)
6. blas T: dgemv('T',n,n,1.0D0,P,n,v,1,0,w,1)

На рисунке 1и на рисунке 2 вы можете сравнить результаты хронометража для вышеупомянутых случаев.В целом мы можем сказать, что использование временного в gfortran и ifort не рекомендуется.Оба компилятора могут оптимизировать MATMUL(TRANSPOSE(P),v) намного лучше.В то время как в gfortran реализация MATMUL быстрее, чем по умолчанию BLAS, ifort ясно показывает, что mkl-blas быстрее.

рисунок 1: Матрично-векторное умножение.Сравнение различных реализаций выполнялось на gfortran.Левые панели показывают абсолютную синхронизацию, деленную на общее время ручного вычисления для системы размером 1000. Правые панели показывают абсолютную синхронизацию, деленную на n2 × δ.Здесь δ - среднее время ручного вычисления размера 1000, деленное на 1000 × 1000.

рисунок 2: Матрицаумножение вектораСравнение различных реализаций выполнялось на однопоточном компиляторе ifort.Левые панели показывают абсолютную синхронизацию, деленную на общее время ручного вычисления для системы размером 1000. Правые панели показывают абсолютную синхронизацию, деленную на n2 × δ.Здесь δ - среднее время ручного вычисления размера 1000, деленное на 1000 × 1000.

матрица-матрица

Сравнивалось шесть различных реализаций:

1. руководство: do l=1,n; do j=1,n; do k=1,n; Q(j,l) = P(j,k)*P(k,l); end do; end do; end do
2. matmul: matmul(P,P)
3. blas N: dgemm('N','N',n,n,n,1.0D0,P,n,P,n,0.0D0,R,n)
4. matmul-transpose: matmul(transpose(P),P)
5. matmul-transpose-tmp: Q=transpose(P); matmul(Q,P)
6. blas T: dgemm('T','N',n,n,n,1.0D0,P,n,P,n,0.0D0,R,n)

На рисунках 3 и 4 вы можете сравнить результаты синхронизации для вышеупомянутых случаев.В отличие от вектора-случая, использование временного рекомендуется только для gfortran.В то время как в gfortran реализация MATMUL быстрее, чем BLAS по умолчанию, ifort ясно показывает, что mkl-blas быстрее.Примечательно, что ручная реализация сравнима с mkl-blas.

рисунок 3: Матрица-матричное умножение.Сравнение различных реализаций выполнялось на gfortran.Левые панели показывают абсолютную синхронизацию, деленную на общее время ручного вычисления для системы размером 1000. Правые панели показывают абсолютную синхронизацию, деленную на n3 × δ.Здесь δ - среднее время ручного вычисления размера 1000, деленное на 1000 × 1000 × 1000.

рисунок 4: Матрица-матричное умножение.Сравнение различных реализаций выполнялось на однопоточном компиляторе ifort.Левые панели показывают абсолютную синхронизацию, деленную на общее время ручного вычисления для системы размером 1000. Правые панели показывают абсолютную синхронизацию, деленную на n3 × δ.Здесь δ - среднее время ручного вычисления размера 1000, деленное на 1000 × 1000 × 1000.

---

Используемый код:

program matmul_test

  implicit none

  double precision, dimension(:,:), allocatable :: P,Q,R
  double precision, dimension(:), allocatable :: v,w

  integer :: n,i,j,k,l
  double precision,dimension(12) :: t1,t2

  do n = 1,1000
     allocate(P(n,n),Q(n,n), R(n,n), v(n),w(n))
     call random_number(P)
     call random_number(v)

     i=0

     i=i+1
     call cpu_time(t1(i))
     do j=1,n; do k=1,n; w(j) = P(j,k)*v(k); end do; end do
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     w=matmul(P,v)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     call dgemv('N',n,n,1.0D0,P,n,v,1,0,w,1)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     w=matmul(transpose(P),v)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     Q=transpose(P)
     w=matmul(Q,v)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     call dgemv('T',n,n,1.0D0,P,n,v,1,0,w,1)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     do l=1,n; do j=1,n; do k=1,n; Q(j,l) = P(j,k)*P(k,l); end do; end do; end do
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     Q=matmul(P,P)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     call dgemm('N','N',n,n,n,1.0D0,P,n,P,n,0.0D0,R,n)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     Q=matmul(transpose(P),P)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     Q=transpose(P)
     R=matmul(Q,P)
     call cpu_time(t2(i))

     i=i+1
     call cpu_time(t1(i))
     call dgemm('T','N',n,n,n,1.0D0,P,n,P,n,0.0D0,R,n)
     call cpu_time(t2(i))

     write(*,'(I6,12D25.17)') n, t2-t1
     deallocate(P,Q,R,v,w)
  end do

end program matmul_test