Ответ 1

Я связался с технической поддержкой Mathworks, и Rylan наконец прояснил эту проблему. (Спасибо, Райлан!) Его полный ответ ниже. Проблема с функцией vs script связана с некоторыми оптимизациями. Matlab автоматически применяется к функциям (но не к скриптам), которые не работают должным образом.

Ответ Rylan:

Благодарим вас за терпение по этому вопросу. Я проконсультировался с разработчиками вычислительных машин MATLAB, чтобы понять это лучше.

Эта проблема вызвана внутренними оптимизациями, выполняемыми MATLAB, когда сталкиваются с некоторыми конкретными операциями, такими как умножение матрицы матриц и транспонирование. Некоторые из этих оптимизаций могут быть активированы специально при выполнении функции MATLAB (или анонимной функции), а не в скрипте.

Когда ваш исходный код выполнялся из сценария, определенная оптимизация транспонирования матрицы не выполняется, что приводит к тому, что выражение "res2" выполняется быстрее, чем выражение "res1":

  n = 2000;
  a=gpuArray(sprand(n,n,0.01)); 
  b=gpuArray(rand(n));

  tic;res1=a*b*a;wait(gpuDevice);toc                                         % Elapsed time is 0.884099 seconds.
  tic;res2=transpose(transpose(a)*transpose(a*b));wait(gpuDevice);toc        % Elapsed time is 0.068855 seconds.

Однако, когда вышеуказанный код помещается в файл функции MATLAB, выполняется дополнительная оптимизация времени транспозиции матрицы, которая заставляет выражение "res2" проходить через другой путь кода (и другой вызов функции библиотеки CUDA) по сравнению с той же строкой, что и вызванный из сценария. Поэтому эта оптимизация генерирует более медленные результаты для строки "res2" при вызове из файла функций.

Чтобы избежать возникновения этой проблемы в файле функций, операции переноса и умножения необходимо разбить так, чтобы MATLAB не применял эту оптимизацию. Разделение каждого предложения в выражении "res2" представляется достаточным для этого:

  tic;i1=transpose(a);i2=transpose(a*b);res3=transpose(i1*i2);wait(gpuDevice);toc      % Elapsed time is 0.066446 seconds.

В приведенной выше строке "res3" генерируется из двух промежуточных матриц: "i1" и "i2". Производительность (в моей системе), по-видимому, соответствует показателю выражения "res2" при выполнении из сценария; кроме того, выражение 'res3' также показывает аналогичную производительность при исполнении из файла функции MATLAB. Однако обратите внимание, что для хранения транспонированной копии исходного массива может использоваться дополнительная память. Пожалуйста, дайте мне знать, если вы видите разные характеристики производительности в своей системе, и я могу исследовать это дальше.

Кроме того, операция "res3" показывает более высокую производительность при измерении с помощью функции "gputimeit". Дополнительную информацию об этом см. В прилагаемом файле testScript2.m. Я также добавил "test_v2.m", который является модификацией функции "test.m" в вашей статье "Переполнение стека".

Спасибо, что сообщили мне об этом. Я хотел бы извиниться за любые неудобства, вызванные этой проблемой. Я создал внутренний отчет об ошибке, чтобы сообщить разработчикам MATLAB об этом поведении. Они могут обеспечить исправление для этого в будущей версии MATLAB.

Поскольку у вас возникли дополнительные вопросы о сравнении производительности графического процессора с использованием "gputimeit" и "tic" и "toc", я просто хотел предложить одно предложение, о котором упомянули ранее разработчики компьютеров MATLAB. Как правило, полезно также "wait (gpuDevice)" перед операциями "tic", чтобы гарантировать, что операции GPU из предыдущих строк не перекрываются при измерении для следующей строки. Например, в следующих строках:

  b=gpuArray(rand(n));
  tic; res1=a*b*a; wait(gpuDevice); toc  

если "wait (gpuDevice)" не вызывается до "tic", некоторое время, затраченное на создание массива "b" из предыдущей строки, может перекрываться и подсчитываться за время, затраченное на выполнение выражения "res1". Это было бы предпочтительнее:

  b=gpuArray(rand(n));
  wait(gpuDevice); tic; res1=a*b*a; wait(gpuDevice); toc  

Помимо этого, я не вижу никаких конкретных проблем в том, как вы используете функции "tic" и "toc". Однако обратите внимание, что использование "gputimeit" обычно рекомендуется с использованием "tic" и "toc" непосредственно для профилирования, связанного с GPU.

Я продолжу и закрываю этот случай, но, пожалуйста, дайте мне знать, если у вас возникнут дополнительные вопросы.

%testscript2.m
n = 2000;
a = gpuArray(sprand(n, n, 0.01)); 
b = gpuArray(rand(n)); 

gputimeit(@()transpose_mult_fun(a, b))
gputimeit(@()transpose_mult_fun_2(a, b))

function out = transpose_mult_fun(in1, in2)

i1 = transpose(in1);
i2 = transpose(in1*in2);

out = transpose(i1*i2);

end

function out = transpose_mult_fun_2(in1, in2)

