Короткий курс комп'ютерної графіки: пишемо спрощений OpenGL своїми руками, стаття 3.14 з 6

Зміст основного курсу
Поліпшення коду

4 Привітання і вступ

Нумерація у минулій статті закінчилася на 3, в цій будемо продовжувати нумерувати наскрізь.
UPD: УВАГА! Розділ, починаючи з номера 3.1, 3.14 і 3.141 і далі, буде про тонкощі реалізації основи основ комп'ютерної графіки — лінійної алгебри та обчислювальної геометрії. Про засади графіки пише haqreu, я ж буду писати про те, як це можна виразно запрограмувати!

Ця стаття є продовженням серії статей про практичної реалізації елементів обчислювальної геометрії, і, зокрема, програмного отрисовщика, з використанням C++98. Ми з haqreu свідомо йдемо на використання попередньої версії стандарту і написання власної геометричної бібліотеки для того, щоб, по-перше, випустити код прикладів, які без особливих труднощів будуть компілюватися більшістю наявних компіляторів, а по-друге, щоб в нашому коді не було нічого, що приховано в надрах бібліотеки. У статті викладено питання реалізації шаблону прямокутної матриці
template<size_t DimRows,size_t DimCols,typename number_t> class mat;


4.1 Подяки
Я висловлюю величезну вдячність haqreu, як основоположнику даного курсу. Так тримати!
Я дуже вдячний lemelisk за попереднє рецензування і рев'ю моїх джерел. Дякуємо за плідні дискусії!

4.2 Невелике філософське зауваження про мову, математику, C++, і всьому такому
В цілому, такі абстракції, як «матриця» і «вектор» та операції над ними, необхідні для того, щоб спростити і формалізувати мову викладу геометрії. Без них нам доводилося б постійно говорити (й писати в програмі) щось на кшталт:

Говоримо:
Щоб отримати нові координати точки Point, яка переміщена в напрямку Direction, потрібно скласти
x точки Point з x переміщення Direction,
біля точки Point з y переміщення Direction,
z точки Point з z переміщення Direction
Пишемо:
Point.x=Point.x+Direction.x;
Point.y=Point.y+Direction.y;
Point.z=Point.x+Direction.z; //(Привіт, "Ефект останнього рядка")

Після того, як поняття «вектор» увійшло в наш мовний зворот, ми вже можемо говорити:
Щоб отримати новий радіус-вектор точки Point після переміщення в напрямку Direction, потрібно скласти вектори Point і Direction
А після включення в наш код шаблону vec, ми можемо писати:
Point=Point+Direction;
Таким чином, ми скорочуємо написане, страхуємося від помилок при копіювання-вставки-правці і отримуємо спосіб міркувати більш ємними термінами.

Ось що про це пише Джефф Елджер у книзі «C++ for Real programmers» (видавництво AP Professional, переклад випущений видавництвом «Пітер» під назвою «C++: Бібліотека програміста»)

С++ -це насправді не стільки мова, скільки інструмент для створення ваших власних мов. Його елегантність полягає аж ніяк не в простоті (слова С++ і простота ріжуть слух своїм явним протиріччям), а в його потенційні можливості. За кожною потворною проблемою ховається якась розумна ідіома, витончений мовної фінт, завдяки якому проблема тане прямо на очах. Проблема вирішується так само елегантно, як це зробив би справжній мову типу Smalltalk або Lisp, але при цьому ваш процесор не димить від напруги, а на Уолл-Стріт не ростуть акції виробників чіпів пам'яті
Саме цим ми і займемося — будемо будувати з використанням C++ мова керування геометричними примітивами, намагаючись робити це якомога ближче до визначень з математики.

5 Шаблон класу для обробки і зберігання матриць

