Здравствуйте %username%, данная статья посвящена проблеме увеличения скорости математических вычислений на основе языка python с использованием scipy.weave и openMP.
Многие могут задаться вопросом: «Зачем вообще использовать python для математических вычислений?», но мы не будем отвечать на «вечные» вопросы, как и не будем рассматривать множество других решений данной проблемы, таких как, например, psyco.
Инструменты
Как описано выше, наш инструмент — это библиотека scipy.weave, а также библиотека openMP.scipy — набор библиотек для вычислений в прикладной математике и науке. openMP — открытый стандарт для распараллеливания программ на языках Си, Си++ и Фортран.
Установка пакетов
В Debian-подобных Linux системах необходимо выполнить:
apt-get install python-scipy
apt-get install libgomp1
Метод
Для увеличения скорости вычислений необходимо реализовать «узкую» часть python кода (обычно это цикл, в котором происходят какие-либо действия с матрицей) на C и добавить директивы openMP для распараллеливания.
Пример
Думаю, что нет ничего лучше, чем убедиться в этом способе на примере решения следующей задачи:есть матрица размера n на n, вектор размера n и целое число;
необходимо вычесть из каждой строки матрицы вектор, умноженный на целое число (из симплекс метода).
Python реализация
В python c использование numpy данная задача, не учитывая различных подготовительных операций, как генерация матрицы и прочего, решает в пару строк кода:
# цикл по строкам матрицы, где i - номер строки
# с - целое число, randRow - случайный вектор
for i in xrange(N):
matrix[i,:] -= c*randRow
Генерация случайной матрицы x на y, в нашем случае x=y:
# генерация случайной матрицы x на y
# элементы матрицы - случайные числа от 0 до 99 включ.
def randMat(x,y):
randRaw = lambda a: [ randint(0, 100) for i in xrange(0,a) ]
randConst = lambda x, y: [ randRaw(x) for e in xrange(0,y)]
return array(randConst(x, y))
Реализация scipy.weave без openMP
scipy.weave — часть библиотеки scipy, которая позволяет использовать C/C++ код внутри кода python.
Происходит это следующим образом:
# C код
codeC =
"""
int i = 0;
for(i = 0; i < N*M; i++) {
matrix[0,i] = matrix[0,i] - (c * randRow[i%M]);
}
"""
weave.inline(codeC, ['matrix','c', 'randRow','N', 'M'], compiler = 'gcc')
т.е. сам С код хранится в виде multiline string, а переменные python кода, передаются в С списком, где элементы — одноименные текстовые константы. Также numpy массивы передаются в C не в виде матрицы, а в виде вектора, именно поэтому в коде один цикл, а не два.
Кстати, получившийся С код можно поискать в /tmp/%user%/python2x_intermediate/compiler_x
Реализация scipy.weave с openMP
Теперь уже к добавленной версии необходимо добавить директивы openMP и в вызове inline добавить недостающие параметры, а именно:
# C и openMP код
codeOpenMP =
"""
int i = 0;
omp_set_num_threads(2);
#pragma omp parallel shared(matrix, randRow, c) private(i)
{
#pragma omp for
for(i = 0; i < N*M; i++) {
matrix[0,i] = matrix[0,i] - (c * randRow[i%M]);
}
}
"""
...
weave.inline(codeOpenMP, ['matrix','c', 'randRow','N', 'M'],
extra_compile_args =['-O3 -fopenmp'],
compiler = 'gcc',
libraries=['gomp'],
headers=[''])
Полный исходный код со всеми реализациями можно скачать здесь
Сравнение результатов
Выше изложенный исходный код можно запустить и убедиться в том, что scipy.weave действительно даёт прирост в скорости:
Test on size: 100x100
Pure python: 0.0725984573364
Pure C: 0.303888320923
C plus OpenMP: 0.109100341797
Test - ok
Test on size: 1000x1000
Pure python: 1.00839138031
Pure C: 0.506997108459
C plus OpenMP: 0.333213806152
Test - ok
Test on size: 2000x2000
Pure python: 3.24151515961
Pure C: 2.10800170898
C plus OpenMP: 1.17690563202
Test - ok
Test on size: 3000x3000
Pure python: 5.54490089417
Pure C: 4.61800098419
C plus OpenMP: 2.56960391998
Test - ok
Литература
В написании кода использовались следующие ресурсы:
Документация по scipy.weave
Сайт openMP
