7. Использование оптимизированных библиотек¶
7.1. Общие сведения¶
Все производители современных процессоров разрабатывают и поставляют пользователям высокопроизводительные библиотеки, обеспечивающие скорость работы, близкую к максимальной для данного типа процессоров. Примерами таких библиотек являются библиотеки IPP/MKL фирмы Intel, библиотеки mediaLib/Perflib фирмы Oracle, библиотеки ACML/APL фирмы AMD.
Для процессоров линии «Эльбрус» также была разработана библиотека EML - высокопроизводительная математическая и мультимедийная библиотека, представляющая из себя набор разнообразных функций для обработки сигналов, изображений, видео, математических вычислений.
7.2. Состав¶
Библиотека eml состоит из следующих разделов:
- Ядро (Core)
выделение и освобождение памяти, номер версии и статус.
- Вектор (Vector)
различные операции над векторами: арифметические, логические, преобразование типов, математические функции, статистика.
- Сигналы (Signal)
цифровая обработка сигналов: конволюция, фильтрация, усиление, генерация, быстрые преобразования Фурье и Хартли.
- Изображение (Image)
создание и заполнение изображений, арифметические операции, фильтрация, геометрические и цветовые преобразования, ДПФ.
- Линейная Алгебра (Algebra)
стандартные пакеты работы с матрицами и векторами BLAS 1/2/3, LAPACK.
- Видео (Video)
обработка видео: интерполяция, усреднение, оценка движения, цветовые преобразования, ДКП, квантизация.
- Графика (Graphics)
рисование/закрашивание точек/линий/треугольников/прямоугольников/ полигонов/дуг/окружностей/эллипсов, закрашивание/перекрашивание области.
- Объем (Volume)
бросание параллельных/произвольных лучей c интерполяцией, линейное масштабирование вокселей, поиск максимальных значений на луче.
7.3. Информационная система¶
Полная документация поставляется вместе с самой библиотекой. Она находится в файле
/opt/mcst/doc/eml/index.html.
Также начинать поиск справочной информации можно в /opt/mcst/doc/eml/annotatid.html.
7.4. Примеры использования¶
Рассмотрим простейшие варианты применения библиотеки eml.
7.4.1. Умножение векторов¶
Самое простое перемножение без eml. Компилируем пример с оптимизацией -O3.
Размер вектора и количество повторений выбраны достаточно большими, чтобы показать вариант, когда расчёты не помещаются в кэш-память.
Пример ниже стоит скомпилировать на ВК «Эльбрус» и на машине с архитектурой х86.
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define N 100000
#define M 1000000
int main()
{
int i,j;
double s=0.0;
double *A;
double *B;
double *C;
A = (double*)malloc(N * sizeof(double));
B = (double*)malloc(N * sizeof(double));
C = (double*)malloc(N * sizeof(double));
srand(time(NULL));
for (i = 0; i < N; i++)
{
A[i] = (double)(rand()%10000)/(double)10000+(double)(rand()%10000);
B[i] = (double)(rand()%10000)/(double)10000+(double)(rand()%10000);
}
for (j=0; j<M; j++)
for(i = 0; i < N; i++)
C[i] = A[i] * B[i];
s+=C[0];
free(A);
free(B);
free(C);
return (int)s;
}
Пример с использованием библиотеки eml будет отличаться двумя строчками, помеченными звездочкой (* после комментариев //). Чтобы скомпилировать пример, надо добавить в сборку линковку библиотеки eml:
gcc -O3 -o eml mul_vector_eml.c -leml
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <eml/eml.h>
#define N 100000
#define M 1000000
int main()
{
int i,j;
double s=0.0;
double *A;
double *B;
double *C;
A = (double*)malloc(N * sizeof(double));
B = (double*)malloc(N * sizeof(double));
C = (double*)malloc(N * sizeof(double));
srand(time(NULL));
for (i = 0; i < N; i++)
{
A[i] = (double)(rand()%10000)/(double)10000+(double)(rand()%10000);
B[i] = (double)(rand()%10000)/(double)10000+(double)(rand()%10000);
}
for (j=0; j<M; j++)
eml_Vector_Mul_64F(A,B,C,N);
s+=C[0];
free(A);
free(B);
free(C);
return (int)s;
}
Проверить корректность работы можно путем вывода вектора и результата на печать:
printf("vector A\n");
for (i = 0; i < N; i++)
{
printf("%F ", A[i]);
}
printf("\nvector B\n");
for (i = 0; i < N; i++)
{
printf("%f ", B[i]);
}
printf("\nthe result\n");
for (i = 0; i < N; i++)
{
printf("%f ", C[i]);
}
printf("\n");
Результаты приведены в таблице Время выполнения программы умножения векторов.
Интел |
Эльбрус |
Эльбрус + eml |
|---|---|---|
real 64.90 |
real 163.65 |
real 63.71 |
user 64.89 |
user 163.54 |
user 63.66 |
sys 0.00 |
sys 0.01 |
sys 0.01 |
7.4.2. Умножение матриц¶
Преимущество использования библиотеки eml более явно проявляется на примере умножения матриц.
Для eml критично, чтобы матрица была объявлена как линейный массив, через A[i*N+j].
Если удобнее обращаться через A[i][j], тогда нужно объявлять массив как A[M][N].
Объявление double **A приведет к нелинеаризованному представлению.
Линеаризованный массив выглядит следующим образом:
A11 |
А12 |
А13 |
А14 |
А21 |
А22 |
А23 |
А24 |
А31 |
… |
Общая функция умножения матриц dgemm¶
Параметры:
Order признак расположения матрицы в памяти по строкам или по столбцам.
TransA признак транспонирования матрицы A.
TransB признак транспонирования матрицы B.
M число строк матрицы A.
N число столбцов матрицы B.
K число столбцов матрицы A.
alpha, beta входные константы.
A входная матрица общего вида.
lda leading Array Dimension - ведущая размерность массива A.
B входная матрица общего вида.
ldb leading Array Dimension - ведущая размерность массива B.
C выходная матрица общего вида.
ldc leading Array Dimension - ведущая размерность массива C.
В EML, Blas, Lapack:
Матрица А¶
Пример умножения матриц без использования eml:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define N 1000
#define M 100
int main()
{
int i,j,k,l;
double s = 0.0;
double *A = (double*)malloc(N * N * sizeof(double));
double *B = (double*)malloc(N * N * sizeof(double));
double *C = (double*)malloc(N * N * sizeof(double));
srand(time(NULL));
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
{
A[i * N + j] = (double)(rand()%10000)/(double)10000+(double)(rand()%10000);
B[i * N +j] = (double)(rand()%10000)/(double)10000+(double)(rand()%10000);
}
for (l=0;l<M;l++)
{
for(i = 0; i < N; i++)
for(j = 0; j < N; j++)
{
for(k = 0; k < N; k++)
C[i*N+j] += A[i*N+k] * B[k*N+j];
}
}
/*
printf("matrix A\n");
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
printf("%f ", A[i * N + j]);
printf("\n");
}
printf("\nmatrix B\n");
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
printf("%f ", B[i * N + j]);
printf("\n");
}
printf("\nthe result of multiplying\n");
for (i = 0; i < N; i++)
{
for (j = 0; j < N; j++)
printf("%3f ", C[i * N + j]);
printf("\n");
}
*/
s+=C[0];
free(A);
free(B);
free(C);
return (int)s;
}
Пример умножения матриц с использованием eml:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <eml.h>
#define N 1000
#define M 100
int main()
{
int i,j,k;
double s=0.0;
double *A = (double*)malloc(N * N * sizeof(double));
double *B = (double*)malloc(N * N * sizeof(double));
double *C = (double*)malloc(N * N * sizeof(double));
srand(time(NULL));
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
{
A[i * N + j] = (double)(rand()%10000)/(double)10000+(double)(rand()%10000);
B[i * N +j] = (double)(rand()%10000)/(double)10000+(double)(rand()%10000);
}
for (i=0; i<M;i++)
eml_Algebra_GEMM_64F(EML_MATRIX_ROW_MAJOR, EML_MATRIX_NO_TRANS,
EML_MATRIX_NO_TRANS, N, N, N, 1, A, N, B, N, 0, C, N);
s+=C[0];
free(A);
free(B);
free(C);
return (int)s;
}
Результаты приведены в таблице Время выполнения программы умножения матриц.
Интел |
Эльбрус |
Эльбрус + eml |
|---|---|---|
real 316.88 |
real 1207.73 |
real 14.72 |
user 316.82 |
user 1206.92 |
user 14.69 |
sys 0.01 |
sys 0.03 |
sys 0.02 |
Рассмотрим шаги, которые могут улучшить результат умножения матриц. Данные приёмы уже реализованы в библиотеке eml для повышения производительности при умножении матриц.
Классическое гнездо циклов для перемножения матриц выглядит следующим образом:
for (i = 0; i < M; i++)
for (j = 0; j < N; j++)
for (k = 0; k < K; k++)
c[i][j] += a[i][k] * b[k][j];
Для того, чтобы избавиться от лишних чтений из памяти, первое, что можно сделать – присваивать значение a[i][k] * b[k][j] не в c[i][j],
а в переменную, хранящуюся на регистре. Тогда наше гнездо будет выглядеть следующим образом:
for (i = 0; i < M; i++)
{
for (j = 0; j < N; j++)
{
s = 0.0;
for (k = 0; k < K; k++)
S += a[i][k] * b[k][j];
c[i][j] = s;
}
}
Это улучшит время приведенного примера с теми же параметрами:
real 448.55
user 448.24
sys 0.02
Для расчета такого гнезда требуется 2*M*N*K вещественных операций, в предположении соизмеримости размеров матриц.
Сложность такого алгоритма O(N3).
Для оптимизации внутреннего цикла применимы следующие оптимизации:
k++иk<Nзаменяются на аппаратный счетчик циклов LSR;операции ld из массивов
a[i][k]иb[i][k]заменяются на mova.
После у нас есть два подхода для оптимизации внутреннего цикла:
- Конвейеризация:
ресурсы: 2 mova, 1 fmul, 1 fadd;
рекуррентность: 4 такта (fadd зависит от fadd на предыдущей итерации);
результат: 0.5 flop/cycle;
- loop unroll + балансировка рекуррентности (
-ffast-math), конвейеризация: ресурсы: 16 mova, 8 fmul, 8 fadd;
рекуррентность: 4 тактов;
результат: 4.0 flop/cycle.
- loop unroll + балансировка рекуррентности (
Второй вариант выглядит намного лучше обычной конвейеризации. Время, полученное в результате такого выполнения, будет следующим:
real 167.63
user 167.52
sys 0.01
Цикл, приведенный в примере, будет выглядеть так:
for (l=0;l<M;l++)
{
for(i = 0; i < N; i++)
for(j = 0; j < N; j++)
{
s0=s1=s2=s3=s4=s5=s6=s7=0.0;
for(k = 0; k < N; k+=8)
{
s0 += A[i*N+k] * B[k*N+j];
s1 += A[i*N+k+1] * B[k+1*N+j];
s2 += A[i*N+k+2] * B[k+2*N+j];
s3 += A[i*N+k+3] * B[k+3*N+j];
s4 += A[i*N+k+4] * B[k+4*N+j];
s5 += A[i*N+k+5] * B[k+5*N+j];
s6 += A[i*N+k+6] * B[k+6*N+j];
s7 += A[i*N+k+7] * B[k+7*N+j];
}
C[i*N+j]=s0;
C[i*N+1+j]=s1;
C[i*N+2+j]=s2;
C[i*N+3+j]=s3;
C[i*N+4+j]=s4;
C[i*N+5+j]=s5;
C[i*N+6+j]=s6;
C[i*N+7+j]=s7;
}
}
После применения второго варианта появятся проблемы при дальнейшей оптимизации внутреннего цикла:
for (k = 0; k < N; k+=8)
{
s0 += a[i][k ] * b[k ][j];
s1 += a[i][k+1] * b[k+1][j];
...
s7 += a[i][k+7] * b[k+7][j];
}
Проблема 1: плохая локальность по памяти для массива b.
Решение: транспонирование матрицы B. Сложность О(N2) относительно мала.
Проблема 2: одно чтение происходит для одной флотовой операции.
Решение: переход от оптимизации цикла к оптимизации гнезда циклов.
После транспонирования цикл будет выглядеть так:
for (i = 0; i < M; i+=2)
{
for (j = 0; j < N; j+=2)
{
s00 = s01 = s10 = s11 = 0.0;
for (k = 0; k < K; k++)
{
s00 += a[i ][k] * b[j ][k];
s01 += a[i+1][k] * b[j ][k];
s10 += a[i ][k] * b[j+1][k];
s11 += a[i+1][k] * b[j+1][k];
}
c[i][j ] = s00; c[i+1][j ] = s01;
c[i][j+1] = s10; c[i+1][j+1] = s11;
}
}
Обращаем внимание на тот факт, что чтение a[i][k] не зависит от j;
применение unroll&fuse (unroll&jam) к циклу по j позволит сделать одно чтение a[i][k] для одной раскрученной итерации j.
Аналогично unroll&fuse по i позволит сделать одно чтение b[j][k] для одной раскрученной итерации i.
При оптимизации гнезда циклов важно не упираться в ресурсы аппаратуры: количество вещественных операций, которое мы можем выполнять в одной ШК, количество доступных регистров. Поэтому для архитектуры e8c было бы правильным подобрать параметры следующим образом:
Подбор параметров unroll&fuse по i на 8 и по j на 6 и дальнейшая конвейеризация:
ресурсы: 8+6 = 14 mova, 8*6 = 48 fmul_add;
рекуррентность: 8 тактов (длительность fmul_add).
Успешная конвейеризация в 8 тактов, и соотношение счет/память стало существенно лучше:
12/14 = 0.857 flop/byte (> 0.375).
Результат: 12.0 flop/cycle (!)
При этом гнездо будет выглядеть следующим образом:
for (i = 0; i < M; i+=8)
{
for (j = 0; j < N; j+=6)
{
s00 = ... = s75 = 0.0;
for (k = 0; k < K; k++)
{
s00 += a[i ][k] * b[j ][k];
...
s75 += a[i+7][k] * b[j+5][k];
}
c[i][j] = s00; ... ; c[i+7][j+5] = s75;
}
}