Сообщество

интересная программа, спасибо! сделаю пару замечаний, если речь идет о перемножении матриц трансформации (наверное, это предположение я сделал из-за упоминания в тексте OpenGL).
Матрицы трансформации перемещают, меняют объект в трехмерном пространстве и имеют следующий вид:

[ R_xx, R_xy, R_xz, T_x ]
[ R_yx, R_yy, R_yz, T_y ]
[ R_zx, R_zy, R_zz, T_z ]
[   0,   0,    0,    1  ]

[ X ]
[ Y ]
[ Z ]
[ 1 ]

[ 0, 0, 0, 1] 

В этом случае, мне кажется, в программе нужно использовать свойства матрицы: Последняя строка имеет всегда один и тот же вид, и соответственно, ее не нужно умножать - ответ известен.
четвертый член в скалярном произведении строки матрицы a на столбец b тоже известен и равен нулю для первых 3х строк матрицы a и единицы для последней.
Т.е. вообще говоря для задания матриц трансформации, достаточно задания 3х4 значений.

Кроме того у матрицы трансформации есть антисимметрия относительно диагонали. Я не знаю как его использовать в 3х мерном случае, однако, в моей задаче, которая сводится к поворотам на двумерной плоскости, я сначала решил применить матрицу трансформации

[  cos(a) sin(a) 0 ]
[ -sin(a) cos(a) 0 ]
[   0       0    1 ]

[  cos(a) -sin(a) 0 ]   [ X ]   [ X cos(a) - Y sin(a) ]
[  sin(a)  cos(a) 0 ] x [ Y ] = [ Y cos(a) + X sin(a) ]
[   0        0    1 ]   [ 1 ]   [          1          ]

(A, B) x (C, D) = (AC - BD, AD + BD)

PS: Наконец сам себе сформулировал (для своей двумерной задачи) что двигаюсь правильно.

Да, речь идёт именно о матрицах трансформации. Понятно, что нижняя строка у таких матриц всегда (0, 0, 0, 1), однако в данном случае ничего от этого факта выиграть не получится, т.к. кол-во операций для скалярного произведения 3-х и 4-х компонентных векторов одинаково.

    madd @va04, @vb01
    madd @va8C, @vb23
    vc_0 := addf @1, @2
    ; против
    madd @va04, @vb01
    addf @va_8, @vb_2
    vc_0 := addf @1, @2

А раз так, то пусть уж честно считает. =) Собственно, мои вопросы касаются архитектурно-специфических оптимизаций, а не алгоритмических.

Вектор трехмерный, хотя и описывается четырьмя координатами, четвертая координата у вектора всегда равна 1.

У вектора (как направления) четвёртый компонент всегда 0, это у точки 1, что, при умножении на матрицу, позволяет ей принимать перемещение, а вектору это не нужно.

Повороты в 3-х мерном пространстве гораздо сложнее и могут быть заданы разными методами (углы Эйлера для всех случаев XYZ, XZY, YXZ и т.д, кватернионы).

А это:

    x = X cos(a) - Y sin(a)
    y = Y cos(a) + X sin(a)

классический способ поворота двумерной точки (http://ru.wikipedia.org/wiki/Поворот). Если у Вас нет сложной иерархии объектов, то его лучше и использовать. Матрицы полезны, чтобы собирать сложные иерархические трансформации в одну сущность.

Однако, произведение последней строки a [0,0,0,1] на столбцы b можно не считать.
это не нужно:

    madd @va37, @vb01
    madd @vaBF, @vb23

    madd @va37, @vb45
    madd @vaBF, @vb67

    madd @va37, @vb89
    madd @vaBF, @vbAB

    madd @va37, @vbCD
    madd @vaBF, @vbEF

и даже не инициализировать

va37 := pack @va76, @va23
vaBF := pack @vaFE, @vaAB

хотя, я понимаю, что основной вопрос в оптимизации. но, смотрите, сейчас в программе около 90 строк, а нужно для параграфа 64. Десять строк я уже убрал.

Еще более радикально, представьте что матрица трансформации инициализирована в порядке

[ 0  3  6  9 ]
[ 1  4  7  A ]
[ 2  5  8  B ]

qewerty, в вашем коде для достижения большего быстродействия можно поставить все команды madd подряд, а затем после них сделать сложение.
Запись программы в один параграф, на мой взгляд, в данном случае не должна принести увеличения скорости, но это лучше проверить с помощью таймера.
С записью в память по 64 или 32 бита тоже можно поэкспериментировать и посмотреть результат по таймеру.
У вас отладочная плата имеется своя или используете удалённый доступ?

krufter_multiclet,

в вашем коде для достижения большего быстродействия можно поставить все команды madd подряд, а затем после них сделать сложение

Да, у меня были такие мысли, но с учётом того, что команды исполняются по готовности потребителя, я сомневался. Спасибо!

У вас отладочная плата имеется своя или используете удалённый доступ?

Я пока только на модели запускал. Плату заказал недавно, сегодня из курьерской службы звонили, сказали завтра привезут. На днях смогу потестировать.

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

trott, ну если допустить, что каждая команда выполняется за один такт, то эта функция на 4-х клеточном процессоре в идеале должна занять около 25 тактов, что при 80 МГц должно дать 3.2 мильёна умножений в секунду. =)

