Московский Физико-Технический Институт Оптимизация методов умножения матриц библиотеки линейной алгебры для ВК Эльбрус-3M1 и Эльбрус-90 микро Выполнил: - презентация

1 Московский Физико-Технический Институт Оптимизация методов умножения матриц библиотеки линейной алгебры для ВК Эльбрус-3M1 и Эльбрус-90 микро Выполнил: Чернилевский Д.А. Научный руководитель: Логинов В.Е г.

2 Разделы библиотеки EML.

3 Структура раздела Algebra (библиотека BLAS): y = α x + y y = α Ax + β y C = α AB + β C A и В – матрицы, х и у – вектора, α и β – числа. BLAS:

4 Цель работы Провести анализ существующей реализации библиотеки EML Выявить возможные способы оптимизации кода с применением архитектурных особенностей ВКЭльбрус-90 микро и Эльбрус-3М1 Провести оптимизацию функций умножения матриц, основываясь на проведенном анализе

5 Оптимизации функций умножения матриц для Эльбрус-90 микро Разбиение матриц на блоки Применение оптимизации unroll к исходным текстам Устранение зависимостей между итерациями циклов Использование последовательного чтения/записи из/в память Вынос нелинейных индексов массивов за пределы цикла

6 Разбиение матриц на блоки Причины оптимизации: - Отсутствие устройства предподкачки - Размер L1$: 32kb Решение: - Поблочное умножение *TRMM – Умножение треуг. матриц, GEMM – Умножение прямоуг. матриц