Ми зберігаємо матрицю як масив вектор-рядків, використовуючи відповідну спеціалізацію шаблону vec. Наш шаблон орієнтований тільки на роботу з матрицями малої розмірності (максимум 10x10 елементів. Хоча це вже буде жах). Але в обчислювальній геометрії, як правило, вся робота йде саме з матрицями не більше, ніж 4x4. Великі ж матриці виходять, як правило, дуже розрідженими — для них заздалегідь відомо, що вони практично повністю заповнені нулями, що дозволяє оптимізувати їх обробку і зберігання.
5.1 Індексні оператори
Ми визначили константный і неконстантный варіанти індексного оператора [] таким чином, щоб вони повертали посилання на відповідну вектор-рядок матриці:
вихідні тексти
Константный варіант
Неконстантный варіант
const vec<DimCols,number_t>& operator[] (const size_t idx) const 
{
assert(idx<DimRows);
return rows[idx];
}
vec<DimCols,number_t>& operator[] (const size_t idx) 
{
assert(idx<DimRows);
return rows[idx];
}
Два варіанти операторів потрібні із-за того, що ми скрізь, де тільки можливо, застосовуємо
const
. Якщо б не було константного варіанти індексного оператора, компілятор не дозволяв би нам зібрати код, в якому, наприклад, оператор викликається з константною посиланням на примірник матриці:
Так лається GCC 4.9.2, коли не знаходить константного варіанти оператора []
geometry.h:182: помилка: no match for 'operator[]' (operand types are 'const mat<2ul, 3ul, float>' and 'size_t {aka long unsigned int}')
for (size_t i=DimRows; i--; ret[i]=lhs[i]*rhs);
^
Крім того, ми тут застосовуємо макрос
assert()
для відлову спроб звернутися до неіснуючої рядку матриці. Докладніше про assert написано в статті 3.1

5.2 Обчислення оберненої матриці
Нам потрібна була підпрограма для обчислення оберненої матриці. Зокрема, використовуючи обернену матрицю можна «красиво» знаходити барицентричні координати (можна подивитися в моїй гілці на гітхабі, функція
filltria
).
Задача обчислення оберненої матриці в загальному випадку є досить витратним за часом і може супроводжуватися накопиченням значної похибки при обчисленнях. Однак для таких малих матриць ці проблеми не так істотні (якщо не будувати такий приклад спеціально). Ми вибрали алгоритм з застосуванням союзної матриці, тому що він здався нам більш прозорим, ніж алгоритм Гауса-Жордана. Використовується наступний факт:
  • якщо для даної матриці скласти союзну матрицю,
  • а потім транспонувати
  • і поділити її на визначник даної,
ми отримаємо обернену до матриці. Це всього три виразних кроку.
Союзна ж матриця складається з алгебраїчних доповнень. Алгебраїчне доповнення — це визначник підматриці, отриманої з основним шляхом видалення рядків і стовпців, що містять елемент, для якого ми знаходимо це додаток.
Пояснювальна картинка з Вікіпедії, автор — Олександр Мехоношинimage

Отримали наступне питання — як знайти визначник.
Використовуємо таку властивість:
визначник дорівнює добутку елементів рядка або стовпця матриці на їх алгебраїчні доповнення.
Вийшла рекурсія — алгебраїчне доповнення — теж визначник матриці, тільки розмірність на одиницю менше. А визначник матриці розмірності 1 — це її єдиний елемент. Рекурсія замкнулася. Ми скористаємося чудовою особливістю шаблонів C++ — розгорнемо цю рекурсію прямо під час компіляції.

Зауваження: Матриця у нас — прямокутна, а визначник можна обчислити тільки для квадратної матриці. Хотілося б вже на етапі компіляції відловлювати спроби знайти визначник неквадратной матриці. Є два способи це зробити:
  • успадкувати від mat клас squareMat, в якому оголосити визначник
  • скористатися тим, що розміри матриці нам відомі ще на етапі компіляції, тому ми можемо порівняти їх і перервати компіляцію, якщо вони не збіглися