qewerty, это вопрос к производителю, сколько тактов у него занимает каждая операция. кстати, очень хотелось бы узнать, cколько?
я вчера оценивал результат. у меня получилось, в предположении что умножение/сложение занимают 10 тактов. 100х100х10^3х10=10^7 или 10 миллионов. поделим на 4 клетки = 2.5 мильёна

Примеры работы с таймером имеются конечно. Могу прикрепить сюда пример.

да, пожалуйста, если можно пример ... и есть ли данные, сколько тактов занимает madd, pack, rdl, и прочие.

С помощью таймера периферии:

.include "HDL50001_pcf.inc" 
.alias ron1 1
.alias ron3 3

.text

initTimer0:
    jmp StartTimer0
    getl 0xFFFFFFFF
    wrl @1, TIM0_CNT ;период
    getl 0x00000001
    wrl @1, TIM0_PSC ;значение предделителя
complete

StartTimer0:
    jmp Get_start_time
    getl 0x00000005
    wrl @1, TIM0_CR
complete

Get_start_time:
    rdl TIM0_CNTVAL
    setl #ron1, @1
    jmp code
complete

code:
    ;code here
    jmp Get_finish_time
complete

Get_finish_time:
    rdl TIM0_CNTVAL
    getl #ron1
    subl @1, @2
    setl #ron3, @1
    jmp initUART0
complete

В текущей версии у таймера периферии необходимо период умножить на два, чтобы получить количество тактов(у системного таймера делать это не нужно).

В потоке на каждом такте выдаётся результат выполнения команды, например madd, rdl, pack... Точное время выполнения команд по отдельности может быть вычислено только с помощью пользовательской потактовой модели(или мы можем проанализировать параграф на своей модели и сказать сколько и что выполняется). Можем только составить чистое время выполнения команд в ALU без учёта выходного диспетчера.

На LDM-MCp0411100101-Q208 Evolution получилось выжать около мильёна умножений в секунду на 80 МГц. Функция вот такая: http://dumpz.org/605010/. Выполняется около 80 тактов, точнее сказать сложно из-за расходов на организацию цикла, кстати, тоже на асме, т.к. на Си слишком большие накладные расходы были, практически равные умножению.

А вот можно посмотреть аналогичный тест на Raspberry Pi: http://thinkingeek.com/2013/05/12/arm-assembler-raspberry-pi-chapter-14/. Лучший результат (асм+неон) смог умножить 2 в 21 степени матриц (2097152) за 1.51 сек. И это на 700 МГц (до 1 ГГц в турбо режиме, но будем брать 700). Для 80 МГц получается 158725 умножений в сек, что в 6 раз меньше чем у Мультиклета на МГц.

Хочу ещё попробовать сделать более жизненные тесты и сравнить с Raspberry Pi, которая у меня тоже имеется.

Спасибо за приятную информацию и проведённые тесты.
А не покажете как цикл организовали? Или весь исходник? Т.е. вызов функции был из Си?

krufter_multiclet, цикл, как выше писал, на асме тоже. Как-то так (воспроизвожу по памяти, т.к. сейчас не могу посмотреть исходники):

; Си-прототип (возвращает замеренное таймером значение):
;     int qf_test_mat44_mul(const float* a, const float* b, float* c);

.global qf_test_mat44_mul
qf_test_mat44_mul
    jmp start_timer
    getl 0xFFFFFFFF
    wrl @1, TIM0_CNT
    getl 0x00000001
    wrl @1, TIM0_PSC
complete

.local start_timer
start_timer
    jmp get_start_time
    getl 0x00000005
    wrl @1, TIM0_CR
complete

.local get_start_time
get_start_time
    rdl TIM0_CNTVAL
    setl #1, @1
    jmp qf_loop_mat44_mul
