Главная
Блог разработчиков phpBB
 
+ 17 предустановленных модов
+ SEO-оптимизация форума
+ авторизация через соц. сети
+ защита от спама

Пишем numpy-модуль для убыстрения математических функций с поддержкой SIMD-инструкций

Anna | 16.06.2014 | нет комментариев
Пакеты numpy и scipy предоставляют красивые вероятности для стремительного решения разных вычислительных задач. Доктрина многофункциональных функций (ufunc), работающих как со скалярными значениями, так и с массивами разных размерностей, разрешает получить высокую эффективность при сохранении принадлежащей языку Python простоты и элегантности. Универсальная функция обыкновенно применяются для выполнения одной операции над огромным массивом данных, что безупречно подходит для оптимизации с поддержкой SIMD-инструкций, впрочем мне не удалось обнаружить готового решения, основанного на свободном программном обеспечении и разрешающего применять SIMD для вычисления в numpy таких математических функций, как синус, косинус и экспонента. Реализовывать алгорифмы вычисления этих функций с нуля вовсе не хотелось, но к счастью в интернете нашлось несколько свободных библиотек на языке «С». Одолев лень сомнения, я решил написать личный numpy-модуль, предлагающий универсальные функции для синуса, косинуса и экспоненты. За подробностями и итогами тестов добродушно пожаловать под кат.

Немножко о SIMD-инструкциях

SIMD-инструкции разрешают единовременно изготавливать один и тот же комплект операций над несколькими числами, записанными в одном регистре. Таким образом, за один такт дозволено обрабатывать сразу несколько чисел и допустимо увеличить продуктивность в несколько раз.
Скажем, комплект SIMD-инструкций Advanced Vector Extensions (AVX) разрешает исполнять операции с 256-битными регистрами, всякий из которых может включать восемь 32-разрядных чисел с плавающей точкой (числа одинарной точности), либо четыре 64-разрядных (числа двойственный точности). Комплект операций достаточно скромный, в основном это сложение, вычитание, умножение и деление. Точная и стремительная реализация тригонометрических функций с поддержкой этих операций — задача нетривиальная и, Дабы не изобретать велосипед, стоит применять какую-нибудь готовую библиотеку. Помимо проприетарной Intel MKL(которая теснее может трудиться с numpy) вариантов нашлось каждого три (раздватри).
1-й вариант представляет собой заголовочный файл с дюже скромным комплектом функций, фактически без документации и тестов. 2-й вариант — С библиотека vecmathlib, которая у меня отчего-то настойчиво отказывалась компилироваться, невзирая на применение рекомендованного компилятора GCC-4.7. Вариант многообещающий, но пока еще схоже сырой. 3-й вариант — библиотеку SLEEF — удалось обнаружить именно вследствие vecmathlib, которая использует его кодовую базу. Вариант мне сразу понравился простотой и ясностью кода, а также обилием тестов.

Маленький тест для мотивации

Дабы получить довольную мотивацию для написания модуля, а заодно разобраться с применением SLEEF, я решил сравнить скорость вычисления синуса на языке «C» при применении SLEEF со стандартной math.h. Безусловно, речь идет о поэлементном вычислении синуса для большого массива данных.
К сожалению, документации и примеров в SLEEF фактически нет, но но есть абсолютно ясно написанные тесты, так что разобраться с применением библиотеки было не трудно. Начальный код SLEEF состоит из четырех директорий: javapurecsimd и tester. Помимо этого, там лежит файл README с коротким изложением библиотеки и всеобщий Makefile, дергающий Makefile из перечисленных директорий. Меня безусловно огромнее каждого заинтересовала директория simd, в которой, как дозволено додуматься из наименования, содержались оптимизированные с поддержкой SIMD-инструкций функции.
Из Makefile в директории simd ясно, что поддерживаются 4 варианта SIMD-инструкций: SSE2AVXAVX2 иFMA4. Прототипы функций определены в файле sleefsimd.h, а необходимый комплект инструкций выбирается при компиляции с поддержкой флагов -DENABLE_SSE2-DENABLE_AVX-DENABLE_AVX2 либо -DENABLE_FMA4. Makefile собирает исполняемые файлы для тестирования функций с применением всякого из комплектов инструкций: iutsse2iutavxiutavx2 либо iutfma4. Эти файлы вызываются из многофункциональной программы tester (из директории tester) и осуществляют выполнение получаемых от tester команд. Реализация команд находится в файле iut.c, откуда становится явственным метод применения библиотеки.

