Нерозв'язана проблема інформатики: Який алгоритм перемножування матриць є найшвидшим? |
Оскільки множення матриць є настільки центральною операцією в багатьох чисельних алгоритмах, у те, щоби зробити алгори́тми перемно́жування ма́триць ефективними, було вкладено чимало праці. Застосування множення матриць в обчислювальних задачах зустрічаються в багатьох областях, включно з [en] та розпізнаванням образів, і в, здавалося би, не пов'язаних задачах, таких як підрахунок шляхів графом. Було розроблено багато різних алгоритмів для перемножування матриць на різних типах апаратного забезпечення, включно з паралельними та розподіленими системами, де обчислювальну працю розподілювано декількома процесорами (можливо, через мережу).
Безпосереднє застосування математичного означення множення матриць дає алгоритм, що для перемножування двох матриць n × n займає часу порядку n3 (в нотації Ландау — Θ(n3)). Кращі асимптотичні межі часу, необхідного для перемножування матриць, були відомі від часів праці Страссена 1960 року, але досі не відомо, яким є оптимальний час (тобто, якою є складність цієї задачі).
Ітеративний алгоритм
За (означенням множення матриць), якщо C = AB для матриці A розміру n × m та матриці B розміру m × p, то C є матрицею розміру n × p з елементами
- .
З цього може бути побудовано простий алгоритм, який обходить циклами індекси i з 1 по n та j з 1 по p, обчислюючи наведене вище із застосуванням вкладеного циклу:
- Вхід: матриці A та B
- Нехай C буде матрицею відповідного розміру
- Для i з 1 по n:
- Для j з 1 по p:
- Нехай sum = 0
- Для k з 1 по m:
- Встановити sum ← sum + Aik × Bkj
- Встановити Cij ← sum
- Для j з 1 по p:
- Повернути C
Цей алгоритм займає час Θ(nmp) (в нотації Ландау). Поширеним спрощенням для цілей аналізу алгоритмів є вважати, що всі входи є квадратними матрицями розміру n × n, в разі чого час виконання становить Θ(n3), тобто, є кубічним.
Поведінка кешу
Ці три цикли в ітеративному перемножуванні матриць можливо довільно переставляти між собою без впливу на правильність чи асимптотичну тривалість виконання. Проте цей порядок може мати значний вплив на практичну продуктивність через [en] алгоритму та використання ним кешу; який порядок є кращим також залежить від того, чи зберігаються матриці в [en], постовпчиковому, чи суміші обох.
Зокрема, в ідеалізованому випадку (повністю асоціативного кешу), складеного з M байтів та b байтів на рядок кешу (тобто, M/b рядків кешу), наведений вище алгоритм є недооптимальним для A та B, що зберігаються в порядко́вому порядку. Коли n > M/b, кожна ітерація внутрішнього циклу (одночасного проходження рядком A та стовпчиком B) зазнає́ невлучання в кеш при отримуванні доступу до елементу B. Це означає, що цей алгоритм у найгіршому випадку зазнає Θ(n3) невлучань у кеш. Станом на 2010 рік швидкість пам'яті в порівнянні зі швидкістю процесорів є такою, що для матриць значного розміру в тривалості виконання домінують невлучання в кеш, а не фактичні обчислення.
Оптимальним варіантом ітеративного алгоритму для A та B, що зберігаються в порядку рядків, є [en] версія, в якій матрицю неявно ділять на квадратні плитки розміру √M на √M:
- Вхід: матриці A та B
- Нехай C буде матрицею відповідного розміру
- Обрати розмір плитки T = Θ(√M)
- Для I з 1 по n кроками по T:
- Для J з 1 по p кроками по T:
- Для K з 1 по m кроками по T:
- Перемножити AI:I+T, K:K+T та BK:K+T, J:J+T до CI:I+T, J:J+T, тобто:
- Для i з I по min(I + T, n):
- Для j з J по min(J + T, p):
- Нехай sum = 0
- Для k з K по min(K + T, m):
- Встановити sum ← sum + Aik × Bkj
- Встановити Cij ← Cij + sum
- Для j з J по min(J + T, p):
- Для K з 1 по m кроками по T:
- Для J з 1 по p кроками по T:
- Повернути C
В ідеалізованій моделі кешу цей алгоритм зазнає лише Θ(n3/b √M) невлучань у кеш; дільник b √M становить декілька порядків на сучасних машинах, тож у тривалості виконання домінують фактичні обчислення, а не невлучання в кеш.
Алгоритм «розділюй та володарюй»
Альтернативою до ітеративного алгоритму є алгоритм «розділюй та володарюй» для множення матриць. Він спирається на блокове розбиття
- ,
яке працює для всіх квадратних матриць, чиї розміри є степенями двійки, тобто, форми 2n × 2n для деякого n. Тепер добутком матриць є
що складається з восьми множень пар підматриць, з подальшим кроком додавання. Алгоритм «розділюй та володарюй» обчислює менші добутки рекурсивно, застосовуючи як основу скалярне множення c11 = a11b11.
Складність цього алгоритму як функції від n задається рекурентно як
- ;
- ,
із враховуванням восьми рекурсивних викликів на матрицях розміру n/2, та Θ(n2) для поелементного підсумовування чотирьох пар отримуваних в результаті матриць. Застосування майстер-методу для рекурсій «розділюй та володарюй» показує, що ця рекурсія має розв'язок Θ(n3), такий же, як й ітеративний алгоритм.
Неквадратні матриці
Варіант цього алгоритму, який працює для матриць довільних форм, і є швидшим на практиці, розбиває матриці на дві замість чотирьох підматриць наступним чином. Розбивання матриці тепер означає поділ її на дві частини рівного розміру, або якомога ближче до рівних розмірів у разі, якщо розміри є непарними.
- Входи: матриці A розміру n × m, B розміру m × p.
- Базовий випадок: якщо max(n, m, p) є нижчим за певний поріг, застосувати розмотану версію ітеративного алгоритму.
- Рекурсивні випадки:
- Якщо max(n, m, p) = n, розбити A горизонтально:
- Інакше, якщо max(n, m, p) = p, розбити B вертикально:
- Інакше, max(n, m, p) = m. Робити A вертикально та B горизонтально:
Поведінка кешу
Рівень невлучання в кеш в рекурсивного матричного множення є таким же, як і в [en] ітеративної версії, але, на відміну від того алгоритму, рекурсивний алгоритм є [en]: він не має параметру налаштування, необхідного для досягнення оптимальної продуктивності кешу, і він добре поводиться в [en] середовищі, де розміри кешу є фактично динамічними через те, що інші процеси займають його простір. (Простий ітеративний алгоритм також є буферо-незалежним, але набагато повільнішим на практиці, якщо компонування матриць не пристосовано до алгоритму.)
Кількість невлучань у кеш, що зазнає́ цей алгоритм на машині з M рядками ідеально кешу розміру b байтів кожен, обмежено
Підкубічні алгоритми
Існують алгоритми, що забезпечують кращу тривалість виконання, ніж прямолінійні. Першим було відкрито алгоритм Штрассена, винайдений Фолькером Штрассеном 1969 року, який часто називають «швидким множенням матриць» (англ. "fast matrix multiplication"). Він ґрунтується на способі множення двох матриць 2 × 2, який вимагає лише 7 множень (замість звичайних 8), ціною декількох додаткових операції додавання та віднімання. Його рекурсивне застосування дає алгоритм з витратами на множення . Алгоритм Штрассена є складнішим та має нижчу чисельну стійкість у порівнянні з наївним алгоритмом, але він є швидшим у випадках, коли n > 100 або близько того, і зустрічається в декількох бібліотеках, таких як BLAS. Він є дуже корисним для великих матриць над точними областями, такими як скінченні поля, де чисельна стійкість не є проблемою.
Поточний алгоритм O(nk) з найнижчим відомим степенем k є узагальненнями алгоритму Копперсміта — Вінограда від Франсуа ле Ґалля, яке має асимптотичну складність O(n2.3728639). Алгоритм ле Ґалля та алгоритм Копперсміта — Вінограда, на якому він ґрунтується, є подібними до алгоритму Штрассена: винайдено спосіб множення двох матриць k × k менше ніж k3 множеннями, і цю методику застосовують рекурсивно. Проте сталий коефіцієнт, прихований нотацією Ландау, є настільки великим, що ці алгоритми є доцільними лише для матриць, що є завеликими для обробки на сучасних комп'ютерах.
Оскільки будь-який алгоритм для перемножування двох матриць n × n має обробити всі 2n2 елементів, існує асимптотична нижня межа в Ω(n2) операцій. Рац довів нижню межу в Ω(n2 log(n)) для арифметичних схем з обмеженими коефіцієнтами над дійсними або комплексними числами.
Кон та ін. помістили такі методи, як алгоритми Штрассена та Копперсміта — Вінограда, до зовсім відмінного контексту теорії груп, використавши трійки підмножин скінченних груп, які задовольняють властивості неперетинності, що називають [en] (ВПД, англ. triple product property, TPP). Вони показали, що якщо сімейства [en] абелевих груп з симетричними групами втілюють сімейства трійок підмножин зі спільною версією ВПД, то існують алгоритми перемножування матриць із по суті квадратичною складністю. Більшість дослідників вважають, що це так і є. Проте Алон, Шпілка та [en] нещодавно показали, що деякі з цих гіпотез, що передбачають швидке множення матриць, є несумісними з іншою правдоподібною гіпотезою, [en].
[en] є простим алгоритмом Монте-Карло, який для заданих матриць A, B та C перевіряє за час Θ(n2), чи AB = C.
Паралельні та розподілені алгоритми
Розпаралелювання зі спільною пам'яттю
Окреслений вище алгоритм «розділюй та володарюй» можливо розпаралелити для багатопроцесорності зі спільною пам'яттю двома способами. Вони ґрунтуються на тім факті, що вісім рекурсивних перемножувань матриць у
можливо виконувати незалежно одне від одного, як і чотири підсумовування (хоча алгоритмові потрібно «об'єднати» перемножування перед виконанням додавань). Використовуючи повну паралельність задачі, отримують алгоритм, який може бути виражено псевдокодом стилю [en]:
Процедура перемножити(C, A, B):
- Базовий випадок: якщо n = 1, встановити c11 ← a11 × b11 (або перемножити маленьку блокову матрицю).
- Інакше, виділити простір для нової матриці T форми n × n, а тоді:
- Розбити A на A11, A12, A21, A22.
- Розбити B на B11, B12, B21, B22.
- Розбити C на C11, C12, C21, C22.
- Розбити T на T11, T12, T21, T22.
- Паралельне виконання:
- Відгалузити перемножити(C11, A11, B11).
- Відгалузити перемножити(C12, A11, B12).
- Відгалузити перемножити(C21, A21, B11).
- Відгалузити перемножити(C22, A21, B12).
- Відгалузити перемножити(T11, A12, B21).
- Відгалузити перемножити(T12, A12, B22).
- Відгалузити перемножити(T21, A22, B21).
- Відгалузити перемножити(T22, A22, B22).
- Об'єднати (дочекатися завершення паралельних відгалужень).
- додати(C, T).
- Звільнити T.
Процедура додати(C, T) додає T до C, поелементно:
- Базовий випадок: якщо n = 1, встановити c11 ← c11 + t11 (або виконати короткий цикл, можливо, розмотаний).
- Інакше:
- Розбити C на C11, C12, C21, C22.
- Розбити T на T11, T12, T21, T22.
- Паралельно:
- Відгалузити додати(C11, T11).
- Відгалузити додати(C12, T12).
- Відгалузити додати(C21, T21).
- Відгалузити додати(C22, T22).
- Об'єднати.
Тут відгалузити є ключовим словом, яке сигналізує, що обчислення може бути виконувано паралельно до решти виклику функції, тоді як об'єднати чекає на завершення всіх попередньо «відгалужених» обчислень. Розбити досягає своєї мети лише маніпулюванням вказівниками.
Цей алгоритм має критичною довжиною шляху Θ(log2 n) кроків, що означає, що йому потрібно стільки часу на ідеальній машині з нескінченним числом процесорів. Отже, на будь-якому реальному комп'ютері він має максимальне можливе прискорення Θ(n3/log2 n). Цей алгоритм не є практичним через витрати на передавання, властиві переміщуванню даних до та з тимчасової матриці T, але практичніший варіант досягає прискорення Θ(n2) без використання тимчасової матриці.
Алгоритми з униканням передавання, та розподілені алгоритми
На сучасних архітектурах з ієрархічною пам'яттю вартість завантажування та зберігання елементів матриць входу має схильність переважати над вартістю арифметики. На одній машині це — кількість даних, передаваних між оперативною пам'яттю та кешем, тоді як на багатовузловій машині з розподіленою пам'яттю це — кількість, що передається між вузлами. В кожному з випадків її називають пропускною здатністю обміну (англ. communication bandwidth). Наївний алгоритм із використанням трьох вкладених циклів використовує пропускну здатністю обміну Ω(n3).
[en], відомий також як двовимірний алгоритм (англ. 2D algorithm), є [en], який розбиває кожну з матриць входу на блокову матрицю, чиї елементи є підматрицями розміру √M/3 на √M/3, де M є розміром швидкої пам'яті. Потім застосовують наївний алгоритм над блоковими матрицями, обчислюючи добутки підматриць цілком у швидкій пам'яті. Це знижує пропускну здатність обміну до O(n3/√M), яка є асимптотично оптимальною (для алгоритмів, що виконують обчислення Ω(n3)).
У розподіленій постановці з p процесорами, розставленими у двовимірній сітці √p на √p, кожному з процесорів можливо призначувати по одній з підматриць результату, і добуток буде обчислювано з передаванням кожним з процесорів O(n2/√p) слів, що є асимптотично оптимальним, виходячи з того, що кожен з вузлів зберігає щонайменше O(n2/p) елементів. Це може бути вдосконалено тривимірним алгоритмом (англ. 3D algorithm), який впорядковує процесори у тривимірну кубічну сітку, призначуючи кожен добуток двох підматриць входу одному процесорові. Підматриці результату потім породжують виконанням зведення над кожним з рядків. Цей алгоритм передає O(n2/p2/3) слів на процесор, що є асимптотично оптимальним. Проте, це вимагає повторювання кожного з елементів матриць входу p1/3 разів, і відтак вимагає в p1/3 разів більше пам'яті, ніж необхідно для зберігання входів. Для подальшого зниження тривалості виконання цей алгоритм може бути поєднувано з алгоритмом Штрассена. «2,5-вимірні» алгоритми (англ. "2.5D" algorithms) забезпечують безперервний компроміс між використанням пам'яті та пропускною здатністю передавання. На сучасних розподілених обчислювальних середовищах, таких як MapReduce, було розроблено спеціалізовані алгоритми перемножування.
Алгоритми для сіток
Існує низка алгоритмів для перемножування на сітках. Для перемножування двох n×n на стандартній двовимірній сітці із застосуванням двовимірного [en] обчислюють перемножування в 3n-2 кроків, хоча це знижується до половини цього числа для повторюваних обчислень. Стандартний масив є неефективним, оскільки дані з двох матриць не надходять одночасно, і його має бути доповнювано нулями.
Результат є ще швидшим на двошаровій перехрещуваній сітці, де потрібно лише 2n-1 кроків. Продуктивність додатково покращується для повторюваних обчислень, ведучи до стовідсоткової ефективності. Перехрещуваний сітковий масив можна розглядати як особливий випадок не планарної (тобто, багатошарової) оброблювальної структури.
Див. також
- [en]
- [en]
- [en]
- [en]
Примітки
- (2008). Sorting and Searching. The Algorithm Design Manual. Springer. с. 45–46, 401—403. doi:10.1007/978-1-84800-070-4_4. ISBN . (англ.)
- Томас Кормен; Чарльз Лейзерсон, Рональд Рівест, Кліфорд Стайн (2009) [1990]. 28 Дії над матрицями. Вступ до алгоритмів (вид. 3rd). MIT Press і McGraw-Hill. ISBN .
- Amarasinghe, Saman; Leiserson, Charles (2010). 6.172 Performance Engineering of Software Systems, Lecture 8. MIT OpenCourseWare. Massachusetts Institute of Technology. оригіналу за 7 жовтня 2019. Процитовано 27 січня 2015. (англ.)
- Lam, Monica S.; Rothberg, Edward E.; Wolf, Michael E. (1991). The Cache Performance and Optimizations of Blocked Algorithms. Int'l Conf. on Architectural Support for Programming Languages and Operating Systems (ASPLOS). (англ.)
- (1999). (PDF) (Master's). MIT. Архів оригіналу (PDF) за 24 листопада 2019. Процитовано 15 грудня 2019. (англ.)
- Miller, Webb (1975), Computational complexity and numerical stability, SIAM News, 4: 97—107, CiteSeerX 10.1.1.148.9947 (англ.)
- Press, William H.; Flannery, Brian P.; ; Vetterling, William T. (2007). Numerical Recipes: The Art of Scientific Computing (вид. 3rd). Cambridge University Press. ISBN . (англ.)
- Le Gall, François (2014), Powers of tensors and fast matrix multiplication, Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation (ISSAC 2014), arXiv:1401.7714, Bibcode:2014arXiv1401.7714L (англ.). Первинний алгоритм, представлений Доном Копперсмітом та [en] 1990 року, має асимптотичну складність O(n2.376). Його було поліпшено 2013 року до O(n2.3729) [en], що дало тривалість лише трошки гіршу, ніж у вдосконалення ле Ґалля: . (PDF). Архів оригіналу (PDF) за 8 жовтня 2013. Процитовано 15 грудня 2019. (англ.)
- Iliopoulos, Costas S. (1989), (PDF), SIAM Journal on Computing, 18 (4): 658—669, CiteSeerX 10.1.1.531.9309, doi:10.1137/0218045, MR 1004789, архів оригіналу (PDF) за 5 березня 2014, процитовано 15 грудня 2019,
The Coppersmith–Winograd algorithm is not practical, due to the very large hidden constant in the upper bound on the number of multiplications required.
(англ.) - Robinson, Sara (2005). (PDF). SIAM News. 38 (9). Архів оригіналу (PDF) за 31 березня 2019. Процитовано 15 грудня 2019. (англ.)
- [en]. On the complexity of matrix product. In Proceedings of the thirty-fourth annual ACM symposium on Theory of computing. ACM Press, 2002. DOI:10.1145/509907.509932. (англ.)
- Henry Cohn, [en], [en], and Chris Umans. Group-theoretic Algorithms for Matrix Multiplication. arXiv:math.GR/0511460. Proceedings of the 46th Annual Symposium on Foundations of Computer Science, 23–25 October 2005, Pittsburgh, PA, IEEE Computer Society, pp. 379–388. (англ.)
- Henry Cohn, Chris Umans. A Group-theoretic Approach to Fast Matrix Multiplication. arXiv:math.GR/0307321. Proceedings of the 44th Annual IEEE Symposium on Foundations of Computer Science, 11–14 October 2003, Cambridge, MA, IEEE Computer Society, pp. 438–449. (англ.)
- Alon, Shpilka, Umans, On Sunflowers and Matrix Multiplication [ 9 грудня 2016 у Wayback Machine.] (англ.)
- Randall, Keith H. (1998). (PDF) (Ph.D.). Massachusetts Institute of Technology. с. 54—57. Архів оригіналу (PDF) за 6 листопада 2020. Процитовано 15 грудня 2019. (англ.)
- Lynn Elliot Cannon, A cellular computer to implement the Kalman Filter Algorithm, Technical report, Ph.D. Thesis, Montana State University, 14 July 1969. (англ.)
- Hong, J. W.; Kung, H. T. (1981). (PDF). STOC '81: Proceedings of the Thirteenth Annual ACM Symposium on Theory of Computing: 326—333. Архів оригіналу (PDF) за 15 грудня 2019. Процитовано 15 грудня 2019. (англ.)
- Irony, Dror; Toledo, Sivan; Tiskin, Alexander (September 2004). Communication lower bounds for distributed-memory matrix multiplication. J. Parallel Distrib. Comput. 64 (9): 1017—1026. CiteSeerX 10.1.1.20.7034. doi:10.1016/j.jpdc.2004.03.021. (англ.)
- Agarwal, R.C.; Balle, S. M.; Gustavson, F. G.; Joshi, M.; Palkar, P. (September 1995). A three-dimensional approach to parallel matrix multiplication. IBM J. Res. Dev. 39 (5): 575—582. CiteSeerX 10.1.1.44.3404. doi:10.1147/rd.395.0575. (англ.)
- Solomonik, Edgar; Demmel, James (2011). Communication-optimal parallel 2.5D matrix multiplication and LU factorization algorithms. Proceedings of the 17th International Conference on Parallel Processing. Part II: 90—109. (англ.)
- Bosagh Zadeh, Reza; Carlsson, Gunnar. (PDF). Архів оригіналу (PDF) за 14 липня 2014. Процитовано 12 липня 2014. (англ.)
- Bae, S.E., T.-W. Shinn, T. Takaoka (2014) A faster parallel algorithm for matrix multiplication on a mesh array [ 18 травня 2021 у Wayback Machine.]. Procedia Computer Science 29: 2230-2240 (англ.)
- Kak, S. (1988) A two-layered mesh array for matrix multiplication. Parallel Computing 6: 383-385 (англ.)
- Kak, S. (2014) Efficiency of matrix multiplication on the cross-wired mesh array. https://arxiv.org/abs/1411.3273 [ 23 березня 2019 у Wayback Machine.] (англ.)
- Kak, S. (1988) Multilayered array computing. Information Sciences 45: 347-365 (англ.)
Література
- Buttari, Alfredo; Langou, Julien; Kurzak, Jakub; (2009). A class of parallel tiled linear algebra algorithms for multicore architectures. Parallel Computing. 35: 38—53. arXiv:0709.1272. doi:10.1016/j.parco.2008.10.002. (англ.)
- Goto, Kazushige; van de Geijn, Robert A. (2008). Anatomy of high-performance matrix multiplication. ACM Transactions on Mathematical Software. 34 (3): 1—25. CiteSeerX 10.1.1.140.3583. doi:10.1145/1356052.1356053. (англ.)
- Van Zee, Field G.; van de Geijn, Robert A. (2015). BLIS: A Framework for Rapidly Instantiating BLAS Functionality. ACM Transactions on Mathematical Software. 41 (3): 1—33. doi:10.1145/2764454. (англ.)
- How To Optimize GEMM [ 15 вересня 2020 у Wayback Machine.] (англ.)
Вікіпедія, Українська, Україна, книга, книги, бібліотека, стаття, читати, завантажити, безкоштовно, безкоштовно завантажити, mp3, відео, mp4, 3gp, jpg, jpeg, gif, png, малюнок, музика, пісня, фільм, книга, гра, ігри, мобільний, телефон, android, ios, apple, мобільний телефон, samsung, iphone, xiomi, xiaomi, redmi, honor, oppo, nokia, sonya, mi, ПК, web, Інтернет
Nerozv yazana problema informatiki Yakij algoritm peremnozhuvannya matric ye najshvidshim Oskilki mnozhennya matric ye nastilki centralnoyu operaciyeyu v bagatoh chiselnih algoritmah u te shobi zrobiti algori tmi peremno zhuvannya ma tric efektivnimi bulo vkladeno chimalo praci Zastosuvannya mnozhennya matric v obchislyuvalnih zadachah zustrichayutsya v bagatoh oblastyah vklyuchno z en ta rozpiznavannyam obraziv i v zdavalosya bi ne pov yazanih zadachah takih yak pidrahunok shlyahiv grafom Bulo rozrobleno bagato riznih algoritmiv dlya peremnozhuvannya matric na riznih tipah aparatnogo zabezpechennya vklyuchno z paralelnimi ta rozpodilenimi sistemami de obchislyuvalnu pracyu rozpodilyuvano dekilkoma procesorami mozhlivo cherez merezhu Bezposerednye zastosuvannya matematichnogo oznachennya mnozhennya matric daye algoritm sho dlya peremnozhuvannya dvoh matric n n zajmaye chasu poryadku n3 v notaciyi Landau 8 n3 Krashi asimptotichni mezhi chasu neobhidnogo dlya peremnozhuvannya matric buli vidomi vid chasiv praci Strassena 1960 roku ale dosi ne vidomo yakim ye optimalnij chas tobto yakoyu ye skladnist ciyeyi zadachi Iterativnij algoritmZa oznachennyam mnozhennya matric yaksho C AB dlya matrici A rozmiru n m ta matrici B rozmiru m p to C ye matriceyu rozmiru n p z elementami c i j k 1 m a i k b k j displaystyle c ij sum k 1 m a ik b kj Z cogo mozhe buti pobudovano prostij algoritm yakij obhodit ciklami indeksi i z 1 po n ta j z 1 po p obchislyuyuchi navedene vishe iz zastosuvannyam vkladenogo ciklu Vhid matrici A ta B Nehaj C bude matriceyu vidpovidnogo rozmiru Dlya i z 1 po n Dlya j z 1 po p Nehaj sum 0 Dlya k z 1 po m Vstanoviti sum sum Aik Bkj Vstanoviti Cij sum Povernuti C Cej algoritm zajmaye chas 8 nmp v notaciyi Landau Poshirenim sproshennyam dlya cilej analizu algoritmiv ye vvazhati sho vsi vhodi ye kvadratnimi matricyami rozmiru n n v razi chogo chas vikonannya stanovit 8 n3 tobto ye kubichnim Povedinka keshu Ilyustraciya poryadko vogo ta postovpchikovogo poryadku Ci tri cikli v iterativnomu peremnozhuvanni matric mozhlivo dovilno perestavlyati mizh soboyu bez vplivu na pravilnist chi asimptotichnu trivalist vikonannya Prote cej poryadok mozhe mati znachnij vpliv na praktichnu produktivnist cherez en algoritmu ta vikoristannya nim keshu yakij poryadok ye krashim takozh zalezhit vid togo chi zberigayutsya matrici v en postovpchikovomu chi sumishi oboh Zokrema v idealizovanomu vipadku povnistyu asociativnogo keshu skladenogo z M bajtiv ta b bajtiv na ryadok keshu tobto M b ryadkiv keshu navedenij vishe algoritm ye nedooptimalnim dlya A ta B sho zberigayutsya v poryadko vomu poryadku Koli n gt M b kozhna iteraciya vnutrishnogo ciklu odnochasnogo prohodzhennya ryadkom A ta stovpchikom B zaznaye nevluchannya v kesh pri otrimuvanni dostupu do elementu B Ce oznachaye sho cej algoritm u najgirshomu vipadku zaznaye 8 n3 nevluchan u kesh Stanom na 2010 rik shvidkist pam yati v porivnyanni zi shvidkistyu procesoriv ye takoyu sho dlya matric znachnogo rozmiru v trivalosti vikonannya dominuyut nevluchannya v kesh a ne faktichni obchislennya Optimalnim variantom iterativnogo algoritmu dlya A ta B sho zberigayutsya v poryadku ryadkiv ye en versiya v yakij matricyu neyavno dilyat na kvadratni plitki rozmiru M na M Vhid matrici A ta B Nehaj C bude matriceyu vidpovidnogo rozmiru Obrati rozmir plitki T 8 M Dlya I z 1 po n krokami po T Dlya J z 1 po p krokami po T Dlya K z 1 po m krokami po T Peremnozhiti AI I T K K T ta BK K T J J T do CI I T J J T tobto Dlya i z I po min I T n Dlya j z J po min J T p Nehaj sum 0 Dlya k z K po min K T m Vstanoviti sum sum Aik Bkj Vstanoviti Cij Cij sum Povernuti C V idealizovanij modeli keshu cej algoritm zaznaye lishe 8 n3 b M nevluchan u kesh dilnik b M stanovit dekilka poryadkiv na suchasnih mashinah tozh u trivalosti vikonannya dominuyut faktichni obchislennya a ne nevluchannya v kesh Algoritm rozdilyuj ta volodaryuj Alternativoyu do iterativnogo algoritmu ye algoritm rozdilyuj ta volodaryuj dlya mnozhennya matric Vin spirayetsya na blokove rozbittya C C 11 C 12 C 21 C 22 A A 11 A 12 A 21 A 22 B B 11 B 12 B 21 B 22 displaystyle C begin pmatrix C 11 amp C 12 C 21 amp C 22 end pmatrix A begin pmatrix A 11 amp A 12 A 21 amp A 22 end pmatrix B begin pmatrix B 11 amp B 12 B 21 amp B 22 end pmatrix yake pracyuye dlya vsih kvadratnih matric chiyi rozmiri ye stepenyami dvijki tobto formi 2n 2n dlya deyakogo n Teper dobutkom matric ye C 11 C 12 C 21 C 22 A 11 A 12 A 21 A 22 B 11 B 12 B 21 B 22 A 11 B 11 A 12 B 21 A 11 B 12 A 12 B 22 A 21 B 11 A 22 B 21 A 21 B 12 A 22 B 22 displaystyle begin pmatrix C 11 amp C 12 C 21 amp C 22 end pmatrix begin pmatrix A 11 amp A 12 A 21 amp A 22 end pmatrix begin pmatrix B 11 amp B 12 B 21 amp B 22 end pmatrix begin pmatrix A 11 B 11 A 12 B 21 amp A 11 B 12 A 12 B 22 A 21 B 11 A 22 B 21 amp A 21 B 12 A 22 B 22 end pmatrix sho skladayetsya z vosmi mnozhen par pidmatric z podalshim krokom dodavannya Algoritm rozdilyuj ta volodaryuj obchislyuye menshi dobutki rekursivno zastosovuyuchi yak osnovu skalyarne mnozhennya c11 a11b11 Skladnist cogo algoritmu yak funkciyi vid n zadayetsya rekurentno yak T 1 8 1 displaystyle T 1 Theta 1 T n 8 T n 2 8 n 2 displaystyle T n 8T n 2 Theta n 2 iz vrahovuvannyam vosmi rekursivnih viklikiv na matricyah rozmiru n 2 ta 8 n2 dlya poelementnogo pidsumovuvannya chotiroh par otrimuvanih v rezultati matric Zastosuvannya majster metodu dlya rekursij rozdilyuj ta volodaryuj pokazuye sho cya rekursiya maye rozv yazok 8 n3 takij zhe yak j iterativnij algoritm Nekvadratni matrici Variant cogo algoritmu yakij pracyuye dlya matric dovilnih form i ye shvidshim na praktici rozbivaye matrici na dvi zamist chotiroh pidmatric nastupnim chinom Rozbivannya matrici teper oznachaye podil yiyi na dvi chastini rivnogo rozmiru abo yakomoga blizhche do rivnih rozmiriv u razi yaksho rozmiri ye neparnimi Vhodi matrici A rozmiru n m B rozmiru m p Bazovij vipadok yaksho max n m p ye nizhchim za pevnij porig zastosuvati rozmotanu versiyu iterativnogo algoritmu Rekursivni vipadki Yaksho max n m p n rozbiti A gorizontalno C A 1 A 2 B A 1 B A 2 B displaystyle C begin pmatrix A 1 A 2 end pmatrix B begin pmatrix A 1 B A 2 B end pmatrix dd Inakshe yaksho max n m p p rozbiti B vertikalno C A B 1 B 2 A B 1 A B 2 displaystyle C A begin pmatrix B 1 amp B 2 end pmatrix begin pmatrix AB 1 amp AB 2 end pmatrix dd Inakshe max n m p m Robiti A vertikalno ta B gorizontalno C A 1 A 2 B 1 B 2 A 1 B 1 A 2 B 2 displaystyle C begin pmatrix A 1 amp A 2 end pmatrix begin pmatrix B 1 B 2 end pmatrix A 1 B 1 A 2 B 2 dd Povedinka keshu Riven nevluchannya v kesh v rekursivnogo matrichnogo mnozhennya ye takim zhe yak i v en iterativnoyi versiyi ale na vidminu vid togo algoritmu rekursivnij algoritm ye en vin ne maye parametru nalashtuvannya neobhidnogo dlya dosyagnennya optimalnoyi produktivnosti keshu i vin dobre povoditsya v en seredovishi de rozmiri keshu ye faktichno dinamichnimi cherez te sho inshi procesi zajmayut jogo prostir Prostij iterativnij algoritm takozh ye bufero nezalezhnim ale nabagato povilnishim na praktici yaksho komponuvannya matric ne pristosovano do algoritmu Kilkist nevluchan u kesh sho zaznaye cej algoritm na mashini z M ryadkami idealno keshu rozmiru b bajtiv kozhen obmezheno 13 8 m n p m n n p m p b m n p b M displaystyle Theta left m n p frac mn np mp b frac mnp b sqrt M right Pidkubichni algoritmiNajnizhche w take sho vidomo sho mnozhennya matric ye v O nw vidobrazhene za chasom Isnuyut algoritmi sho zabezpechuyut krashu trivalist vikonannya nizh pryamolinijni Pershim bulo vidkrito algoritm Shtrassena vinajdenij Folkerom Shtrassenom 1969 roku yakij chasto nazivayut shvidkim mnozhennyam matric angl fast matrix multiplication Vin gruntuyetsya na sposobi mnozhennya dvoh matric 2 2 yakij vimagaye lishe 7 mnozhen zamist zvichajnih 8 cinoyu dekilkoh dodatkovih operaciyi dodavannya ta vidnimannya Jogo rekursivne zastosuvannya daye algoritm z vitratami na mnozhennya O n log 2 7 O n 2 807 displaystyle O n log 2 7 approx O n 2 807 Algoritm Shtrassena ye skladnishim ta maye nizhchu chiselnu stijkist u porivnyanni z nayivnim algoritmom ale vin ye shvidshim u vipadkah koli n gt 100 abo blizko togo i zustrichayetsya v dekilkoh bibliotekah takih yak BLAS Vin ye duzhe korisnim dlya velikih matric nad tochnimi oblastyami takimi yak skinchenni polya de chiselna stijkist ne ye problemoyu Potochnij algoritm O nk z najnizhchim vidomim stepenem k ye uzagalnennyami algoritmu Koppersmita Vinograda vid Fransua le Gallya yake maye asimptotichnu skladnist O n2 3728639 Algoritm le Gallya ta algoritm Koppersmita Vinograda na yakomu vin gruntuyetsya ye podibnimi do algoritmu Shtrassena vinajdeno sposib mnozhennya dvoh matric k k menshe nizh k3 mnozhennyami i cyu metodiku zastosovuyut rekursivno Prote stalij koeficiyent prihovanij notaciyeyu Landau ye nastilki velikim sho ci algoritmi ye docilnimi lishe dlya matric sho ye zavelikimi dlya obrobki na suchasnih komp yuterah Oskilki bud yakij algoritm dlya peremnozhuvannya dvoh matric n n maye obrobiti vsi 2n2 elementiv isnuye asimptotichna nizhnya mezha v W n2 operacij Rac doviv nizhnyu mezhu v W n2 log n dlya arifmetichnih shem z obmezhenimi koeficiyentami nad dijsnimi abo kompleksnimi chislami Kon ta in pomistili taki metodi yak algoritmi Shtrassena ta Koppersmita Vinograda do zovsim vidminnogo kontekstu teoriyi grup vikoristavshi trijki pidmnozhin skinchennih grup yaki zadovolnyayut vlastivosti neperetinnosti sho nazivayut en VPD angl triple product property TPP Voni pokazali sho yaksho simejstva en abelevih grup z simetrichnimi grupami vtilyuyut simejstva trijok pidmnozhin zi spilnoyu versiyeyu VPD to isnuyut algoritmi peremnozhuvannya matric iz po suti kvadratichnoyu skladnistyu Bilshist doslidnikiv vvazhayut sho ce tak i ye Prote Alon Shpilka ta en neshodavno pokazali sho deyaki z cih gipotez sho peredbachayut shvidke mnozhennya matric ye nesumisnimi z inshoyu pravdopodibnoyu gipotezoyu en en ye prostim algoritmom Monte Karlo yakij dlya zadanih matric A B ta C pereviryaye za chas 8 n2 chi AB C Paralelni ta rozpodileni algoritmiRozparalelyuvannya zi spilnoyu pam yattyu Okreslenij vishe algoritm rozdilyuj ta volodaryuj mozhlivo rozparaleliti dlya bagatoprocesornosti zi spilnoyu pam yattyu dvoma sposobami Voni gruntuyutsya na tim fakti sho visim rekursivnih peremnozhuvan matric u A 11 B 11 A 12 B 21 A 11 B 12 A 12 B 22 A 21 B 11 A 22 B 21 A 21 B 12 A 22 B 22 displaystyle begin pmatrix A 11 B 11 A 12 B 21 amp A 11 B 12 A 12 B 22 A 21 B 11 A 22 B 21 amp A 21 B 12 A 22 B 22 end pmatrix mozhlivo vikonuvati nezalezhno odne vid odnogo yak i chotiri pidsumovuvannya hocha algoritmovi potribno ob yednati peremnozhuvannya pered vikonannyam dodavan Vikoristovuyuchi povnu paralelnist zadachi otrimuyut algoritm yakij mozhe buti virazheno psevdokodom stilyu en Procedura peremnozhiti C A B Bazovij vipadok yaksho n 1 vstanoviti c11 a11 b11 abo peremnozhiti malenku blokovu matricyu Inakshe vidiliti prostir dlya novoyi matrici T formi n n a todi Rozbiti A na A11 A12 A21 A22 Rozbiti B na B11 B12 B21 B22 Rozbiti C na C11 C12 C21 C22 Rozbiti T na T11 T12 T21 T22 Paralelne vikonannya Vidgaluziti peremnozhiti C11 A11 B11 Vidgaluziti peremnozhiti C12 A11 B12 Vidgaluziti peremnozhiti C21 A21 B11 Vidgaluziti peremnozhiti C22 A21 B12 Vidgaluziti peremnozhiti T11 A12 B21 Vidgaluziti peremnozhiti T12 A12 B22 Vidgaluziti peremnozhiti T21 A22 B21 Vidgaluziti peremnozhiti T22 A22 B22 Ob yednati dochekatisya zavershennya paralelnih vidgaluzhen dodati C T Zvilniti T Procedura dodati C T dodaye T do C poelementno Bazovij vipadok yaksho n 1 vstanoviti c11 c11 t11 abo vikonati korotkij cikl mozhlivo rozmotanij Inakshe Rozbiti C na C11 C12 C21 C22 Rozbiti T na T11 T12 T21 T22 Paralelno Vidgaluziti dodati C11 T11 Vidgaluziti dodati C12 T12 Vidgaluziti dodati C21 T21 Vidgaluziti dodati C22 T22 Ob yednati Tut vidgaluziti ye klyuchovim slovom yake signalizuye sho obchislennya mozhe buti vikonuvano paralelno do reshti vikliku funkciyi todi yak ob yednati chekaye na zavershennya vsih poperedno vidgaluzhenih obchislen Rozbiti dosyagaye svoyeyi meti lishe manipulyuvannyam vkazivnikami Cej algoritm maye kritichnoyu dovzhinoyu shlyahu 8 log2 n krokiv sho oznachaye sho jomu potribno stilki chasu na idealnij mashini z neskinchennim chislom procesoriv Otzhe na bud yakomu realnomu komp yuteri vin maye maksimalne mozhlive priskorennya 8 n3 log2 n Cej algoritm ne ye praktichnim cherez vitrati na peredavannya vlastivi peremishuvannyu danih do ta z timchasovoyi matrici T ale praktichnishij variant dosyagaye priskorennya 8 n2 bez vikoristannya timchasovoyi matrici Blokove peremnozhuvannya matric U dvovimirnomu algoritmi kozhen iz procesoriv ye vidpovidalnim za odnu pidmatricyu C U trivimirnomu algoritmi kozhnu z par peremnozhuvanih pidmatric z A ta B priznachayut odnomu procesorovi Algoritmi z unikannyam peredavannya ta rozpodileni algoritmi Na suchasnih arhitekturah z iyerarhichnoyu pam yattyu vartist zavantazhuvannya ta zberigannya elementiv matric vhodu maye shilnist perevazhati nad vartistyu arifmetiki Na odnij mashini ce kilkist danih peredavanih mizh operativnoyu pam yattyu ta keshem todi yak na bagatovuzlovij mashini z rozpodilenoyu pam yattyu ce kilkist sho peredayetsya mizh vuzlami V kozhnomu z vipadkiv yiyi nazivayut propusknoyu zdatnistyu obminu angl communication bandwidth Nayivnij algoritm iz vikoristannyam troh vkladenih cikliv vikoristovuye propusknu zdatnistyu obminu W n3 en vidomij takozh yak dvovimirnij algoritm angl 2D algorithm ye en yakij rozbivaye kozhnu z matric vhodu na blokovu matricyu chiyi elementi ye pidmatricyami rozmiru M 3 na M 3 de M ye rozmirom shvidkoyi pam yati Potim zastosovuyut nayivnij algoritm nad blokovimi matricyami obchislyuyuchi dobutki pidmatric cilkom u shvidkij pam yati Ce znizhuye propusknu zdatnist obminu do O n3 M yaka ye asimptotichno optimalnoyu dlya algoritmiv sho vikonuyut obchislennya W n3 U rozpodilenij postanovci z p procesorami rozstavlenimi u dvovimirnij sitci p na p kozhnomu z procesoriv mozhlivo priznachuvati po odnij z pidmatric rezultatu i dobutok bude obchislyuvano z peredavannyam kozhnim z procesoriv O n2 p sliv sho ye asimptotichno optimalnim vihodyachi z togo sho kozhen z vuzliv zberigaye shonajmenshe O n2 p elementiv Ce mozhe buti vdoskonaleno trivimirnim algoritmom angl 3D algorithm yakij vporyadkovuye procesori u trivimirnu kubichnu sitku priznachuyuchi kozhen dobutok dvoh pidmatric vhodu odnomu procesorovi Pidmatrici rezultatu potim porodzhuyut vikonannyam zvedennya nad kozhnim z ryadkiv Cej algoritm peredaye O n2 p2 3 sliv na procesor sho ye asimptotichno optimalnim Prote ce vimagaye povtoryuvannya kozhnogo z elementiv matric vhodu p1 3 raziv i vidtak vimagaye v p1 3 raziv bilshe pam yati nizh neobhidno dlya zberigannya vhodiv Dlya podalshogo znizhennya trivalosti vikonannya cej algoritm mozhe buti poyednuvano z algoritmom Shtrassena 2 5 vimirni algoritmi angl 2 5D algorithms zabezpechuyut bezperervnij kompromis mizh vikoristannyam pam yati ta propusknoyu zdatnistyu peredavannya Na suchasnih rozpodilenih obchislyuvalnih seredovishah takih yak MapReduce bulo rozrobleno specializovani algoritmi peremnozhuvannya Algoritmi dlya sitok Matrichne peremnozhuvannya vikonane za 2n 1 krokiv dlya dvoh matric n n na perehreshuvanij sitci Isnuye nizka algoritmiv dlya peremnozhuvannya na sitkah Dlya peremnozhuvannya dvoh n n na standartnij dvovimirnij sitci iz zastosuvannyam dvovimirnogo en obchislyuyut peremnozhuvannya v 3n 2 krokiv hocha ce znizhuyetsya do polovini cogo chisla dlya povtoryuvanih obchislen Standartnij masiv ye neefektivnim oskilki dani z dvoh matric ne nadhodyat odnochasno i jogo maye buti dopovnyuvano nulyami Rezultat ye she shvidshim na dvosharovij perehreshuvanij sitci de potribno lishe 2n 1 krokiv Produktivnist dodatkovo pokrashuyetsya dlya povtoryuvanih obchislen veduchi do stovidsotkovoyi efektivnosti Perehreshuvanij sitkovij masiv mozhna rozglyadati yak osoblivij vipadok ne planarnoyi tobto bagatosharovoyi obroblyuvalnoyi strukturi Div takozh en en en en Primitki 2008 Sorting and Searching The Algorithm Design Manual Springer s 45 46 401 403 doi 10 1007 978 1 84800 070 4 4 ISBN 978 1 84800 069 8 angl Tomas Kormen Charlz Lejzerson Ronald Rivest Kliford Stajn 2009 1990 28 Diyi nad matricyami Vstup do algoritmiv vid 3rd MIT Press i McGraw Hill ISBN 0 262 03384 4 Amarasinghe Saman Leiserson Charles 2010 6 172 Performance Engineering of Software Systems Lecture 8 MIT OpenCourseWare Massachusetts Institute of Technology originalu za 7 zhovtnya 2019 Procitovano 27 sichnya 2015 angl Lam Monica S Rothberg Edward E Wolf Michael E 1991 The Cache Performance and Optimizations of Blocked Algorithms Int l Conf on Architectural Support for Programming Languages and Operating Systems ASPLOS angl 1999 PDF Master s MIT Arhiv originalu PDF za 24 listopada 2019 Procitovano 15 grudnya 2019 angl Miller Webb 1975 Computational complexity and numerical stability SIAM News 4 97 107 CiteSeerX 10 1 1 148 9947 angl Press William H Flannery Brian P Vetterling William T 2007 Numerical Recipes The Art of Scientific Computing vid 3rd Cambridge University Press ISBN 978 0 521 88068 8 angl Le Gall Francois 2014 Powers of tensors and fast matrix multiplication Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation ISSAC 2014 arXiv 1401 7714 Bibcode 2014arXiv1401 7714L angl Pervinnij algoritm predstavlenij Donom Koppersmitom ta en 1990 roku maye asimptotichnu skladnist O n2 376 Jogo bulo polipsheno 2013 roku do O n2 3729 en sho dalo trivalist lishe troshki girshu nizh u vdoskonalennya le Gallya PDF Arhiv originalu PDF za 8 zhovtnya 2013 Procitovano 15 grudnya 2019 angl Iliopoulos Costas S 1989 PDF SIAM Journal on Computing 18 4 658 669 CiteSeerX 10 1 1 531 9309 doi 10 1137 0218045 MR 1004789 arhiv originalu PDF za 5 bereznya 2014 procitovano 15 grudnya 2019 The Coppersmith Winograd algorithm is not practical due to the very large hidden constant in the upper bound on the number of multiplications required angl Robinson Sara 2005 PDF SIAM News 38 9 Arhiv originalu PDF za 31 bereznya 2019 Procitovano 15 grudnya 2019 angl en On the complexity of matrix product In Proceedings of the thirty fourth annual ACM symposium on Theory of computing ACM Press 2002 DOI 10 1145 509907 509932 angl Henry Cohn en en and Chris Umans Group theoretic Algorithms for Matrix Multiplication arXiv math GR 0511460 Proceedings of the 46th Annual Symposium on Foundations of Computer Science 23 25 October 2005 Pittsburgh PA IEEE Computer Society pp 379 388 angl Henry Cohn Chris Umans A Group theoretic Approach to Fast Matrix Multiplication arXiv math GR 0307321 Proceedings of the 44th Annual IEEE Symposium on Foundations of Computer Science 11 14 October 2003 Cambridge MA IEEE Computer Society pp 438 449 angl Alon Shpilka Umans On Sunflowers and Matrix Multiplication 9 grudnya 2016 u Wayback Machine angl Randall Keith H 1998 PDF Ph D Massachusetts Institute of Technology s 54 57 Arhiv originalu PDF za 6 listopada 2020 Procitovano 15 grudnya 2019 angl Lynn Elliot Cannon A cellular computer to implement the Kalman Filter Algorithm Technical report Ph D Thesis Montana State University 14 July 1969 angl Hong J W Kung H T 1981 PDF STOC 81 Proceedings of the Thirteenth Annual ACM Symposium on Theory of Computing 326 333 Arhiv originalu PDF za 15 grudnya 2019 Procitovano 15 grudnya 2019 angl Irony Dror Toledo Sivan Tiskin Alexander September 2004 Communication lower bounds for distributed memory matrix multiplication J Parallel Distrib Comput 64 9 1017 1026 CiteSeerX 10 1 1 20 7034 doi 10 1016 j jpdc 2004 03 021 angl Agarwal R C Balle S M Gustavson F G Joshi M Palkar P September 1995 A three dimensional approach to parallel matrix multiplication IBM J Res Dev 39 5 575 582 CiteSeerX 10 1 1 44 3404 doi 10 1147 rd 395 0575 angl Solomonik Edgar Demmel James 2011 Communication optimal parallel 2 5D matrix multiplication and LU factorization algorithms Proceedings of the 17th International Conference on Parallel Processing Part II 90 109 angl Bosagh Zadeh Reza Carlsson Gunnar PDF Arhiv originalu PDF za 14 lipnya 2014 Procitovano 12 lipnya 2014 angl Bae S E T W Shinn T Takaoka 2014 A faster parallel algorithm for matrix multiplication on a mesh array 18 travnya 2021 u Wayback Machine Procedia Computer Science 29 2230 2240 angl Kak S 1988 A two layered mesh array for matrix multiplication Parallel Computing 6 383 385 angl Kak S 2014 Efficiency of matrix multiplication on the cross wired mesh array https arxiv org abs 1411 3273 23 bereznya 2019 u Wayback Machine angl Kak S 1988 Multilayered array computing Information Sciences 45 347 365 angl LiteraturaButtari Alfredo Langou Julien Kurzak Jakub 2009 A class of parallel tiled linear algebra algorithms for multicore architectures Parallel Computing 35 38 53 arXiv 0709 1272 doi 10 1016 j parco 2008 10 002 angl Goto Kazushige van de Geijn Robert A 2008 Anatomy of high performance matrix multiplication ACM Transactions on Mathematical Software 34 3 1 25 CiteSeerX 10 1 1 140 3583 doi 10 1145 1356052 1356053 angl Van Zee Field G van de Geijn Robert A 2015 BLIS A Framework for Rapidly Instantiating BLAS Functionality ACM Transactions on Mathematical Software 41 3 1 33 doi 10 1145 2764454 angl How To Optimize GEMM 15 veresnya 2020 u Wayback Machine angl