complete

.local check_timer_and_quit
check_timer_and_quit:
    rdl TIM0_CNTVAL
    getl #1
    subl @1, @2
    rdl #SP
    jmp @1
    wrl @3, #SP,4
complete

.global qf_loop_mat44_mul
qf_test_mat44_mul:
    jmp  qf_loop_mat44_mul.CHECK
    getl 1000000
    setl #0, @1
complete

.local qf_loop_mat44_mul.CHECK
qf_test_mat44_mul:
    c := getl #0
    je   @c, check_timer_and_quit ; как ваше Get_finish_time, вычисляет разницу таймера, но выходит по #SP
    jne  @c, qf_mat44_mul
    subl @c, 1
    setl #0, @1
complete

И в самой функции qf_mat44_mul, в параграфе P0, возврат делается на в qf_test_mat44_mul.CHECK, а не по #SP:

.global qf_mat44_mul
qf_mat44_mul:
    .alias qf_mat44_mul.a #SP,8
    .alias qf_mat44_mul.b #SP,12
    .alias qf_mat44_mul.c #SP,16

    jmp qf_mat44_mul.P0
    ; Сохранение адресов аргументов в регистрах
    rdl qf_mat44_mul.a
    rdl qf_mat44_mul.b
    rdl qf_mat44_mul.c
    setl #32, @3
    setl #33, @3
    setl #34, @3
complete

.local qf_mat44_mul.P0
qf_mat44_mul.P0:
    ; вот так:
    jmp qf_loop_mat44_mul.CHECK ; лучше было бы в регистре сохранить адрес возврата
    ; вместо:
    ; rdl #SP
    ; jmp @1
    ...

Всё это не по-сишному конечно, но в частных случаях вполне возможно, к тому же оптимизирующий компилятор также может не оформлять стековый кадр если функции он не нужен.

А qf_test_mat44_mul вызывалась уже из Си.

Чего-то я напутал в первый раз. Пост выше поправил. Вечером до дома доберусь, выложу исходники.

qewerty, я протестировал свой микроконтроллер stm32f4-discovery(87мГц) на умножение матриц 4х4. результат: 5 секунд для 1 мильона умножений. код на C выложен на https://github.com/trot-t/stm32-matrix . плавающая точка включена.

Интригующие тесты. Во внутричиповой памяти Мультиклета памяти больше чем у Спектрума. Если можно будет подключить монитор, то возможно написать демку с 3D сценой, на подобии из тех, что были на спектруме.

Да на VGA монитор выводили красный цвет, дальше времени не было.

Прошлый тест совсем уж наивный был, выкладывать не стал. Кое-что доработал, исходники можно взять тут: https://bitbucket.org/qewerty/qf/overview. Функция умножения теперь сохраняет и восстанавливает значения регистров и судя по всему должна работать корректно. Да и цикл теперь более честный, с пробросом аргументов через стек.

Кому интересно может собрать, проверить и убедиться:

hg clone https://bitbucket.org/qewerty/qf
cd qf
make
mc-ploader qf_mat44_mul_TEST

Будет мигать диодами с интервалом в миллион матричных умножений. =)

Перемножение матриц - задача с явным параллелизмом, так что это не супер аргумент в пользу мультклета. Потому как гипотетический микроконтроллер с 4 обычными ядрами так же покажет пропорциональный рост производительности на этой задаче. И тут уже надо мерить не число умножений/мегагерц, а число умножений/площадь кристалла или число умножений/ватт. и так далее. Т.е. приводить задачу к объему вентелей, требуемых для ее решения.

Так же очень интересно, а что будет с перемножением матриц хотябы 100 х 100 например, что бы понять как работает подсистема памяти.

Sergei.S,

Перемножение матриц - задача с явным параллелизмом

Да неужели? В том то и суть. Важно КАК оно распараллеливается. А задачи где параллелизма нет, нельзя распараллелить никак.

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

Да, только на потоки код должен разбивать программист, а это гораздо сложней чем написать рядом несколько независимых инструкций. А если матрицы нужно умножать строго последовательно, в виду их зависимости, то никакого прироста не будет (в графе сцены, например, они частично зависимы). А если одно умножение на несколько потоков раздробить, то может быть даже ещё медленней из-за расходов на синхронизацию.

Хочу добавить свои пять копеек. Пишу реализацию delay-функции на асме. Пробы пера. Первый вариант имеет оверхед с лишними вызовами двух команд в пределах одного параграфа:

.include "HDL50001_pcf.inc" 