Функция тестирования синуса из файла simd/iut.c начального кода SLEEF

double xxsin(double d) { 
  double s[VECTLENDP];
  int i;
  for(i=0;i<VECTLENDP;i  ) {
    s[i] = random()/(double)RAND_MAX*20000-10000;
  }
  int idx = random() & (VECTLENDP-1);
  s[idx] = d; 

  vdouble a = vloadu(s);
  a = xsin(a);
  vstoreu(s, a);

  return s[idx];
}

Для чисел двойственный точности (double) необходимо определить массив длиной VECTLENDP, заполнить доводами волнующей нас функции и передать в функцию vloadu, которая скопирует их в нужное для применения SIMD место и вернет значение типа vdouble. Значение vdouble мы передаем функции xsin, которая вычисляет значение синуса сразу для всех VECLENDP доводов и возвращает вновь же vdouble. Итог распаковывается в массив из double с поддержкой функции vstoreu.
Для желающих проверить SLEEF на своей машине привожу полный начальный код программы, которую написал для оценки возможного убыстрения от применения SIMD с поддержкой SLEEF.

Программа для оценки скорости вычисления синуса с поддержкой SLEEF

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "sleefsimd.h"

#define TESTSIZE (VECTLENDP*10000000)
double s[TESTSIZE];
double r1[TESTSIZE];
double r2[TESTSIZE];
#define COUNT 10
int main(int argc, char *argv[])
{
	int k, i;
    clock_t t1, t2;
    double time1, time2;
    double max, rmax;

    srandom(time(NULL));
	for(i = 0; i < TESTSIZE; i  ) 
    {
		s[i] = random()/(double)RAND_MAX*20000-10000;
	}

	printf("Testing sin, %d valuesn", TESTSIZE*COUNT);
	t1 = clock();
	for(k = 0; k < COUNT; k  )
	{
		for(i = 0; i < TESTSIZE; i  ) 
        {
			r1[i] = sin(s[i]);
		}
	}
    t2 = clock();
    time1 = (double)(t2 - t1)/CLOCKS_PER_SEC;
	printf("Finish sin, spent time = %lg secnn", time1);

	printf("Testing xsinn");
	t1 = clock();
	for(k = 0; k < COUNT; k  )
	{
		for(i = 0; i < TESTSIZE; i  = VECTLENDP) 
        {
			vdouble a = vloadu(s i);
			a = xsin(a);
			vstoreu(r2 i, a);
		}
	}
    t2 = clock();
    time2 = (double)(t2 - t1)/CLOCKS_PER_SEC;
	printf("Finish xsin, spent time = %lg secnn", time2);

	printf("Speed ratio: %lfn", time1/time2);

	max = r1[0] - r2[0];
    rmax = (r1[0] - r2[0])/r1[0];
	for(i = 0; i < TESTSIZE; i  ) 
    {
		double delta = (r1[i] - r2[i]);
		if(abs(delta) > abs(max)) max = delta;
        delta = (r1[i] - r2[i])/r1[i];
		if(abs(delta) > abs(max)) rmax = delta;
	}

	printf("Max absolute delta: %lg, relative delta %lgn", max, rmax);
	return 0;
}

Самый продвинутый из поддерживаемых на моем компьютере комплект команд — это AVX, следственно я компилировал программу (записанную в файл simd/speedtest.c в исходниках SLEEF) дальнейшей командой:

gcc -O3 -Wall -Wno-unused -Wno-attributes -DENABLE_AVX -mavx speedtest.c sleefsimddp.c sleefsimdsp.c -o speedtest -lm

Я ждал убыстрения приблизительно в 4 раза, но итог превзошел все мои ожидания:

Testing sin, 400000000 values
Finish sin, spent time = 14.95 sec

Testing xsin
Finish xsin, spent time = 1.31 sec

Speed ratio: 11.412214
Max absolute delta: 5.55112e-17, relative delta 1.58441e-16