У першому випадку ми отримаємо окремий тип матриць і купу проблем з перетворенням матриць туди-сюди, бо як у нас може утворитися квадратна матриця в рамках типу mat, але щоб обчислити визначник, нам доведеться конструювати з неї екземпляр типу squareMat.
Тому ми виберемо другий спосіб — при розгортанні шаблонів виявляти розбіжність розмірів і ламати компіляцію. В С++11 є стандартне засіб це зробити — static_assert(). В С++98 такого засобу немає, тому ми винайдемо сурогатне:
template<bool>
struct static_assertion;

template<>
struct static_assertion<true> {};
Якщо в цей шаблон підставляти true, просто утворюється порожня структура static_assertion{}. А от якщо поставити false, утворюється помилка (спроба використовувати неповне оголошення структури static_assertion), що і зупиняє компіляцію:
geometry.h:125: помилка: invalid use of incomplete type 'struct static_assertion<false>'
static_assertion<DimCols<=10>();
^
Чим більше помилок зміщений на етап компіляції, тим краще, адже вони просто не допускають поганий код до роботи.
Повернемося до нашої рекурсії:
Вихідні тексти шаблону, який відповідає за рекурсію
///////////////////////////////////////////// повний шаблон структури dt
1. template<size_t Dim,typename number_t> struct dt 
2. {
3. static number_t det(const mat<Dim,Dim,number_t>& src) 
4. {
5. number_t ret=0;
6. for (size_t i=Dim; i--;)
7. {
8. ret += src[0][i] * src.cofactor(0,i);
9. }
10. return ret;
11. }
12. };
///////////////////////////////////////////// часткова спеціалізація шаблону для випадку, коли розмірність дорівнює 1.
13. template < typename number_t> struct dt<1,number_t> 
14. {
15. static number_t det(const mat<1,1,number_t>& src) 
16. {
17. return src[0][0];
18. }
19. };


Отже, для всіх ситуацій, коли перший параметр шаблону не дорівнює 1, буде застосовуватися верхній шаблон (рядки 1-11). Для ситуації, перший параметр дорівнює 1, визначена спеціальна версія шаблону (рядки 12-18). Зверніть увагу на рядок 13: Доповнення
<1,number_t>
якраз і вказує, що це оголошення — окрема сутність. З загальної версією шаблону dt її нічого не пов'язує. Зустрічається помилка, коли плутають поведінка часткової спеціалізації і спадкування, чекаючи, що все оголошене в загальній версії виявиться і в частковій. Цього не буде — на програміста покладається відтворення часткової спеціалізації все, що йому потрібно.
  • Обернення в окрему структуру виконано тому, що стандарт не дозволяє робити часткову спеціалізацію функцій.
  • Ми винесли цю функцію з mat тому, що вступивши інакше, були б змушені дописувати в mat<1,1,number_t> частина методів із загальної версії, що призвело б до деякого раздутию коду.
Часткова спеціалізація, точно за визначенням визначника матиці 1x1, віддає нам єдиний елемент (рядок 16).

Загальна ж версія шаблону виконує розкладання матриці за нульовою рядку — вона бере елемент з індексом
0, i
з першого рядка, множить його на алгебраїчне доповнення —
cofactor(0,i)
а результат накопичує у змінній
ret
. Все в повній згоді з формулою цього розкладання:image

5.3 Як відбувається рекурсія
ВАЖЛИВО: Строго кажучи, механізм обробки шаблонів — особиста справа компілятора. Нижче я розповідаю про те, як це могло б відбуватися. Це робиться тільки для ілюстрації того, як рекурсивно розгортаються декларації шаблонів.