.alias SP 39

.global MCp_delay

.text

MCp_delay:
    rdl #SP, 8
    setl #32, @1
    jmp loop
complete

loop:
    getl #32
    je @1, return
    jne @2, loop
    subl @3, 1
    setl #32, @1
complete

return:
    rdl #SP
    jmp @1
complete

.include "HDL50001_pcf.inc" 

.alias SP 39

.global MCp_delay

.text

MCp_delay:
    rdl #SP, 8
    setl #32, @1
    jmp loop
complete

loop:
    getl #32
    je @1, return
    jne @2, decrement
complete

decrement:
    getl #32
    subl @1, 1
    setl #32, @1
    jmp loop
complete

return:
    rdl #SP
    jmp @1
complete

.include "HDL50001_pcf.inc" 

.alias SP 39

.global MCp_delay

.text

MCp_delay:
    rdl #SP, 8
    setl #32, @1
    setq #33, 0
    setl #MODR, 2  ; пробовал маску вплоть до 15
    jmp loop
complete

loop:
    jmp loop
    getl #32
    xor @1, #33    ; в модели регистр 0x21 принимал значение 0xf_0000_000f и не менялся от вызова к вызову параграфа
                   ; какой есть альтернативный способ сравнения двух значений?
    je @1, return
complete

return:
    rdl #SP
    jmp @1
complete

    jmp loop
    […]
    je @1, return

mouse, в вашем примере с автоинкрементом установка значения в #MODR разрешает инкрементировать индексный регистр #33, а маска, которая записывается в сам регистр #33, запрещает автоинремент.

MCp_delay:
    rdl #SP, 8
    setl #32, @1
    setq #33, 0  ; <= здесь выставили базу = 0, маску = 0, индекс = 0. если маска во всех разрядах равна 0, то инкремента не будет, 
; из-за того, что он запрещен маской, 0000000000000000 - маска запретила изменять все разряды у индекса.  
; 11111111111111 - такая маска разрешит менять все разряды индекса в регистре   
    setl #MODR, 2  ; пробовал маску вплоть до 15 <== здесь говорим, что будем автоинкрементировать единственный индексный регистр #32, это не маска, а резрешение инкремента
    jmp loop
complete

mouse, такая программа в модели 3 раза делает автоинкремент регистра #32

.include "HDL50001_pcf.inc" 
.text

.alias SP 39

.global MCp_delay

.text

MCp_delay:
;    rdl #SP, 8
    getl 2 ; сделаем 3 шага автоинкремента для проверки в модели
    mask := getl 65535        ; готовим маску 1111111111111111
    ind_aux := subl @mask, @2   
    ind := slll @ind_aux, 16  ; индекс  
    addl @ind, @mask ; заполняем индекс и маску
    patch @1, 0      ; склеиваем маску индекс и базу в 64 разряда
    setq #32, @1     ; записываем значение в индексный регистр #32
    setl #MODR, 1   ; разрешаем модификацию первого индексного регистра, он имеет #32
    jmp loop
complete

loop:
    exa #32
    je @1, return ; проверка на 0
    jne @2, loop
complete

return:
    getl 0x1000 ; остановимся, компилятор не любит пустых параграфов
;    rdl #SP
;    jmp @1
complete

trott wrote:

mouse, в вашем примере с автоинкрементом установка значения в #MODR разрешает инкрементировать индексный регистр #33, а маска, которая записывается в сам регистр #33, запрещает автоинремент.
[...]
setl #MODR, 2 ; пробовал маску вплоть до 15 <== здесь говорим, что будем автоинкрементировать единственный индексный регистр #32, это не маска, а резрешение инкремента

Если обратиться к документации , то выставленный бит разрешает инкремент. Что значит "маска / не маска"? Здесь я вас не понимаю. Про регистр #32, думаю, опечатка, т.к. 2-й бит отвечает за регистр #33.

Давайте вернёмся к моему начальному сообщению:

mouse, […]Вариант с авто-инкрементом индексных регистров не подходит[…]

Мы имеем 16-ти битный автоинкрементирующийся счётик . Как мне сделать точную паузу более 5 мс при частоте процессора 100MHz? Есть ещё один тонкий момент с автоинкрементом. В момент вызова MCp_delay из Си мне придётся сохранять значение как оригинального MODR, так и индексного регистра, дабы после завершения функции вернуть всё на место (я же не знаю, кем этот регистр может быть ещё использован). Т.е. функция откушает 8 (а то и все 12) байт на ровном месте + ненужные операции над памятью.