Блочное LU-разложение

Добрый день. Меня интересует, как видно из темы, как правильно сделать блочное LU-разложение.

История такая: в известной книжке "Матричные вычисления" была найдена блочная версия сабжа:
1) Матрица делится на 4 квадрата
2) Для верхнего левого (A11) выполняется LU-разложение обычным алгоритмом
3) Для правого верхнего (A12) и левого нижнего (A21) решаются треугольные системы
4) Из правого нижнего (A22) вычитается произведение A21*A12
5) Повторяем алгоритм для A22

Я взял библиотеку ViennaCL, и чуть-чуть поменял их kernel-ы - так, чтобы вычисления проводились на нужных мне индексах матрицы. Запускаю - а работает, во-первых, медленнее, чем даже процессорная версия, во-вторых, с увеличением размера матрицы скорость очень сильно падает. Получается так, что наилучшая производительность выходит при размере блока = размеру матрицы, т.е. когда выполняется обычное, не блочное разложение.
Результат работы алгоритма, кстати, верный.

Хочу спросить, имеет ли смысл делать именно блочное разложение, и, если да (а думается, что все-таки да), то какой должен быть правильный алгоритм? Как считать, чтобы был выигрыш в производительности?
Спасибо

Forums: 

Системы каких размеров вы

Системы каких размеров вы тестировали?

2) Для верхнего левого (A11) выполняется LU-разложение обычным алгоритмом
Каким именно?

4) Из правого нижнего (A22) вычитается произведение A21*A12
Какими средствами вычисляется произведение? Сделайте измерение времени вычисления этого произведения, и скажите результаты и размеры перемножаемых матриц

1) Матрица делится на 4 квадрата
Может всё-таки A21 и A12 могут быть прямоугольниками?

5) Повторяем алгоритм для A22
Как именно повторяем? Рекурсивно с уменьшением размера блока или уже не блочный алгоритм?
Или опять же A21 и A12 могут быть прямоугольниками?

Хочу спросить, имеет ли смысл делать именно блочное разложение, и, если да (а думается, что все-таки да), то какой должен быть правильный алгоритм? Как считать, чтобы был выигрыш в производительности?

В книге говорится, что при большом отношении n/r (r - размер блока, n - размер системы), почти все арифметические операции выполняются в контексте умножения матриц.
Как известно умножение матриц(gemm) на GPU может достаточно эффективно выполнятся, то есть по-идеи есть смысл в попытках использовать блочный алгоритм(обязательно нужно проверить каким образом ViennaCL делает LU разложение - может там уже блочный вариант?). Также известно что для разных размеров матриц gemm может выполнятся с разной эффективностью(которая обычно оценивается flop/s).
В итоге вам нужно следить за двумя вещами:
1) чтобы на вход gemm поступали матрицы с такими размерностями, для которых ваше железо показывает приемлемые результаты gemm flop/s.
2) отношение n/r должно быть достаточно большим, чтобы bottleneck'ом было gemm - в книге приводится выражение для вычисления доли gemm в общих flop'ах, в зависимости от n/r