out = transpose(transpose(in1)*transpose(in1*in2));

end

,

function test_v2

clear

%% transposed expression
n = 2000;
rng('default');rng(1);
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);

tic;
c1 = gather(transpose( transpose(a) * transpose(a * b) ));

disp(['time for (a''*(a*b)'')'': ' , num2str(toc),'s'])

clearvars -except c1

%% non-transposed expression
rng('default');
rng(1)
n = 2000;
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);

tic;
c2 = gather(a * b * a);

disp(['time for a*b*a: ' , num2str(toc),'s'])
disp(['error = ',num2str(max(max(abs(c1-c2))))])

%% sliced equivalent
rng('default');
rng(1)
n = 2000;
a = sprand(n, n, 0.1);
b = rand(n, n);
a = gpuArray(a);
b = gpuArray(b);

tic;
intermediate1 = transpose(a);
intermediate2 = transpose(a * b);
c3 = gather(transpose( intermediate1 * intermediate2 ));

disp(['time for split equivalent: ' , num2str(toc),'s'])
disp(['error = ',num2str(max(max(abs(c1-c3))))])

end

Ответ 2

Кажется, существует проблема с умножением двух разреженных матриц на GPU. время для разреженной полной матрицы более чем в 1000 раз быстрее, чем разреженное по разреженной. Простой пример:

str={'sparse*sparse','sparse*full'};
for ii=1:2
    rng(1);
    a=sprand(3000,3000,0.1);
    b=sprand(3000,3000,0.1);
    if ii==2
        b=full(b);
    end
    a=gpuArray(a);
    b=gpuArray(b);
    tic
    c=a*b;
    disp(['time for ',str{ii},': ' , num2str(toc),'s'])
end

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

str={'a*b*a','a*b*full(a)'};
for ii=1:2
    %rng('default');
    rng(1)
    a=sprand(3000,3000,0.1);
    b=rand(3000,3000);
    rng(1)
    c=sprand(3000,3000,0.1);
    if ii==2
        c=full(c);
    end
    a=gpuArray(a);
    b=gpuArray(b);
    c=gpuArray(c);
    tic;
    c1{ii}=a*b*c;
    disp(['time for ',str{ii},': ' , num2str(toc),'s'])
end
disp(['error = ',num2str(max(max(abs(c1{1}-c1{2}))))])

Возможно, я ошибаюсь, но мой вывод состоит в том, что a * b * a включает в себя умножение двух разреженных матриц (a и a снова) и плохо обрабатывается, тогда как использование метода transpose() делит процесс на двухэтапное умножение ни в одном из которые есть две разреженные матрицы.

Ответ 3

EDIT 2 Возможно, я был прав, см. Другой ответ

EDIT: Они используют MAGMA, который является основным столбцом. Мой ответ не выполняется, однако я оставлю его здесь некоторое время, если он может помочь взломать это странное поведение.

Следующий ответ неверен

Это моя догадка, я не могу сказать вам 100%, не зная кода под капотом MATLAB.

Гипотеза: параллельный компьютерный код MATLAB использует библиотеки CUDA, а не их собственные.

Важная информация

• MATLAB является основным столбцом, а CUDA - основной.
• Нет таких вещей, как 2D-матрицы, только 1D-матрицы с двумя индексами

Почему это имеет значение? Хорошо, потому что CUDA - это высоко оптимизированный код, который использует структуру памяти, чтобы максимизировать кеширование на ядро (самая медленная работа на графических процессорах - чтение памяти). Это означает, что стандартный код умножения матрицы CUDA будет использовать порядок чтения в памяти, чтобы убедиться, что они смежны. Однако, что соседняя память в строке-майне не находится в столбце.

Итак, есть 2 решения для этого, поскольку кто-то пишет программное обеспечение

1. Напишите свои собственные библиотеки алгебр-столбцов в CUDA
2. Возьмите каждый вход/выход из MATLAB и перенесите его (т.е. Преобразовать из столбца в майор)

Они сделали пункт 2 и предположили, что есть умный JIT-компилятор для MATLAB параллельной обработки инструментария (разумное предположение), для второго случая он принимает a и b, переносит их, выполняет математику и переносит вывод при gather,

В первом случае, однако, вам уже не нужно транспонировать вывод, поскольку он внутренне уже транспонирован и JIT ловит его, поэтому вместо вызова gather(transpose( XX )) он просто пропускает выходную транспозицию сбоку. То же самое с transpose(a*b). Обратите внимание, что transpose(a*b)=transpose(b)*transpose(a), поэтому внезапно не требуется транспонирование (все они внутренне пропущены). Транспозиция - дорогостоящая операция.

Действительно, здесь есть странная вещь: превращение функции в код внезапно замедляет работу. Мое лучшее предположение заключается в том, что, поскольку JIT ведет себя по-разному в разных ситуациях, он не улавливает все эти вещи transpose внутри и просто выполняет все операции в любом случае, теряя скорость.

Интересное наблюдение: в процессоре требуется то же самое время, что и GPU, чтобы сделать a*b*a на моем ПК.