#cuda #fortran
Вопрос:
Я выполняю базовое умножение матрицы с использованием CUDA Fortran и C без каких-либо оптимизаций. И Fortran, и C делают одно и то же, но время выполнения для Fortran медленнее.
Ядро C
#define idx(x,y,z) x*y z
__global__ void matrixMultiply(double *d_mat, double *d_matT, double *d_matSym) {
//Global inidices
int tx = blockIdx.y * blockDim.y threadIdx.y;
int ty = blockIdx.x * blockDim.x threadIdx.x;
int k;
if (tx < SIZE amp;amp; ty < SIZE) {
double accum = 0.0;
//Accumulation for (tx,ty) position
for (k=0; k<SIZE; k) {
accum = d_mat[idx(tx,SIZE,k)] * d_matT[idx(k,SIZE,ty)];
}
d_matSym[idx(tx,SIZE,ty)] = accum;
}
}
//Call
dim3 grid_dim(SIZE/32, SIZE/32, 1);
dim3 blk_dim(32, 32, 1);
matrixMultiply<<<grid_dim, blk_dim>>>(d_mat, d_matT, d_matSym);
cudaDeviceSynchronize();
Ядро Фортрана
attributes(global) subroutine matrixMultiply(d_mat, d_matT, d_matSym)
integer :: tx, ty, k
real*8 :: accum
real*8, dimension(:,:) :: d_mat, d_matT, d_matSym
tx = threadIdx%x (blockIdx%x - 1) * blockDim%x
ty = threadIdx%y (blockIdx%y - 1) * blockDim%y
if (tx <= SIZE_ .and. ty <=SIZE_) then
accum = 0.0
do k=1, SIZE_
accum = accum d_mat(tx,k) * d_matT(k,ty)
end do
d_matSym(tx,ty) = accum
end if
end subroutine matrixMultiply
!Call
type(dim3) :: grid_dim, blk_dim
grid_dim = dim3(SIZE_/32, SIZE_/32, 1)
blk_dim = dim3(32, 32, 1)
call matrixMultiply<<<grid_dim, blk_dim>>>(d_mat, d_matT, d_matSym)
err = cudaDeviceSynchronize()
Разница в том, что C использует массив 1D, тогда как Fortran использует 2D. Но это не должно быть проблемой, так как под памятью будет непрерывная память.
Если это доступ к памяти, то в обоих случаях цикл K обращается к одной матрице последовательно, а другой доступ увеличивается по РАЗМЕРУ.
И то, и другое дает одинаковые результаты.
Для матрицы 16384 x 16384, C : 5,4 сек Фортран : 6,3 сек
Графический процессор: Tesla V100 32 ГБ
Я не уверен, что я здесь делаю не так?
Комментарии:
1. Почему вы считаете, что делаете что-то не так? Почему код, написанный на разных языках и скомпилированный с помощью разных компиляторов, должен обеспечивать одинаковую производительность? Разберите получившуюся наглость и поищите отличия. Возможно, ваши неописанные настройки компилятора в обоих случаях отличаются….
2. @talonmies Я согласен с вами, что здесь, возможно, нет ничего плохого. Но этот вопрос справедлив для того факта, что Fortran, как правило, во многих случаях работает так же быстро, как или быстрее, чем C, хотя разница часто минимальна, как и здесь. @abhishekpurandare1297 : У меня нет опыта работы с Cuda Fortran. Но если вы не получите ответа на свой вопрос, предполагая , что ваш компилятор есть
nvfortran, пожалуйста, задайте свой вопрос также на форуме NVIDIA Fortran (и опубликуйте свой ответ здесь тоже): forums.developer.nvidia.com/c/accelerated-computing/…3. @ abhishekpurandare1297 Вы несколько раз просчитывали примеры и усредняли результаты?
4. «Фортран использует 2D» : У меня есть опыт, что большее количество измерений массива часто замедляет код (возможно, из-за большего количества арифметики для вычисления адресов).
5. Просто из любопытства, прикрепление
contiguousатрибута к этой строке «реальный*8, измерение(:,:) :: d_mat, d_matT, d_matSym» может иметь какой-то эффект? (и я также думаю, что форум nvidia может быть полезен)
Ответ №1:
Прежде всего, я предлагаю, чтобы вопросы о производительности включали полные коды. Как правило, мне нужно уметь что-то запускать, и вы можете сэкономить мне время на наборе текста. Конечно, вы можете оставить все как есть. Конечно, я, наверное, смогу понять, что это такое. Но у меня меньше шансов помочь вам таким образом, и я подозреваю, что я не одинок в этом мнении. Мой совет: Сделайте так, чтобы другим было легко вам помогать. Ниже я привел примеры того, что было бы полезно.
Переходим к вопросу:
Разница в том, что C использует массив 1D, тогда как Fortran использует 2D. Но это не должно быть проблемой, так как под памятью будет непрерывная память.
TL;DR: Ваше утверждение («это не должно быть проблемой») , очевидно, не может быть поддержано. Разница между распределением 1D и распределением 2D имеет значение не только с точки зрения хранения, но и с точки зрения расчета индекса. Если вы чувствительны к длине этого ответа, перейдите к примечанию D внизу этого поста.
Подробные сведения:
Когда у нас есть такая петля, как эта:
for (k=0; k<SIZE; k) {
accum = d_mat[idx(tx,SIZE,k)] * d_matT[idx(k,SIZE,ty)];
}
Если компилятор сможет обнаружить, что вычисления индексации могут быть просто «увеличены» по мере выполнения цикла, это будет довольно просто. ОТО, если индексирование должно быть решено с учетом 2D-вычисления (например row*width column ), на которое ранее ссылались на каждой итерации цикла, это будет сложнее. Важным фактором здесь также является то, может ли массив width использоваться/использоваться компилятором во время компиляции.
Похоже, именно это и происходит, и я предполагаю, что это разница во времени. Мы рассмотрим 2 доказательства этого.
- Анализ SASS
Если мы сравним SASS между реализациями, мы увидим некоторые важные различия. Для этого нам понадобятся полные примеры:
CUDA Фортран:
$ cat t11.cuf
module mmk
USE cudafor
!
! Definition of symbols for real types (RP)
!
IMPLICIT NONE
!
INTEGER, PARAMETER :: SP = SELECTED_REAL_KIND(6, 37) ! REAL32
INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15, 307) ! REAL64
INTEGER, PARAMETER :: SIZE_ = 16384
!
Contains
attributes(global) subroutine matrixMultiply(d_mat, d_matT, d_matSym)
integer :: tx, ty, k
REAL(DP) :: accum
REAL(DP), dimension(:,:) :: d_mat, d_matT, d_matSym
tx = threadIdx%x (blockIdx%x - 1) * blockDim%x
ty = threadIdx%y (blockIdx%y - 1) * blockDim%y
if (tx <= SIZE_ .and. ty <=SIZE_) then
accum = 0.0
do k=1, SIZE_
accum = accum d_mat(tx,k) * d_matT(k,ty)
end do
d_matSym(tx,ty) = accum
end if
end subroutine matrixMultiply
end module mmk
PROGRAM Test
!
! This is the main program for Test
!
USE cudafor
USE mmk
!
IMPLICIT NONE
!
REAL(DP), ALLOCATABLE, DEVICE, DIMENSION(:,:) :: d_mat, d_matT, d_matSym
!
INTEGER :: err, i1, i2
type(dim3) :: grid_dim, blk_dim
!
!
! Allocate storage for the arrays
!
Allocate(d_mat(SIZE_,SIZE_),d_matT(SIZE_,SIZE_),d_matSym(SIZE_,SIZE_))
!
! invoke the kernel
!
!Call
grid_dim = dim3(SIZE_/32, SIZE_/32, 1)
blk_dim = dim3(32, 32, 1)
call matrixMultiply<<<grid_dim, blk_dim>>>(d_mat, d_matT, d_matSym)
call matrixMultiply<<<grid_dim, blk_dim>>>(d_mat, d_matT, d_matSym)
err = cudaDeviceSynchronize()
!
! Free storage for the arrays
!
Deallocate(d_mat,d_matT,d_matSym)
!
END PROGRAM Test
$ nvfortran t11.cuf -o t11
$ nvprof ./t11
==38508== NVPROF is profiling process 38508, command: ./t11
==38508== Profiling application: ./t11
==38508== Profiling result:
Type Time(%) Time Calls Avg Min Max Name
GPU activities: 100.00% 12.7168s 2 6.35838s 6.35229s 6.36447s mmk_matrixmultiply_
0.00% 9.2480us 6 1.5410us 1.3120us 2.1760us [CUDA memcpy HtoD]
API calls: 97.08% 12.7260s 9 1.41400s 10.596us 6.36561s cudaFree
2.85% 374.24ms 9 41.583ms 3.5780us 369.55ms cudaMalloc
0.03% 4.5190ms 4 1.1297ms 1.0551ms 1.2608ms cuDeviceTotalMem
0.02% 2.5060ms 404 6.2020us 111ns 1.2057ms cuDeviceGetAttribute
0.01% 1.1501ms 4 287.52us 33.014us 1.0465ms cuDeviceGetName
0.00% 120.00us 6 20.000us 5.1530us 62.523us cudaMemcpy
0.00% 99.593us 2 49.796us 44.275us 55.318us cudaLaunchKernel
0.00% 43.414us 1 43.414us 43.414us 43.414us cudaDeviceSynchronize
0.00% 11.563us 4 2.8900us 1.4400us 5.5870us cuDeviceGetPCIBusId
0.00% 1.7030us 8 212ns 124ns 580ns cuDeviceGet
0.00% 901ns 3 300ns 163ns 534ns cuDeviceGetCount
0.00% 855ns 4 213ns 183ns 263ns cuDeviceGetUuid
$
NVIDIA HPC SDK, 2021.2, Tesla V100
Мы видим, что продолжительность выполнения составляет около 6,3 секунды, в соответствии с вашим отчетом.
CUDA C :
$ cat t12.cu
#define idx(x,y,z) x*y z
const int SIZE = 16384;
__global__ void matrixMultiply(double *d_mat, double *d_matT, double *d_matSym) {
//Global inidices
int tx = blockIdx.y * blockDim.y threadIdx.y;
int ty = blockIdx.x * blockDim.x threadIdx.x;
int k;
if (tx < SIZE amp;amp; ty < SIZE) {
double accum = 0.0;
//Accumulation for (tx,ty) position
for (k=0; k<SIZE; k) {
accum = d_mat[idx(tx,SIZE,k)] * d_matT[idx(k,SIZE,ty)];
}
d_matSym[idx(tx,SIZE,ty)] = accum;
}
}
int main(){
//Call
dim3 grid_dim(SIZE/32, SIZE/32, 1);
dim3 blk_dim(32, 32, 1);
double *d_mat, *d_matT, *d_matSym;
cudaMalloc(amp;d_mat, SIZE*SIZE*sizeof(double));
cudaMalloc(amp;d_matT, SIZE*SIZE*sizeof(double));
cudaMalloc(amp;d_matSym, SIZE*SIZE*sizeof(double));
matrixMultiply<<<grid_dim, blk_dim>>>(d_mat, d_matT, d_matSym);
matrixMultiply<<<grid_dim, blk_dim>>>(d_mat, d_matT, d_matSym);
cudaDeviceSynchronize();
}
$ nvcc -o t12 t12.cu -arch=sm_70
$ nvprof ./t12
==40364== NVPROF is profiling process 40364, command: ./t12
==40364== Profiling application: ./t12
==40364== Profiling result:
Type Time(%) Time Calls Avg Min Max Name
GPU activities: 100.00% 11.0379s 2 5.51893s 5.50503s 5.53282s matrixMultiply(double*, double*, double*)
API calls: 97.80% 11.0379s 1 11.0379s 11.0379s 11.0379s cudaDeviceSynchronize
2.15% 242.74ms 3 80.914ms 1.9747ms 238.78ms cudaMalloc
0.04% 4.3105ms 4 1.0776ms 1.0728ms 1.0890ms cuDeviceTotalMem
0.01% 1.4584ms 404 3.6090us 110ns 163.21us cuDeviceGetAttribute
0.00% 142.28us 4 35.571us 32.475us 41.980us cuDeviceGetName
0.00% 51.695us 2 25.847us 7.0790us 44.616us cudaLaunchKernel
0.00% 11.198us 4 2.7990us 1.6260us 5.5950us cuDeviceGetPCIBusId
0.00% 1.7930us 3 597ns 149ns 1.2200us cuDeviceGetCount
0.00% 1.6590us 8 207ns 120ns 704ns cuDeviceGet
0.00% 733ns 4 183ns 161ns 215ns cuDeviceGetUuid
$
CUDA 11.2, та же машина V100
Снова мы видим время работы ядра ~5,5 секунд, совпадающее с вашим отчетом. Теперь давайте сравним НАХАЛЬСТВО. Мы используем cubojdump утилиту для его извлечения, и в обоих случаях очевидно, что компилятор каким-то образом разворачивает цикл ядра. Ключевая операция DFMA -сложение с двойной точностью с плавким умножением , которое вызывается один раз за итерацию цикла, для вычисления одного продукта (и сохранения текущей суммы). Поэтому я буду извлекать только короткую последовательность операций DFMA в каждом случае, в середине развернутого цикла, чтобы получить представление о «стоимости» и «накладных расходах»на каждом флопе.
CUDA Фортран:
/*1670*/ LDG.E.64.SYS R18, [R12 0x1b8] ; /* 0x0001b8000c127381 */
/* 0x000ea200001eeb00 */
/*1680*/ IADD3 R20, P0, R16, R29, RZ ; /* 0x0000001d10147210 */
/* 0x001fc60007f1e0ff */
/*1690*/ LDG.E.64.SYS R22, [R12 0x1b0] ; /* 0x0001b0000c167381 */
/* 0x000e6400001eeb00 */
/*16a0*/ IMAD.X R21, R17, 0x1, R27, P0 ; /* 0x0000000111157824 */
/* 0x000fe400000e061b */
/*16b0*/ LDG.E.64.SYS R24, [R16] ; /* 0x0000000010187381 */
/* 0x00006c00001eeb00 */
/*16c0*/ LDG.E.64.SYS R16, [R20] ; /* 0x0000000014107381 */
/* 0x001ea200001eeb00 */
/*16d0*/ DFMA R22, R24, R22, R14 ; /* 0x000000161816722b */
/* 0x002084000000000e */
/*16e0*/ IADD3 R14, P0, R20, R29, RZ ; /* 0x0000001d140e7210 */
/* 0x001fca0007f1e0ff */
/*16f0*/ IMAD.X R15, R21, 0x1, R27, P0 ; /* 0x00000001150f7824 */
/* 0x000fe200000e061b */
/*1700*/ DFMA R16, R16, R18, R22 ; /* 0x000000121010722b */
Описанная выше последовательность повторяется в области развернутого цикла. Мы отмечаем, что существует 2 DFMA операции, 4 LDG операции (имеет смысл — по две на умножение), а остальные инструкции (4) предположительно представляют собой обновления индексации по мере прохождения цикла через его итерации.
CUDA C :
/*02d0*/ LDG.E.64.SYS R16, [R6 0x40000] ; /* 0x0400000006107381 */
/* 0x001f2800001eeb00 */
/*02e0*/ LDG.E.64.SYS R24, [R4 0x10] ; /* 0x0000100004187381 */
/* 0x000f2200001eeb00 */
/*02f0*/ DFMA R10, R20, R10, R22 ; /* 0x0000000a140a722b */
/* 0x0200860000000016 */
/*0300*/ LDG.E.64.SYS R22, [R6 0x60000] ; /* 0x0600000006167381 */
/* 0x001f6800001eeb00 */
/*0310*/ LDG.E.64.SYS R20, [R4 0x18] ; /* 0x0000180004147381 */
/* 0x000f6200001eeb00 */
/*0320*/ DFMA R14, R8, R14, R10 ; /* 0x0000000e080e722b */
И снова эта последовательность повторяется. Мы видим 2 DFMA и 4 LDG инструкции, но никаких других инструкций в повторяющейся последовательности. Отметим, что LDG инструкции, похоже, взаимозачетов, которые могут быть включены во время компиляции в инструкцию, и эти смещения исключить необходимость каких-либо дополнительных инструкций, ведь они только «приращение» ранее вычисленных смещений (Отметим также, что «с шагом» по-видимому, на уровне столбца смещения в одном случае, и построчный смещения в друга, как можно было бы ожидать, чтобы забрать умножение операндов в ядре петли.)
- Рефакторинг
Можем ли мы реорганизовать код Fortran, чтобы он работал так же, как код C ? ДА:
$ cat t11a.cuf
module mmk
USE cudafor
!
! Definition of symbols for real types (RP)
!
IMPLICIT NONE
!
INTEGER, PARAMETER :: SP = SELECTED_REAL_KIND(6, 37) ! REAL32
INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15, 307) ! REAL64
INTEGER, PARAMETER :: SIZE_ = 16384
!
Contains
attributes(global) subroutine matrixMultiply(d_mat, d_matT, d_matSym)
integer :: tx, ty, k
REAL(DP) :: accum
REAL(DP), dimension(:) :: d_mat, d_matT, d_matSym
tx = threadIdx%x (blockIdx%x - 1) * blockDim%x
ty = threadIdx%y (blockIdx%y - 1) * blockDim%y
if (tx <= SIZE_ .and. ty <=SIZE_) then
accum = 0.0
do k=1, SIZE_
accum = accum d_mat((ty-1)*SIZE_ k) * d_matT((k-1)*SIZE_ tx)
end do
d_matSym((ty-1)*SIZE_ tx) = accum
end if
end subroutine matrixMultiply
end module mmk
PROGRAM Test
!
! This is the main program for Test
!
USE cudafor
USE mmk
!
IMPLICIT NONE
!
REAL(DP), ALLOCATABLE, DEVICE, DIMENSION(:) :: d_mat, d_matT, d_matSym
!
INTEGER :: err, i1, i2
type(dim3) :: grid_dim, blk_dim
!
! Allocate storage for the arrays
!
Allocate(d_mat(SIZE_*SIZE_),d_matT(SIZE_*SIZE_),d_matSym(SIZE_*SIZE_))
!
! invoke the kernel
!
!Call
grid_dim = dim3(SIZE_/32, SIZE_/32, 1)
blk_dim = dim3(32, 32, 1)
call matrixMultiply<<<grid_dim, blk_dim>>>(d_mat, d_matT, d_matSym)
call matrixMultiply<<<grid_dim, blk_dim>>>(d_mat, d_matT, d_matSym)
err = cudaDeviceSynchronize()
!
! Free storage for the arrays
!
Deallocate(d_mat,d_matT,d_matSym)
!
END PROGRAM Test
$ nvfortran t11a.cuf -o t11a
$ nvprof ./t11a
==45544== NVPROF is profiling process 45544, command: ./t11a
==45544== Profiling application: ./t11a
==45544== Profiling result:
Type Time(%) Time Calls Avg Min Max Name
GPU activities: 100.00% 10.8169s 2 5.40847s 5.39118s 5.42576s mmk_matrixmultiply_
0.00% 9.2160us 6 1.5360us 1.3120us 2.1440us [CUDA memcpy HtoD]
API calls: 96.72% 10.8240s 9 1.20266s 12.102us 5.42594s cudaFree
3.22% 360.34ms 9 40.038ms 3.2760us 355.64ms cudaMalloc
0.04% 4.3488ms 4 1.0872ms 1.0598ms 1.1593ms cuDeviceTotalMem
0.02% 2.3633ms 404 5.8490us 111ns 869.22us cuDeviceGetAttribute
0.00% 144.50us 4 36.125us 33.254us 41.317us cuDeviceGetName
0.00% 106.43us 6 17.737us 4.7910us 51.453us cudaMemcpy
0.00% 101.46us 2 50.732us 44.479us 56.985us cudaLaunchKernel
0.00% 20.926us 1 20.926us 20.926us 20.926us cudaDeviceSynchronize
0.00% 7.2850us 4 1.8210us 491ns 4.6210us cuDeviceGetPCIBusId
0.00% 1.7650us 8 220ns 118ns 672ns cuDeviceGet
0.00% 926ns 4 231ns 218ns 264ns cuDeviceGetUuid
0.00% 614ns 3 204ns 130ns 310ns cuDeviceGetCount
$
И мы видим, что производительность примерно такая же, как у версии CUDA C .
Примечания:
О. Я опустил анализ SASS для окончательного варианта. Упражнение оставлено на усмотрение читателя. Однако развернутая часть цикла представляет собой повторяющуюся последовательность LDG-LDG-DFMA, как это видно в версии CUDA C .
Б. У вас могут возникнуть два вопроса: 1. Почему это происходит именно так? Я пытался ответить на этот вопрос. 2. Должно ли это быть так/должно ли это быть так? Я не ответил на этот вопрос, я полагаю, что это комбинация некоторых утверждений Fortran, которые я не хочу делать, а также вопрос для инженеров-компиляторов. Если вы считаете, что 2D-реализация должна дублировать 1D-реализацию, я предлагаю сообщить об ошибке.
C. Я немного поскреб в затылке, прежде чем заметил, что вы сделали это в случае C , но не в случае fortran:
int tx = blockIdx.y * blockDim.y threadIdx.y;
^ ^ ^ ^
D. Одна из причин, по которой я считаю, что это несоответствие не поддается тривиальному разрешению во время компиляции, заключается в том, что обработка 2D-массива компилятором в CUDA Fortran (независимо от того, нужна она или нет), по-видимому, связана с созданием пакета метаданных во время выполнения, который включает ширину массива, возможно, среди другой информации. Этот пакет метаданных передается на устройство отдельно от полезной нагрузки массива данных-об этом можно догадаться из приведенного выше вывода профилировщика или путем запуска, --print-gpu-trace чтобы увидеть его более четко. В вашем случае CUDA C , по сути, ширина массива ( SIZE ) известна во время компиляции. В случае Фортрана это не так (или, в качестве альтернативы, если это так, компилятор, похоже, не распространяет эти знания через генерацию кода ядра). Вы можете видеть, что «улучшенный» режим адресации выше предполагает смещения, которые известны во время компиляции. Можно ли решить это «должно» или нет при более тщательном обращении с Фортраном, я не могу сказать. Однако я бы ожидал, что какое-то contiguous объявление/украшение само по себе недостаточно для устранения этого несоответствия для компилятора. Одно это объявление не устранит необходимость в вычислениях индексации во время выполнения на лету, основанных на ширине массива.
Комментарии:
1. Когда компилятор поддерживает ОБЪЕДИНЕНИЕ/СОПОСТАВЛЕНИЕ, обычно может быть удобно работать с матрицей как с 2D. Но для операций ЭЛЕМЕНТАРНОГО типа обратитесь к нему как к 1D массиву/вектору.
2. Спасибо, что нашли для этого время. Я действительно рад, что у вас все еще есть энтузиазм по поводу такой глубины анализа. Обычно я этого не делаю. Как обычно, многие из новых функций Fortran (указатели, интеллектуальное индексирование массивов, нарезка, широковещательная передача и арифметика и т. Д.) Имеют удивительно жадные реализации, которые ослабили большинство традиционных преимуществ производительности устаревшего Fortran по сравнению с «новыми» языками….
3. Проблема заключалась в том, что компилятор не знал размер массива. Я все еще новичок в Fortran, поэтому я ошибочно проигнорировал определение размера массивов параметров ядра.
real*8, dimension(SIZE_,SIZE_) :: d_mat, d_matT, d_matSymдекларация устранила проблему. Я ценю ваше время и глубокое понимание!