nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

UDUT-разложение (оно же LTDL)

Это как LDLT-разложение симметричной положительно определённой матрицы A, но только UDUT, т.е верхняя унитреугольная матрица (треугольная с единицами на главной диагонали) умножается на диагональную и на себя транспонированную. Его можно было бы также назвать LTDL-разложением, тут устоявшейся терминологии нет.

Такое разложение ничем принципиально не отличается от LDLT, но для этого разложения и всех последующих операций при инвертировании матрицы A, мы будем идти по элементам "в обратном порядке", с правого нижнего угла к левому верхнему. Для системы команд QuatCore, где циклы "вниз" организуются проще всего с помощью "команд" iLOOP, jLOOP и kLOOP (если индексная переменная равна нулю - идти дальше, иначе уменьшить её на единицу и прыгнуть в начало цикла), это предпочтительнее.


Итак, у нас опять есть матрица A - симметричная, положительно определённая:


В этот раз мы ведём нумерацию индексов с нуля, чтобы стать "ближе к железу". Верхний треугольник мы не записываем, чтобы напомнить: в памяти хранится только нижний треугольник. Это не обязательно означает треугольную адресацию:


когда нижний треугольник занимал бы последовательные ячейки памяти


тоже была идея реализовать такую адресацию "аппаратно", просто в виде таблицы:
itri[i]
00
11
23
36
410
515


(у меня матриц крупнее 6×6 здесь не ожидается)
Это весьма экономично.
[Код на верилоге]
module TriangleNum (input [2:0] D, output [3:0] Q);

assign Q = 	(D == 0)? 4'd0:
		(D == 1)? 4'd1:
		(D == 2)? 4'd3:
		(D == 3)? 4'd6:
		(D == 4)? 4'd10:
		(D == 5)? 4'd15:
			  4'bxxxx;
endmodule


Как ни странно, этот код синтезируется в 3 ЛЭ, так как выходит Q[3] = D[2].


Есть в такой адресации определённая красота: в формуле, преобразующей индексы в адрес в памяти, нигде не фигурируют размеры матрицы! Но можно и самой обычной адресацией обойтись, а в верхний треугольник запихать прочие величины, как в таблице памяти данных ВИПС. Буду ещё громко думать, что лучше.

Далее мы хотим эту матрицу A представить в виде произведения:


Две эти матрицы умещаются в памяти, выделенной под матрицу A:

и это разложение может быть выполнено "на месте", не требуя временного хранилища.

Чтобы понять, как это сделать, помножим 3 матрицы и посмотрим, как их коэффициенты соотносятся с коэффициентами исходной матрицы:




Разумеется, и здесь мы получили симметричную матрицу, так что уберём ненужный верхний треугольник:


Как и ожидалось, в этот раз самые простые выражения расположены в правом нижнем углу и усложняются по мере движения к левому верхнему (в прошлый раз всё было наоборот).

Можно идти в разном порядке. В прошлый раз мы двигались по столбцам, но как оказалось, для эффективных вычислений нужно было вводить временные переменные. Будем двигаться по столбцам и в этот раз - упорство должно быть вознаграждено!

Начинаем с правого столбца.


При работе "на месте" делать это присвоение не требуется - d3 расположен в точности там же, где a33.

Берёмся за следующий (предпоследний) столбец, причём двигаемся снизу вверх:


откуда:




Берёмся за следующий столбец:



откуда:





Наконец, обрабатываем левый столбец:





откуда:





Вот и всё. Глядя на эти формулы, можно обнаружить, что введя вспомогательный массив из N-1 значений (N×N-размер матрицы, в рассмотренном случае N=4), мы обойдёмся без львиной доли умножений.

Покажем, как это делается. Запишем выражения для предпоследнего столбца:




Для столбца с индексом 1:






И наконец, для левого столбца (с индексом 0):








Выпишем все эти выкладки в таблицу. Вычисления здесь ведутся сверху вниз и только потом слева направо:

 
 
  
  


