Прошу помощи по обсчёту математики

AXT> Если выкинуть всякую малозначащую труху, то мне надо умножать вектор на матрицу. Много раз.

Миллион раз в секунду?

Тупо: а мультимедийные расширения процессора нельзя дёрнуть? SSE?

AXT> Кто что посоветует?
Попробуй использовать LAPACK - там вроде должны быть алгоритмы оптимизированные.
Библиотека есть для Си и Фортрана точно.

Серокой> Тупо: а мультимедийные расширения процессора нельзя дёрнуть? SSE?

Да я как раз и думаю в эту сторону, но есть non-intel платформы, плюс (самое страшное) данные лежат не подряд.

Massaraksh> Миллион раз в секунду?

Около 20 млн в идеале. Матрица 4*4 — немного, но надо же синхроизироваться, и вот тут залежи граблей ...

(правка) 4*4 она после беспощандной редукции... похоже дальше никак.

Massaraksh>> Миллион раз в секунду?
AXT> Около 20 млн в идеале. Матрица 4*4 — немного, но надо же синхроизироваться, и вот тут залежи граблей ...
AXT> (правка) 4*4 она после беспощандной редукции... похоже дальше никак.

Крайгинг-интерполяция?
А если на матрицу сделать маску, и не проводить лишних умножений?

AXT> Но матрица имеет совсем немного значений элементов и почти регулярную структуру, так что есть подозрения, что там скрыт более простой алгоритм. Ранг точно меньше 32, и сильно.
AXT> Кто что посоветует?
Первое что сделать - выделить линейно-независимые элементы, поискать в матрице какие-то симметрии, пусть неочевидные (симметрия с перестановками столбцов/строк - тоже симметрия). Ну и т.п. - поиграться с математикой, короче.
Можно даже втупую: выписать (машиной, конечно ) полные выражения для элементов результирующего тензора, подставив матрицу как константы и упростить их (поэлементно) машиной.

На математике можно выиграть порядки.
А на библиотеках/SIMD - в лучшем случае только разы.

Если математика не помогает, то я б смотрел на CUDA/OpenCL (в зависимости от железа). Там такие задачи с массивным перемножением всего на всё бегут замечательно.

AXT,

Пардон, не видел. Татарин уже ответил по существу. Я только его ответ прокомментирую.

Татарин> Первое что сделать - выделить линейно-независимые элементы, поискать в матрице какие-то симметрии, пусть неочевидные (симметрия с перестановками столбцов/строк - тоже симметрия). Ну и т.п. - поиграться с математикой, короче.

Это в первую очередь. Согласен. Идеально бы было найти ранг, но это может быть довольно сложная задача.

Татарин> Можно даже втупую: выписать (машиной, конечно ) полные выражения для элементов результирующего тензора, подставив матрицу как константы и упростить их (поэлементно) машиной.

Даже матрица 32 на 32 и вектор на 32 — уже можно сделать 32 линейные записи по 32 элемента в каждом — распрямить циклы. Если в матрице много 0, 1, 2, то эти выражения вместо умножения легко оптимизируются ручками (для 2 в целочисленной арифметике). Если не делал преобразований матрицы, а есть одинаковые строчки (для вектора столбца), то имеет смысл сравнить значения соответствующие в векторе и не перевычислять.

Татарин> На математике можно выиграть порядки.

Я бы посоветовал погонять умножения и разные запросы в Maple. Она может предоставить совершенно нетривиальный вариант упрощения умножения для конкретной матрицы (матрица, как я понял, константная).




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

AXT, а какие данные (целочисленные? с плавающей? с фиксированной? диапазон значений?), какая точность (если с фиксированной, то может легче смоделировать целочисленной?)?

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

Mishka> Даже матрица 32 на 32 и вектор на 32 — уже можно сделать 32 линейные записи по 32 элемента в каждом — распрямить циклы.
Там возможно проще N раз клонировать входной вектор в вектор с дубликатами. А константную матрицу также в виде вектора представить - NxN длиной. Потом просто один пробег попарного умножения плюс прибавления к результату с попутным вычислением целевой ячейки - чтобы она на N+1 элементе меняла значение.
Ну или просто попарное умножение, а потом N вызовов ф-ии сложения элементов вектора с указанием диапазона.

Mishka>> Даже матрица 32 на 32 и вектор на 32 — уже можно сделать 32 линейные записи по 32 элемента в каждом — распрямить циклы.
yacc> Там возможно проще N раз клонировать входной вектор
Для общего случая есть алгоритм Штрассена.
Но АХТ утверждает, что у него матрица какая-то неполноценная - с нулевыми элементами и т.п, и вообще константа. И тут на самом деле может быть, где гульнуть.

101> Автор, кинь матрицу и вектора. Посмотрят наши, что нибудь посоветуют.

