Московский Физико-Технический Институт

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

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; } Цели: - Устранение зависимостейза счет применения оптимизации 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; } 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]; } a0 = (eml_32f *) A; b0 = (eml_32f *) B; b1 = b0 + ldb; b2 = b1 + ldb; b3 = b2 + ldb; 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]; } Итог:вместо 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 < N >> 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 < N >> 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”