Убыстрение больше чем в 10 раз, при относительной погрешности вычисления менее 2·10-16 (порядка точности самого double), на одном ядре процессора! Безусловно, в настоящем приложении убыстрение будет поменьше, но мотивации для написания своего numpy-модуля теснее достатоПри регистрации многофункциональной функции мы указываем поддерживаемые типы и для всякой комбинации входных и выходных типов передаем указатель на С-функцию, которая будет с этой комбинацией трудиться:

static PyUFuncGenericFunction funcs[] = {&double_xsin};
static char types[] = {NPY_DOUBLE, NPY_DOUBLE};
static void *data[] = {NULL};

В массиве types указываются как входные, так и выходные типы, следственно он длиннее массивов funcs иdata. Массив указателей data разрешает для всякой комбинации типов указать свой добавочный параметр, тот, что будет передан в C-фунцкию в качестве последнего довода void* data. В частности, это дозволено применять для реализации различных многофункциональных функций с поддержкой одной C-функции.
Дабы зарегистрировать нашу универсальную функцию, необходимо вызвать PyUFunc_FromFuncAndData и передать ей описанные выше массивы (funcsdata и types), а также число входных и выходных доводов многофункциональной функции, число поддерживаемых комбинаций типов, имя функции в модуле и строку документации.

Полный начальный код модуля

#include "Python.h"
#include "numpy/ndarraytypes.h"
#include "numpy/ufuncobject.h"
#include "numpy/npy_3kcompat.h"
#include "sleef/sleefsimd.h"

/* The loop definition must precede the PyMODINIT_FUNC. */
static void double_xsin(char **args, npy_intp *dimensions,
                            npy_intp* steps, void* data)
{
    npy_intp i;
    npy_intp n = dimensions[0];
    char *in = args[0], *out = args[1];
    npy_intp in_step = steps[0], out_step = steps[1];
    double tmp[VECTLENDP];
    vdouble a;
    int slow_n = n % VECTLENDP;
    if(in_step != sizeof(double) || out_step != sizeof(double))
        slow_n = n;
    for(i = 0; i < slow_n; i  = VECTLENDP)
    {
        int j;
        for(j = 0; j < VECTLENDP && i   j < slow_n; j  )
        {
            tmp[j] = *(double *)in;
            in  = in_step;            
        }
        a = vloadu(tmp);
        a = xsin(a);
        vstoreu(tmp, a);
        for(j = 0; j < VECTLENDP && i   j < slow_n; j  )
        {
            *(double *)out = tmp[j];
            out  = out_step;
        }        
    }
    if(n > slow_n)
    {
        double *in_array = (double *)in;
        double *out_array = (double *)out;
        for(i = 0; i < n - slow_n; i  = VECTLENDP)
        {
            a = vloadu(in_array   i);
        	a = xsin(a);
	        vstoreu(out_array   i, a);    
        }
    }
}

static PyMethodDef AvxmathMethods[] = {
        {NULL, NULL, 0, NULL}
};

static PyUFuncGenericFunction funcs[1] = {&double_xsin};
static char types[] = {NPY_DOUBLE, NPY_DOUBLE};
static void *data[] = {NULL};

void register_xsin(PyObject *module)
{
    PyObject *xsin, *d;
    import_array();
    import_umath();

    xsin = PyUFunc_FromFuncAndData(funcs, data, types, 1, 1, 1,
                                    PyUFunc_None, "sin",
                                    "AVX-accelerated sine calculation", 0);
    d = PyModule_GetDict(module);
    PyDict_SetItemString(d, "sin", xsin);
    Py_DECREF(xsin);
}

#if PY_VERSION_HEX >= 0x03000000
static struct PyModuleDef moduledef = {
    PyModuleDef_HEAD_INIT,
    "avxmath",
    NULL,
    -1,
    AvxmathMethods,
    NULL,
    NULL,
    NULL,
    NULL
};

PyMODINIT_FUNC PyInit_avxmath(void)
{
    PyObject *m;
    m = PyModule_Create(&moduledef);
    if (!m) {
        return NULL;
    }
    register_xsin(m);
    return m;
}
#else
PyMODINIT_FUNC initavxmath(void)
{
    PyObject *m;
    m = Py_InitModule("avxmath", AvxmathMethods);
    if (m == NULL) {
        return;
    }
    register_xsin(m);
}
#endif

