Первая проблема заключается в том, что ваши метрики синхронизации не точны.
Linalg компилирует модули cuda на лету, и вы можете увидеть, что код компилируется при его запуске . Я сделал несколько небольших изменений в вашем коде, чтобы уменьшить размер умножаемых массивов, но, несмотря на это, после двух запусков без каких-либо других улучшений я увидел огромный прирост производительности, например:
Time taken with GPU:2.5476348400115967
Time taken without GPU:0.16627931594848633
против
Time taken with GPU:0.8741757869720459
Time taken without GPU:0.15836167335510254
Однако это все еще намного медленнее, чем версия процессора. Следующее, что я сделал, это дал более точную синхронизацию, основанную на том, где происходят реальные вычисления. Вы не используете свою версию numpy, так что не рассчитывайте ее в своей версии cuda:
REAL Time taken with GPU:0.6461708545684814
Вы также копируете в графический процессор и включаете это в вычисления, но это само по себе занимает нетривиальное время, поэтому давайте удалим это:
t1 = time.time()
zz_Hypg = linalg.multiply(zzg,xxg)
t2 = time.time()
...
REAL Time taken with GPU:0.3689603805541992
Ух ты, это очень помогло. Но мы все еще медленнее, чем NumPy версия? Зачем?
Помните, когда я сказал, что numpy не укладывается? Он не копирует память вообще для широкого вещания. Чтобы получить реальную скорость, вам нужно:
- не плитка
- вещательные размеры
- реализовать это в ядре.
Pycuda предоставляет утилиты для реализации ядра, но его массив GPU не обеспечивает вещание. По сути, вам нужно сделать следующее (ОТКАЗ ОТ ОТВЕТСТВЕННОСТИ: я не проверял это, возможно, есть ошибки, это просто для демонстрации того, как должно выглядеть ядро):
#include <pycuda-complex.hpp>
//KERNEL CODE
constexpr unsigned work_tile_dim = 32
//instruction level parallelism factor, how much extra work to do per thread, may be changed but effects the launch dimensions. thread group size should be (tile_factor, tile_factor/ilp_factor)
constexpr unsigned ilp_factor = 4
//assuming c order:
// x axis contiguous out,
// y axis contiguous in zz,
// x axis contiguous in xx
// using restrict because we know that all pointers will refer to different parts of memory.
__global__
void element_wise_multiplication(
pycuda::complex<double>* __restrict__ array_zz,
pycuda::complex<double>* __restrict__ array_xx,
pycuda::complex<double>* __restrict__ out_array,
unsigned array_zz_w, /*size of w,z,y, dimensions used in zz*/
unsigned array_zz_z,
unsigned array_zz_xx_y,/*size of y,x, dimensions used in xx, but both have same y*/
unsigned array_xx_x){
// z dimensions in blocks often have restrictions on size that can be fairly small, and sometimes can cause performance issues on older cards, we are going to derive x,y,z,w index from just the x and y indicies instead.
unsigned x_idx = blockIdx.x * (work_tile_dim) + threadIdx.x
unsigned y_idx = blockIdx.y * (work_tile_dim) + threadIdx.y
//blockIdx.z stores both z and w and should not over shoot, and aren't used
//shown for the sake of how to get these dimensions.
unsigned z_idx = blockIdx.z % array_zz_z;
unsigned w_idx = blockIdx.z / array_zz_z;
//we already know this part of the indexing calculation.
unsigned out_idx_zw = blockIdx.z * (array_zz_xx_y * array_xx_z);
// since our input array is actually 3D, this is a different calcualation
unsigned array_zz_zw = blockIdx.z * (array_zz_xx_y)
//ensures if our launch dimensions don't exactly match our input size, we don't
//accidently access out of bound memory, while branching can be bad, this isn't
// because 99.999% of the time no branch will occur and our instruction pointer
//will be the same per warp, meaning virtually zero cost.
if(x_idx < array_xx_x){
//moving over y axis to coalesce memory accesses in the x dimension per warp.
for(int i = 0; i < ilp_factor; ++i){
//need to also check y, these checks are virtually cost-less
// because memory access dominates time in such simple calculations,
// and arithmetic will be hidden by overlapping execution
if((y_idx+i) < array_zz_xx_y){
//splitting up calculation for simplicity sake
out_array_idx = out_idx_zw+(y_idx+i)*array_xx_x + x_idx;
array_zz_idx = array_zz_zw + (y_idx+i);
array_xx_idx = ((y_idx+i) * array_xx_x) + x_idx;
//actual final output.
out_array[out_array_idx] = array_zz[array_zz_idx] * array_xx[array_xx_idx];
}
}
}
}
Вам нужно будет сделать размеры запуска примерно такими:
thread_dim = (work_tile_dim, work_tile_dim/ilp_factor) # (32,8)
y_dim = xx0.shape[0]
x_dim = xx0.shape[1]
wz_dim = zz0.shape[0] * zz0.shape[1]
block_dim = (x_dim/work_tile_dim, y_dim/work_tile_dim, wz_dim)
И есть еще несколько способов оптимизации, которыми вы можете воспользоваться:
хранит глобальные обращения к памяти в рабочей плитке в разделяемой памяти внутри ядра, это гарантирует, что обращения к zz0s "y", но на самом деле измерения x объединяются, когда помещаются в разделяемую память, увеличивая производительность, а затем получая доступ из разделяемой памяти (где объединение не имеет значения, но конфликты банков имеют значение). См. здесь о том, как бороться с такого рода банковским конфликтом.
вместо вычисления формулы Эйлера и расширения двойного числа до сложного двойного, разверните его внутри самого ядра, используйте sincos(-x, &out_sin, &out_cos) для достижения того же результата, но используя намного меньшую пропускную способность памяти (см. здесь ).
Но обратите внимание, что даже выполнение этого, скорее всего, не даст вам желаемой производительности (хотя, скорее всего, будет быстрее), если вы не используете GPU более высокого уровня с модулями с двойной точностью, которых нет на большинстве GPU (большинство время эмулируется). Модули с плавающей запятой двойной точности занимают много места, а поскольку графические процессоры используются для графики, они практически не используют двойную точность. Если вам нужна более высокая точность, чем с плавающей запятой, но вы хотите использовать преимущества аппаратного обеспечения с плавающей запятой без удвоения пропускной способности от 1/8 до 1/32, вы можете использовать методы, используемые в в этом ответе , для достижения это на GPU, приближая вас к пропускной способности от 1/2 до 1/3.