Сейчас уже не имеет смысла — ещё летом всё было заборото оптимизацией "в лоб", так что интересно только в смысле поиграться с математикой. Матрицы в стандарте ITU-T H.264, если нет доступа к, то могу выложить матрицы. Они там реально странные, вроде как закономерности есть даже на глаз, но единого принципа я так и не нашёл.

А сейчас попал на похожую задачу в рамках хобби Там тоже 4*4, но источник другой.

И всем, разумеется, спасибо за советы

Татарин> И тут на самом деле может быть, где гульнуть.
Навскидку linear algebra - Recommendations for a usable, fast C++ matrix library? - Computational Science Stack Exchange

Mishka> Massaraksh, я думаю, что ручками оптимизоровать выражения при константной матрице будет куда более эффективно. При этом даже можно врубить для данноё ф-ции максимальную оптимизацию выражений при компилировании данной ф-ции по быстродействию в ущерб размеру. Она всё равно уложится в кэш у современных процессоров. Это будет лучше маски, которую надо прикладывать динамически.
Просто из контекста не совсем понятно, насколько она константная. У меня недавно была задача, как раз связанная с крайгинг-интерполяцией, где тоже приходилось много раз умножать вектор на константную матрицу, которая была константной аж 1/24 секунды.

yacc> Там возможно проще N раз клонировать входной вектор в вектор с дубликатами.

Зачем? Какой смысл? Если предпологать многопоточные процы, то векторы входной и на момент выполнения операцыии является тоже константным. По чтению нет никакой конкунренции. Или ты предлагаешь анализ вектора на дубликаты? Если последнее, то при разных строчках матрицы (для вертикального вектора, для горизонтального — стобцы матрицы), выход будет другим. 32 не настолко большая цифра, чтобы можно было бы выиграть на динамическом анализе, если вообще можно выиграть.

yacc> А константную матрицу также в виде вектора представить - NxN длиной. Потом просто один пробег попарного умножения плюс прибавления к результату с попутным вычислением целевой ячейки - чтобы она на N+1 элементе меняла значение.

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

yacc> Ну или просто попарное умножение, а потом N вызовов ф-ии сложения элементов вектора с указанием диапазона.
Блин, зачем попарные умножения? Вызов стоит дорого. Ты представляешь, что такое динамический поиск общих подвыражений, которые ты пытаешься сделать?

Massaraksh> Просто из контекста не совсем понятно, насколько она константная. У меня недавно была задача, как раз связанная с крайгинг-интерполяцией, где тоже приходилось много раз умножать вектор на константную матрицу, которая была константной аж 1/24 секунды.

Я думаю, что даже в этом случае распрямлённая матрица даст выигрышь по сравнению с поппыткой проанализировать матрицу с целью нахождения масок. Да и проверка маски денег стоит. Но тут надо эксперементировать.

yacc>> Там возможно проще N раз клонировать входной вектор в вектор с дубликатами.
Mishka> Зачем? Какой смысл?
Затем, что матрица последовательно лежит в памяти - дать ей такой же вектор для попарного умножения типа a[i] * b[i}

Mishka> Распрямение умножение — это 32 выражения вычисления
Когда я возился картами рельефа основная просадка была ... на ветвлении. Умножение выполняется быстро.

Mishka> Блин, зачем попарные умножения? Вызов стоит дорого. Ты представляешь, что такое динамический поиск общих подвыражений, которые ты пытаешься сделать?
Тогда самый быстрый способ - динамическая подлинковка.
Т.е. файл, который содержит ф-ию умножения на константную матрицу - отдельный.
Например mult.c содержит ф-ию f( value_type_ptr _inPtr, value_type_ptr _outPtr ) ( чисто к примеру )
А сам файл формируется динамически - т.е. некоторый препроцессор его формирует по константной матрице. сама ф-ия просто содержит набор вида _out[10] = _in[5]*10 + _in[6]*11 + ...
никаких ветвлений.
самое быстрое вычисление.

yacc> Затем, что матрица последовательно лежит в памяти - дать ей такой же вектор для попарного умножения типа a[i] * b[i}

Ы? Что это даст? Какая выгода?

И последовательно лежит — это как? В том же FORTRAN-e матрица лежит по столбцам, а в PL/I, C, C++ — по строкам. В Алголе 68 она вообще может быть разбросана по памяти поэлементно. Ну и мы не знаем, что представляет собой вектор — столбец или строку. Но и это всё не ничего не даёт, т.к. вектор нет необходимости дублировать — пофигк умножению и сложению. И записывать обратно в вектор не можем до окончания операции полностью, т.к.о значение последнего элементо вектора-результата зависит от каждого элемента входного вектора.

yacc> Когда я возился картами рельефа основная просадка была ... на ветвлении. Умножение выполняется быстро.
Какое умножение? На каком процессоре? Мы этого тоже не знаем. И ветвления тоже зависят от архитектуры. Есть несколько водопроводов с предсказанием переходаи или нет?


yacc> Тогда самый быстрый способ - динамическая подлинковка.

