Python / Python, scipy.weave и openMP — разгоняем код

в 20:04, , рубрики: Новости, метки: , ,

Здравствуйте %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

Поделиться

* - обязательные к заполнению поля