Уявімо, що ми самостійно виконуємо розгортання шаблонів. Написаний код:
...
using namespace std;
int main()
{
mat<3,3,int> m;
....
cout << m.det() << endl;
return 0;
}
Ми зустрічаємо в ньому виклик det() з шаблону mat, а отже, повинні отримати текст функції det() для конкретного випадку mat<3,3,int>.
Шаблон
Результат підстановки DimCols=3; DimRows=3; number_t=int
number_t det() const 
{
static_assertion<DimCols==DimRows>();
return dt<DimCols,number_t>::det(*this);
}
int det() const<3,int> 
{
static_assertion<3==3>();
return dt<3,int>::det(*this);
}
Відмінно, тепер нам потрібно з шаблонів template<size_t Dim,typename number_t> struct dt отримати конкретний текст для випадку dt<3,int>:
Шаблон
Результат підстановки Dim=3; number_t=int
template<size_t Dim,typename number_t> struct dt 
{
static number_t det(const mat<Dim,Dim,number_t>& src) 
{
number_t ret=0;
for (size_t i=Dim; i--;)
{
ret += src[0][i] * src.cofactor(0,i);
}
return ret;
}
};
struct dt<3,int>
{
static int det(const mat<3,3,int>& src) 
{
int ret=0;
for (size_t i=3; i--;)
{
ret += src[0][i] * src.cofactor(0,i);
}
return ret;
}
};
Ми зустріли виклик оператора [] і cofactor() з mat<3,3,int>. Тіло [] зараз для нас інтересу не представляє, розбираємося з cofactor():
Шаблон
Результат підстановки DimCols=3; DimRows=3; number_t=int
number_t cofactor(size_t row, size_t col) const 
{
return get_minor(row,col).det()*((row+col)%2 ? -1 : 1);
}
int cofactor(size_t row, size_t col) const 
{
return get_minor(row,col).det() *( (row + col) % 2 ? -1 : 1);
}
Тут викликається get_minor(). Отримаємо для неї текст:
Шаблон
Результат підстановки DimCols=3; DimRows=3; number_t=int
mat<DimRows-1,DimCols-1,number_t> get_minor(size_t row, size_t col) const 
{
mat<DimRows-1,DimCols-1,number_t> ret;
for (size_t i=DimRows-1; i--; )
{
for (size_t j=DimCols-1;j--;)
{
ret[i][j]=rows[ i<row ? i : i+1 ] [ j <col ? j : j+1 ]; 
// конструкція j <col ? j : j+1 забезпечує пропуск стовпця з номером col
//аналогічно i<row ? i : i+1 - забезпечує пропуск рядка з номером row
} 
}
return ret;
}
mat<2,2,int> get_minor(size_t row, size_t col) const 
{
mat<2,2,int> ret;
for (size_t i=2; i--; )
{
for (size_t j=2;j--;)
{
ret[i][j]=rows[ i<row?i:i+1 ] [ j<col?j:j+1 ];
}
}
return ret;
}
Ми отримали, що цей get_minor() повертає вже матрицю 2х2. Тепер ми, для матриці 2x2 будемо викликати mat<2,2,int>::det(). Але у нас є тільки тіло mat<3,3,int>::det(). Тому ми робимо заглиблення в рекурсію на один крок і будуємо ще один набір методів:
Методи, що з'явилися при розгортанні шаблону 2x2Ми повинні отримати текст функції det() для конкретного випадку mat<2,2,int>:
Шаблон
Результат підстановки Dim=2; number_t=int
number_t det() const 
{
static_assertion<DimCols==DimRows>();
return dt<DimCols,number_t>::det(*this);
}
int det() const<2,int> 
{
static_assertion<2==2>();
return dt<2,int>::det(*this);
}
Відмінно, тепер нам потрібно з шаблонів template<size_t Dim,typename number_t> struct dt отримати конкретний текст для випадку dt<2,int>:
Шаблон
Результат підстановки Dim=2; number_t=int
template<size_t Dim,typename number_t> struct dt 
{
static number_t det(const mat<Dim,Dim,number_t>& src) 
{
number_t ret=0;
for (size_t i=Dim; i--;)
{
ret += src[0][i] * src.cofactor(0,i);
}
return ret;
}
};
struct dt<2,int>
{
static int det(const mat<2,2,int>& src) 
{
int ret=0;
for (size_t i=2; i--;)
{
ret += src[0][i] * src.cofactor(0,i);
}
return ret;
}
};
Ми зустріли виклик оператора [] і cofactor() з mat<2,2,int>. Тіло [] зараз для нас інтересу не представляє, розбираємося з cofactor():
Шаблон
Результат підстановки DimCols=2; DimRows=2; number_t=int
number_t cofactor(size_t row, size_t col) const 
{
return get_minor(row,col).det()*((row+col)%2 ? -1 : 1);
}
int cofactor(size_t row, size_t col) const 
{
return get_minor(row,col).det() *( (row + col) % 2 ? -1 : 1);
}
Тут викликається get_minor(). Отримаємо для неї текст:
Шаблон
Результат підстановки DimCols=2; DimRows=2; number_t=int
mat<DimRows-1,DimCols-1,number_t> get_minor(size_t row, size_t col) const 
{
mat<DimRows-1,DimCols-1,number_t> ret;
for (size_t i=DimRows-1; i--; )
{
for (size_t j=DimCols-1;j--;)
{
ret[i][j]=rows[ i<row?i:i+1 ] [ j<col?j:j+1 ];
}
}
return ret;
}
mat<1,1,int> get_minor(size_t row, size_t col) const 
{
mat<1,1,int> ret;
for (size_t i=1; i--; )
{
for (size_t j=1;j--;)
{
ret[i][j]=rows[ i<row?i:i+1 ] [ j<col?j:j+1 ];
}
}
return ret;
}

