Увеличиваем скорость вычислений в Matlab

Добрый день, уважаемые хабровчане!
Предлагаю вашему вниманию статью, на которую меня натолкнул недавний пост Matlab кластер своими руками. В свое время я столкнулся с теми же проблемами что и автор топика, однако для увеличения производительности я не стал уходить в дебри распределенных вычислений, а просто решил максимально оптимизировать свой код. Как это сделать в Matlab, прошу под кат за подробностями.
Статья является описанием жизненного опыта и вольным переводом англоязычной публикации.
Не секрет, что крупномасштабные численные расчеты предъявляют к компьютерному железу немалые требования. Сидя за своей диссертацией, я понял что мой домашний компьютер перестал удовлетворять мои потребности в вычислениях. Так как два месяца назад я сменил компьютер и собрал топовое железо, то вариант со сменой аппаратной части для меня не подходил. Существует мнение, что оптимизация кода — дело не сильно важное, проще просто заняться увеличением расходов на приобретение более современного железа, и это мнение имеет право на существование, но я его не разделяю и поэтому полез в свой код и в специальную литературу.

Matlab относится к интерпретируемым языкам, это казалось бы делает его не подходящим языком для крупномасштабных научных вычислений. Все потенциальные возможности языка реализуются благодаря обширному набору библиотек и операции векторизации. Понятие концепции векторизации занимает центральное место в понимании того, как следует писать эффективный код на Matlab. Все данные хранятся в RAM в виде векторов, поэтому и скорость работы численных алгоритмов в Matlab напрямую зависит от использования операции векторизации.

Ниже будут представлены некоторые основные выдержки из официальной документации Matlab, которые помогут Вам в написании эффективных процедур в Matlab. Можно стать одержимым и тратить возмутительно много человеко-часов рабочего времени на оптимизацию кода, которая поможет выиграть лишь несколько секунд время выполнения программы, но повторюсь, сказав что это была моя производственная необходимость.

Производительность ваших программ может значительно вырасти, если вы будете использовать нижепреведенные приемы:

• Используйте векторные операции, а не циклы;
• Используйте предварительное выделение памяти для векторов и матриц;

Рассмотрим каждый из этих пунктов поподробнее.

Использование векторных операция вместо циклов

Рассмотрим следующий пример:
dx = pi/30;
nx = 1 + 2*pi/dx;
for i = 1:nx
x(i) = (i-1)*dx;
y(i) = sin(3*x(i));
end
Совершенно очевидно, что в нем происходит создание векторов x и y. Однако Matlab выделяет память для переменных «на лету». Во время первого прохода цикла создаются два вектора по одной строке длиной 1. На каждом следующем шаге происходит увеличение длины векторов. Во время каждого прогона цикла будет происходить выделение памяти под переменные Matlab'ом, что отнимает процессорное время.
Предпочтительнее использовать следующий метод задания векторов:
x = 0:pi/30:2*pi
y = sin(3*x);
Первая строка кода создает вектор-строку x, под которую выделяется размер памяти соответствующий ее длине и элементы вектора хранятся в смежных ячейках RAM. Вторая строка кода инициирует создание второй вектор-строки, под которую Matlab сразу же выделяет необходимую память, не возвращаясь больше к этому вопросу.
Помните главное: Matlab больше приспособлен к векторным и матричным вычислениям, именно для этого он и разрабатывался.

Предварительное выделение памяти для векторов и матриц

Хотя Matlab автоматически регулирует размеры векторов и матриц, все-таки лучше предварительно выделить матрицы. Предварительное распределение берет на себя расходы на выделение памяти только один раз, и это гарантирует, что матричные элементы будут храниться в смежных ячейках оперативной памяти (по столбцам).
Рассмотрим следующий код:
dx = pi/30;
nx = 1 + 2*pi/dx;
nx2 = nx/2;

for i = 1:nx2
x(i) = (i-1)*dx;
y(i) = sin(3*x(i));
end

for i = nx2+1:nx
x(i) = (i-1)*dx;
y(i) = sin(5*x(i));
end
Это неоптимизированный код, так как мы люди думающие, то мы знаем размер векторов x и y, и можем воспользоваться этим знанием:
dx = pi/30;
nx = 1 + 2*pi/dx;
nx2 = nx/2;

x = zeros(1,nx);
y = zeros(1,nx);

for i = 1:nx2
x(i) = (i-1)*dx;
y(i) = sin(3*x(i));
end

for i = nx2+1:nx
x(i) = (i-1)*dx;
y(i) = sin(5*x(i));
end
Этот код был преобразован, однако имеются еще кусочки, к которым можно применить правило предварительного распределение памяти. Вектора x (i) =..., и y (i) =… не используют преимущества векторизации.
Оптимальнее всего будет написать вот так:
dx = pi/30;
nx = 1 + 2*pi/dx;
nx2 = nx/2;
x = zeros(1,nx);
y = x;

for i = 1:nx2
x(i) = (i-1)*dx;
y(i) = sin(3*x(i));
end

for i = nx2+1:nx
x(i) = (i-1)*dx;
y(i) = sin(5*x(i));
end
Надеюсь данная статья будет полезна еще кому-либо, в свое время она сэкономила мне немало времени ожидания и помогла успешно выполнить поставленную задачу.