В такой записи получается, что значения dn нам нужны только в знаменателе. Если мы делаем это разложение, чтобы обратить матрицу "на месте", то мы можем прямо "на лету" заносить вместо самих значений dn обратные, как-то так:

 
 
  
  
  


В такой записи нам нужно всего лишь N взятий обратной величины, так что даже при отсутствии аппаратного деления нас это не шибко затормозит (метод Ньютона - наше всё!). Размер вспомогательного массива возрос до N. По идее, можно обращать dj "сразу же", не запихивая их сначала в tj, но так получается более логичный код.

Все остальные выражения очень хорошо ложатся на наше "умножение с накоплением" (Fused Multiply-Add, FMA, и Fused Multiply-Subtract, FMS), тогда как исходные выражения этому не удовлетворяли, поскольку требовали перемножить 3 значения, и лишь затем результат прибавлять к результату. В результате мы не только получаем компактный и быстрый код, он ещё и обладает хорошей точностью за счёт того, что накопление производится во всех 32 битах, и лишь результат записывается в 16 бит.

Не будем подсчитывать отдельно сложения и умножения, нас интересует число операций FMA/FMS.

При обработке столбца под номером n, нам требуется n-1 умножений, чтобы временные значения tj помножить на диагональные элементы. Ещё мы вычисляем n временных значений, причём на первое не требуется арифметических операций вообще, на второе - одна операция FMS, на третье - 2 операции, и на значение n: n-1 операций. В общей сложности это получается n(n-1)/2, и ещё n-1, это получается (n-1)(n+2)/2 = n2/2+n/2-1.

Просуммировав n от 1 до N, получаем количество операций FMA/FMS: N3/6+N2/2-(2/3)N. Асимптотически совпадает с тем, что мы посчитали в LDLT-разложении (т.е коэф 1/6 при N3), но на малых значениях N чувствуется существенный выигрыш:

При N=4 в прошлый раз у нас вышло 6 делений, 16 умножений и 10 сложений. В этот раз: 4 взятия обратной величины и 16 операций FMA/FMS.
При N=6 в прошлый раз получалось 15 делений, 50 умножений и 35 сложений, а теперь: 6 взятий обратной величины и 50 операций FMA/FMS.
Причём, это ЕДИНСТВЕННЫЕ 6 взятий обратной величины, которые нужны для обращения матрицы A, больше не будет!

Понятно, что разница здесь не в том, используется ли LDLT или UDUT, а в последовательности выполнения операций и способа введения вспомогательного массива.
Tags: ПЛИС, математика, работа, странные девайсы
Subscribe

  • Нахождение двух самых отдалённых точек

    Пока компьютер долго и упорно мучал симуляцию, я пытался написать на ассемблере алгоритм захвата на ближней дистанции. А сейчас на этом коде можно…

  • Слишком общительный счётчик

    Вчера я чуть поторопился отсинтезировать проект,параметры не поменял: RomWidth = 8 вместо 7, RamWidth = 9 вместо 8, и ещё EnableByteAccess=1, чтобы…

  • Балансируем конвейер QuatCore

    В пятницу у нас всё замечательно сработало на симуляции, первые 16 миллисекунд полёт нормальный. А вот прошить весь проект на ПЛИС и попробовать "в…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 11 comments

  • Нахождение двух самых отдалённых точек

    Пока компьютер долго и упорно мучал симуляцию, я пытался написать на ассемблере алгоритм захвата на ближней дистанции. А сейчас на этом коде можно…

  • Слишком общительный счётчик

    Вчера я чуть поторопился отсинтезировать проект,параметры не поменял: RomWidth = 8 вместо 7, RamWidth = 9 вместо 8, и ещё EnableByteAccess=1, чтобы…

  • Балансируем конвейер QuatCore

    В пятницу у нас всё замечательно сработало на симуляции, первые 16 миллисекунд полёт нормальный. А вот прошить весь проект на ПЛИС и попробовать "в…