Розгорнувши ланцюжок шаблонів, ми знову зіштовхнулися з викликом get_minor() в тілі cofactor(), але вже для випадку DimCols=2, DimRows=2. Отримана версія get_minor() повертає матрицю 1x1. Далі, для неї в тілі cofactor() запитується визначник. Ми знову повинні отримати текст mat::det(), але вже для випадку DimCols=1; DimRows=1.
Шаблон
Результат підстановки DimCols=1; DimRows=1; number_t=int
number_t det() const 
{
static_assertion<DimCols==DimRows>();
return dt<DimCols,number_t>::det(*this);
}
int det() const<1,int> 
{
static_assertion<1==1>();
return dt<1,int>::det(*this);
}
Відмінно, тепер нам потрібно з шаблонів template<size_t Dim,typename number_t> struct dt отримати конкретний текст для випадку dt<1,int>. СТОП, це ж потрапляє під часткову спеціалізацію шаблону dt<1,number_t>! Потрібно обробляти саме її, а не загальну версію шаблону!
Шаблон
Результат підстановки Dim=1; number_t=int
template < typename number_t> struct dt<1,number_t> 
{
static int det(const mat<1,1,number_t>& src) 
{
return src[0][0];
}
};
template < typename number_t> struct dt<1,int> 
{
static number_t det(const mat<1,1,int>& src) 
{
return src[0][0];
}
};
Бачимо, що тут вже немає звернення до cofactor(). На цьому розкрутка рекурсії закінчена — ми отримали тіла функцій для обчислення визначника матриці 3x3, 2x2, 1x1.
Важливе зауваження: Може виникнути бажання написати так:
number_t det()
{
if(DimCols==1 && DimRows==1)
{
return rows[0][0];
}
number_t ret=0;
for (size_t i=DimCols; i--;)
{
ret += src[0][i] * src.cofactor(0,i);
}
return ret;
}

Мотивація проста — це ніби як звичайна рекурсія, так і часткова спеціалізація не потрібна. Однак під час розгортання шаблонів, компілятор нічого не буде робити з умовним оператором, а почне будувати шаблони(3, 2, 1, 0, std::numeric_limits<size_t>::max, std::numeric_limits<size_t>::max-1, ...), поки його не зупинить внутрішній лічильник
ftemplate-depth=
(за замовчуванням ця глибина занурення дорівнює 900).
Можливість використовувати звичні оператори управління потоком з'являється в C++11 для функцій, оголошених c ключовим словом
constexpr
. Вони дозволяють організувати ще один вид обчислень під час компіляції, відмінний від рекурсії на шаблонах.