Зобудь о ней. Это тоже вызов. Дорогой, с запоминанием точки возравта, со сбросом pipeline (привет ветвлениям).

yacc> самое быстрое вычисление.
Нет.
Самые быстрые — это линейные, когда нет переходов вообще, нет вызовов, линейное и водопровод разогнана, всё в кэше сидит. Это не считая математических преоброзований.

AXT> Сейчас уже не имеет смысла — ещё летом всё было заборото оптимизацией "в лоб", так что интересно только в смысле поиграться с математикой. Матрицы в стандарте ITU-T H.264, если нет доступа к, то могу выложить матрицы.

А, в общем, я ступил, надо было сразу выложить. Умножаем столбец на матрицу.

29	55	74	84
74	74	0	-74
84	-29	-74	55
55	-84	74	-29

Это самая частая, есть ещё:

64	64	64	64
83	36	-36	-83
64	-64	-64	64
36	-83	83	-36

Ну тут структура очевидна ...

PS: Я монстр, я таки обманул тег csv Но только в предварительном просмотре
PPS: Таки есть режим auto

Mishka> И последовательно лежит — это как?
Это с точки зрения адресации процессором.
По столбцам или по строкам - это уже компилятор преобразует [i][j] в соответствующий код

Mishka> И записывать обратно в вектор не можем до окончания операции полностью, т.к.о значение последнего элементо вектора-результата зависит от каждого элемента входного вектора.
В выходной? Это почему? Там операции вида out[k] += in[i] * matix[i][j] так и просятся.
Но быстрее использовать регистровую переменную.

Mishka> Зобудь о ней. Это тоже вызов. Дорогой, с запоминанием точки возравта, со сбросом pipeline (привет ветвлениям).
Ты меня не понял - я имею ввиду загрузку DLL с отображением ф-ии calc в память. А в ф-ии матрица самим набором вычислений представлена. Ее код генерирует отдельная программа-препроцессор.

Mishka> Самые быстрые — это линейные, когда нет переходов вообще, нет вызовов, линейное и водопровод разогнана, всё в кэше сидит. Это не считая математических преоброзований.
Там и нет переходов.
Там операции вида ( псевдо ассемблер, заглавными буквами - КОНСТАНТЫ )
- обнулить регистр С
- загрузить в регистр Б из памяти in + OFFSET1 ( например mov bx, _in + 25 )
- загрузить в регистр А константу MULT1
- умножить регистр А на регистр Б с результатом в А
- сложить регистр А с регистром С с результатов в А
- копировать А в С
...
- копировать С в область памяти out + OFFSET2

и НИКАКИХ ветвлений.

Надо несколько константных матриц? - загрузи несколько DLL вначале в ф-ии calc1, ... calcN
а потом дергай по-усмотрению.
Загрузка делается один раз.
Если вообще особо не меняется - можно и статически подлинковать эти ф-ии в исполняемый файл.

Mishka> Самые быстрые — это линейные, когда нет переходов вообще,

У некоторых DSP для этого имеется мощная фишка — аппаратные циклы. То есть, например (бывают варианты) цикл устроен так (псевдоассемблер):

R1 = (P0)+;
R2 = (P1)+;
LOOP (R0);
R3 = R1 + R2 || R1 = (P0)+;
MNOP || R2 = (P1)+ || (P2)+ = R3;
ENDLOOP;

В цикле автоматом выбираются и исполняются только строчки между LOOP и ENDLOOP, т.е. складываем значения за два такта. Счётчик цикла обрабатывается забесплатно.

yacc> Это с точки зрения адресации процессором.
Копирование займёт время (и опять двойной вложенный цикл). Пробег по этим двум векторам — ещё цикл. И ты говоришь о том, что самое дорогое — это ветвления. Код, написанный ручками для каждого элемента — ни одного ветвления, столько же умножений и сложений с возможным уменьшеием слагаемых (если были нули в матрице), умножения (если были единицы и минус единицы), заменой умножения на сложение для небольших значений (если умножение медленное).

yacc> В выходной? Это почему? Там операции вида out[k] += in[i] * matix[i][j] так и просятся.
В исходный, как было одно время очень популярно на том FORTAN-е.

yacc> Ты меня не понял - я имею ввиду загрузку DLL с отображением ф-ии calc в память. А в ф-ии матрица самим набором вычислений представлена. Ее код генерирует отдельная программа-препроцессор.

Это динамический вызов, который дороже статического. Если ты только не будешь код на лету править, чтобы подставить адреса загруженных процедур.

yacc> Там и нет переходов.
Вызовы. Циклы.

AXT> В цикле автоматом выбираются и исполняются только строчки между LOOP и ENDLOOP, т.е. складываем значения за два такта. Счётчик цикла обрабатывается забесплатно.

Я работал с DSP очень мало, но цикл и там не самая дешёвая вещь. Даже аппаратно поддержанный. Выпрямленный цикл, если память позволяет всё вместить, быстрее.