Для сборки модуля я применял типовой setup.py из документации, заменив наименование модуля и добавив C-файлы библиотеки SLEEF, флаги компилятора и линковщика. Приведенный выше начальный код модуля я сберег рядом с setup.py в файле с именем avxmath.c, директорию simd из начального кода SLEEF переименовал в sleef и также положил рядом с setup.py.

Файл setup.py для сборки модуля avxmath

def configuration(parent_package='', top_path=None):
    import numpy
    from numpy.distutils.misc_util import Configuration

    config = Configuration('',
                           parent_package,
                           top_path)
    config.add_extension('avxmath', ['avxmath.c', 'sleef/sleefsimddp.c', 'sleef/sleefsimdsp.c'], 
                         extra_compile_args=['-O3','-Wall', '-Wno-unused', '-Wno-attributes', '-DENABLE_AVX','-mavx'], 
                         extra_link_args=['-lm'])

    return config

if __name__ == "__main__":
    from numpy.distutils.core import setup
    setup(configuration=configuration)

Для компиляции без установки в систему необходимо исполнить команду python setup.py build_ext --inplace, итогом удачного выполнения которой должен стать готовый модуль в файле avxmath.so. Сейчас дозволено проверить работоспособность нашей функции. Запускаем python в той же директории, где находится файл avxmath.so, и проверяем:

>>> from numpy import array, pi
>>> import avxmath
>>> avxmath.sin(0)
0.0
>>> avxmath.sin(pi)
1.2246467991473532e-16
>>> avxmath.sin([0, pi/2, pi, 3*pi/2, 2*pi])
array([  0.00000000e 00,   1.00000000e 00,   1.22464680e-16,
        -1.00000000e 00,  -2.44929360e-16])
>>> 

Удостоверясь, что модуль avxmath импортируется и работает без ошибок, дозволено сделать маленький тест продуктивности и точности новой функции.

Программа проверки функции sin модуля avxmath и итог ее выполнения

import numpy 
import avxmath
import time
from numpy import random, pi

COUNT=10

x = 2e4*random.random(40000000) - 1e4
t = time.clock()
for i in xrange(COUNT):
    y1 = numpy.sin(x)
duration1 = time.clock() - t
print "numpy.sin %f sec" % duration1

t = time.clock()
for i in xrange(COUNT):
    y2 = avxmath.sin(x)
duration2 = time.clock() - t
print "avxmath.sin %f sec" % duration2

delta = y2 - y1
rdelta = delta/y1
print "max absolute difference is %lg, relative %lg" % (
        delta[abs(delta).argmax()], rdelta[abs(rdelta).argmax()])
print "speedup is %lg" % (duration1/duration2)
numpy.sin 15.510000 sec
avxmath.sin 2.260000 sec
max absolute difference is 2.22045e-16, relative 2.63873e-16
speedup is 6.86283

Выходит, мы получили убыстрение больше чем в 6 раз при точности вычислений не дрянней 3·10-16! Заменив вызовы функции xsin на примитивное копирование памяти, несложно удостовериться, что убыстрение в 10 раз не получилось из-за того, что около 1 секунды из полученных нами 2.26 секунд выполнения ушло на убыточные расходы. Подобно, заменив функцию xsin на обыкновенный синус из math.h, мы найдем, что времена вычислений с поддержкой avxmath.sin и numpy.sin в нашем тесте с высокой точностью совпадут.
Таким образом, с поддержкой SIMD инструкций дозволено добиться существенного убыстрения исполняемых с поддержкой numpy и scipy вычислений, легко заменив импорты обыкновенных функций на оптимизированные. Ну и безусловно начальные коды несколько расширенного по сопоставлению с данной статьей модуля avxmathдоступны на Github по ссылке.

Оригинальная статья на http://habrahabr.ru/post/198916/

Источник: programmingmaster.ru
Оставить комментарий
Форум phpBB, русская поддержка форума phpBB
Рейтинг@Mail.ru 2008 - 2017 © BB3x.ru - русская поддержка форума phpBB