Отже, ми описали в термінах C++ відшукання визначника матриці, а також супутні операції — побудова мінора матриці та алгебраїчні доповнення. Для відшукання оберненої матриці залишилося тільки побудувати союзну матрицю:
Метод, який будує союзну матрицю:
mat<DimRows,DimCols,number_t> adjugate() const 
{
mat<DimRows,DimCols,number_t> ret;
for (size_t i=DimRows; i--;)
{
for (size_t j=DimCols; j--;)
{
ret[i][j]=cofactor(i,j)
}
} 
return ret;
}
діє точно за визначенням — заповнює матрицю алгебраїчними доповненнями відповідних елементів.
Залишилося виконати головний етап побудов — знайти зворотну матрицю (вона буде отримана транспонированном вигляді):
mat<DimRows,DimCols,number_t> invert_transpose() 
{
mat<DimRows,DimCols,number_t> ret = adjugate(); //отримали союзну матрицю

number_t det = ret[0]*rows[0]; //вирахували визначник вихідної матриці, шляхом скалярного
//множення першої вектор-рядка
//елементів на вектор рядок їх алгебраїчних доповнень

return ret/det; //поділили союзну матрицю на визначник вихідної,
//отримавши тим самим транспоновану обернену матрицю
}

5.4 Множення матриць
Операція множення матриць також досить затратна по часу і може супроводжуватися накопиченням похибки. У випадках матрицями більшої розмірності (більше, ніж 64x64), як правило, застосовують спеціальні алгоритми, на зразок алгоритму Штрассена. Асимптотичний рекордсмен — Алгоритм Копперсмита-Винограду (з поліпшеннями Вірджинії Вільямс), зараз не застосовується, тому як починає обганяти алгоритм Штрассена на матрицях настільки величезні, що в сучасній обчислювальній практиці ще не зустрічаються. У нас матриці дрібні, тому нам для множення підходить просто визначення операції множення матриць:
Добуток матриць AB складається з усіх можливих комбінацій скалярних творів вектор-рядків матриці A і вектор-стовпців матриці B
image
Визначення в такому вигляді дозволяє нам використовувати операцію скалярного множення з нашого шаблону vec. Але для цього нам потрібен метод, извлекающий з матриці її стовпець у вигляді вектора:
Прихований текст
vec<DimRows,number_t> col(const size_t idx) const 
{
assert(idx<DimCols);
vec<DimRows,number_t> ret;
for (size_t i=DimRows; i--;)
{
ret[i]=rows[i][idx]; 
}
return ret;
}
Так, тут нібито має місце бути зайве копіювання (пізніше ми подивимося, чи впорається компілятор з його викиданням), але головна мета — написати таке коротке і ємне множення матриць:
template<size_t R1,size_t C1,size_t C2,typename number_t>mat<R1,C2,number_t> operator*(const mat<R1,C1,number_t>& lhs, const mat<C1,C2,number_t>& rhs) 
{
mat<R1,C2,number_t> result;
for (size_t i=R1; i--; )
{
for (size_t j=C2; j--;)
{
result[i][j]=lhs[i]*rhs.col(j); //перемножити скалярно i-тий рядок лівого операнда
// на j-тий стовпець правого
}
}
return result;
}

Приділимо увагу до цього тексту. Тут на C++ викладено математичне визначення операції множення прямокутних матриць, причому зроблено це дуже близько до тексту[«заповнення матриці-результату всілякими скалярними добутками вектор-рядків лівої матриці на вектор-стовпці правою»].
При цьому контроль розмірностей (число стовпців лівої має бути одно числу рядків правої) виконується відразу на етапі компіляції — спроба перемножити не підходять за розміром матриці навіть не откомпилируется.
А це дуже важливо, адже помилка часу виконання може бути плаваючою, тою, що важко вхопити і так далі. А помилка компіляції вилазить практично завжди. (Треба врахувати божевілля IDE, яка може частину проекту закешувати, частина забути оновити і так далі. Але від цього є вірне ліки — Rebuild All).
5.5 Генератор одиничних матриць
В його тексті:
static mat<DimRows,DimCols,number_t> identity() 
{
mat<DimRows,DimCols,number_t> ret;
for (size_t i=DimRows; i--; )
{
for (size_t j=DimCols;j--;)
{
ret[i][j]=(i==j) ;
}
}
return ret;
}