7 for (i = 0; i < M - 3; i += 4) { a0 = (eml_32f *) A; b1 = b0 + ldb; b2 = b1 + ldb; b3 = b2 + ldb; for (j = 0; j < N; j++) { t0 = t1 = t2 = t3 = 0; l = j; #pragma unroll (2) for (; l < K; l++) { t0 += b0[l] * a0[l * lda]; t1 += b1[l] * a0[l * lda]; t2 += b2[l] * a0[l * lda]; t3 += b3[l] * a0[l * lda]; } C[j] += alpha * t0; C[ldc + j] += alpha * t1; C[2 * ldc + j] += alpha * t2; C[3 * ldc + j] += alpha * t3; a0 += 1; } C += 4 * ldc; b0 += 4 * ldb; } for (i = 0; i < M - 3; i += 4) { a0 = (eml_32f *) A; b1 = b0 + ldb; b2 = b1 + ldb; b3 = b2 + ldb; for (j = 0; j < N; j++) { t0 = t1 = t2 = t3 = 0; l = j; #pragma unroll (2) for (; l < K; l++) { t0 += b0[l] * a0[l * lda]; t1 += b1[l] * a0[l * lda]; t2 += b2[l] * a0[l * lda]; t3 += b3[l] * a0[l * lda]; } C[j] += alpha * t0; C[ldc + j] += alpha * t1; C[2 * ldc + j] += alpha * t2; C[3 * ldc + j] += alpha * t3; a0 += 1; } C += 4 * ldc; b0 += 4 * ldb; } Применение оптимизации unroll к исходным текстам и последовательное обращение в память Цели: - Устранение зависимостей за счет применения оптимизации unroll к внешнему циклу - Улучшение загрузки конвейера за счет независимости операций в цикле после оптимизации unroll - Уменьшение количества запросов в память за счет многократного использования элементов матрицы - Эффективное использование поблочной загрузки и записи в память for (i = 0; i < M - 3; i += 4) { a0 = (eml_32f *) A; b1 = b0 + ldb; b2 = b1 + ldb; b3 = b2 + ldb; for (j = 0; j < N; j++) { t0 = t1 = t2 = t3 = 0; l = j; #pragma unroll (2) for (; l < K; l++) { t0 += b0[l] * a0[l * lda]; t1 += b1[l] * a0[l * lda]; t2 += b2[l] * a0[l * lda]; t3 += b3[l] * a0[l * lda]; } C[j] += alpha * t0; C[ldc + j] += alpha * t1; C[2 * ldc + j] += alpha * t2; C[3 * ldc + j] += alpha * t3; a0 += 1; } C += 4 * ldc; b0 += 4 * ldb; } for (i = 0; i < M - 3; i += 4) { a0 = (eml_32f *) A; b1 = b0 + ldb; b2 = b1 + ldb; b3 = b2 + ldb; for (j = 0; j < N; j++) { t0 = t1 = t2 = t3 = 0; l = j; #pragma unroll (2) for (; l < K; l++) { t0 += b0[l] * a0[l * lda]; t1 += b1[l] * a0[l * lda]; t2 += b2[l] * a0[l * lda]; t3 += b3[l] * a0[l * lda]; } C[j] += alpha * t0; C[ldc + j] += alpha * t1; C[2 * ldc + j] += alpha * t2; C[3 * ldc + j] += alpha * t3; a0 += 1; } C += 4 * ldc; b0 += 4 * ldb; } for (i = 0; i < M - 3; i += 4) { a0 = (eml_32f *) A; b1 = b0 + ldb; b2 = b1 + ldb; b3 = b2 + ldb; for (j = 0; j < N; j++) { t0 = t1 = t2 = t3 = 0; l = j; #pragma unroll (2) for (; l < K; l++) { t0 += b0[l] * a0[l * lda]; t1 += b1[l] * a0[l * lda]; t2 += b2[l] * a0[l * lda]; t3 += b3[l] * a0[l * lda]; } C[j] += alpha * t0; C[ldc + j] += alpha * t1; C[2 * ldc + j] += alpha * t2; C[3 * ldc + j] += alpha * t3; a0 += 1; } C += 4 * ldc; b0 += 4 * ldb; }

8 Вынос нелинейных индексов за пределы массивов - Сокращение накладных расходов на пересчет индексов for (l = 0; l < j; l++) { t0 += B[i * ldb + l] * A[j * lda + l]; t1 += B[(i + 1) * ldb + l] * A[j * lda + l]; t2 += B[(i + 2) * ldb + l] * A[j * lda + l]; t3 += B[(i + 3) * ldb + l] * A[j * lda + l]; } for (l = 0; l < j; l++) { t0 += b0[l] * a0[l]; t1 += b1[l] * a0[l]; t2 += b2[l] * a0[l]; t3 += b3[l] * a0[l]; } a0 = (eml_32f *) A; b0 = (eml_32f *) B; b1 = b0 + ldb; b2 = b1 + ldb; b3 = b2 + ldb; Итог: вместо 3х операций умножения выполняется 3 операции сложения

9 Особенности архитектуры Эльбрус-3М1, существенные при проведении оптимизаций Наличие векторных инструкций Конвейеризация циклов Асинхронный доступ к массивам (наличие APB) АЛУ: 2 кластера по 3 канала

10 Оптимизации функций умножения матриц для Эльбрус-3М1 Разбиение матриц на блоки – нет необходимости из-за наличия APB Применение оптимизации unroll к исходным текстам Устранение зависимостей между итерациями циклов Использование последовательного чтения/записи из/в память – нет необходимости из-за наличия APB Использование векторных инструкций Подбор состава операций для максимальной загрузки АЛУ

11 Использование векторных инструкций и подбор состава операций Цели: - Обработка нескольких элементов матрицы за такт - Максимальная загрузка АЛУ Проблемы: - Необходимость выравнивания данных в памяти по границе 64х битных слов … t10 = alpha * A[(i + 1) * lda + l]; t11 = alpha * A[(i + 1) * lda + l + 1]; t12 = alpha * A[(i + 1) * lda + l + 2]; t13 = alpha * A[(i + 1) * lda + l + 3]; … tmp10 = e3m_pshufh (*((eml_32s *) & t10), 0x44); tmp11 = e3m_pshufh (*((eml_32s *) & t11), 0x44); tmp12 = e3m_pshufh (*((eml_32s *) & t12), 0x44); tmp13 = e3m_pshufh (*((eml_32s *) & t13), 0x44); for (j = 0; j > 1; j++) { … sC1[j] = e3m_pfadds (sC1[j], e3m_pfmuls (tmp10, sB0[j])); sC1[j] = e3m_pfadds (sC1[j], e3m_pfmuls (tmp11, sB1[j])); sC1[j] = e3m_pfadds (sC1[j], e3m_pfmuls (tmp12, sB2[j])); sC1[j] = e3m_pfadds (sC1[j], e3m_pfmuls (tmp13, sB3[j])); … } … t10 = alpha * A[(i + 1) * lda + l]; t11 = alpha * A[(i + 1) * lda + l + 1]; t12 = alpha * A[(i + 1) * lda + l + 2]; t13 = alpha * A[(i + 1) * lda + l + 3]; … tmp10 = e3m_pshufh (*((eml_32s *) & t10), 0x44); tmp11 = e3m_pshufh (*((eml_32s *) & t11), 0x44); tmp12 = e3m_pshufh (*((eml_32s *) & t12), 0x44); tmp13 = e3m_pshufh (*((eml_32s *) & t13), 0x44); for (j = 0; j > 1; j++) { … sC1[j] = e3m_pfadds (sC1[j], e3m_pfmuls (tmp10, sB0[j])); sC1[j] = e3m_pfadds (sC1[j], e3m_pfmuls (tmp11, sB1[j])); sC1[j] = e3m_pfadds (sC1[j], e3m_pfmuls (tmp12, sB2[j])); sC1[j] = e3m_pfadds (sC1[j], e3m_pfmuls (tmp13, sB3[j])); … }

12 Результаты оптимизаций Тактов/эл-т Кол-во элементов матрицы Тактов/эл-т Кол-во элементов матрицы

13 Заключение Оптимизированы функции умножения матриц: умножение треугольных матриц различных конфигураций (слева – справа, верхнетреугольные – нижнетреугольные, транспонированные – не транспонированные) для различных типов данных (32F, 32FC, 64F, 64FC). Всего около 60 функций. Проведено функциональное тестирование Проведено измерение производительности Получен прирост производительности в среднем в 2-4 раза по сравнению с изначальными алгоритмами библиотеки BLAS, не использующими особенности ВК Эльбрус-90 микро и Эльбрус-3М1