Є цікава особливість — запис результату типу bool в змінну апріорно невідомого типу number_t. Тут може трапитися кошмарик у людини, яка захотіла б використовувати нашу геометрію в своєму проекті, причому разом з саморобним типом для зберігання чисел (це може бути тип, який зберігає числа з фіксованою комою). Він може не чекати, що наші методи будуть намагатися записувати в його тип значення bool (оператор == для size_t повертає bool), і нагородити власний конструктор, bool приймає і робить це перетворення не так, як описано в стандарті ([§4.7/4] вимагає, щоб при приведенні bool у числовий тип, false перетворювалося в 0, true — 1). Так як в даному випадку, від стандарту відступає гіпотетичний користувач, нам не слід намагатися підстилати йому соломку (ми могли б створити таке:
ret[i][j]=static_cast<number_t>(static_cast<int>(i==j))
, і передавати йому хоча б int) ціною низки перетворень — порушив стандарт — сам винен.

5.6 Множення матриці на вектор
Для множення матриці (ліворуч) на вектор (праворуч), ми розглянемо вектор як матрицю з одного стовпця. Далі, ми можемо виконати вже відому нам операцію множення матриць, яка в підсумку також дасть матрицю з одного стовпця — цей вектор-стовпець і буде результатом множення.
Нам знадобиться метод set_col() для заповнення стовпця матриці з вектора:
void set_col(size_t idx, vec<DimRows,number_t> v) 
{
assert(idx<DimCols);
for (size_t i=DimRows; i--;)
{
rows[i][idx]=v[i];
}
}

генератор матриці по заданому вектор-стовпцем
static mat<DimRows,1,number_t> make_from_vec_col(const vec<DimRows,number_t>& v)
{
static_assertion<DimCols==1>();
mat<DimRows,1,number_t> ret;
ret.set_col(0,v);
return ret;
}

Тепер ми можемо записати оператор множення матриці на вектор:
template<size_t DimRows,size_t DimCols,typename number_t> vec<DimRows,number_t> operator*(const mat<DimRows,DimCols,number_t>& lhs, const vec<DimCols,number_t>& rhs)
{
return (lhs * mat<DimCols,1,number_t>::make_from_vec_col(rhs) ).col(0);
}
Спрацьовано точно за визначенням отримали з вектора матрицю з одним стовпцем, перемножили, повернули стовпець-результат.
5.7 Розподіл матриці на скаляр, висновок ostream
Виконуються з використанням того факту, що матриця у нас зберігається як масив вектор — рядків, а для вектора аналогічні операції вже написані:
Тексти операцій
template<size_t DimRows,size_t DimCols,typename number_t>mat<DimCols,DimRows,number_t> operator/(mat<DimRows,DimCols,number_t> lhs, const number_t& rhs) 
{
for (size_t i=DimRows; i--; )
{
lhs[i]=lhs[i]/rhs;
}
return lhs;
}

template <size_t DimRows,size_t DimCols,class number_t> std::ostream& operator<<(std::ostream& out, mat<DimRows,DimCols,number_t>& m) 
{
for (size_t i=0; i<DimRows; i++) 
{
out << m[i] << std::endl;
}
return out;
}


6 Висновок

У цій статті я постарався максимально доступно описати побудову шаблону класу mat. У наступній статті я торкнуся питання оптимізації та прискорення даного коду на архітектурах x86_64 і ARM (на прикладі STM32F3). До швидкої статті!

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.