Я создаю шаблон класса матрицы с плавающей запятой. Объявление класса показано ниже только с соответствующими функциями и членами.
// columns, rows
template <unsigned int c, unsigned int r>
class Matrix {
public:
Matrix(float value);
float& At(unsigned int x, unsigned int y);
float const& At(unsigned int x, unsigned int y) const;
template <unsigned int p> Matrix<p, r> MultipliedBy(Matrix<p, c> const& other);
private:
// column-major ordering
float data_[c][r];
}
Ниже приведены реализации каждой из перечисленных функций.
template <unsigned int c, unsigned int r>
Matrix<c, r>::Matrix(float value) {
std::fill(&data_[0][0], &data_[c][r], value);
}
template <unsigned int c, unsigned int r>
float& Matrix<c, r>::At(unsigned int x, unsigned int y) {
if (x >= c || y >= r) {
return data_[0][0];
}
return data_[x][y];
}
template <unsigned int c, unsigned int r>
float const& Matrix<c, r>::At(unsigned int x, unsigned int y) const {
if (x >= c || y >= r) {
return data_[0][0];
}
return data_[x][y];
}
template <unsigned int c, unsigned int r>
template <unsigned int p>
Matrix<p, r> Matrix<c, r>::MultipliedBy(Matrix<p, c> const& other) {
Matrix<p, r> result(0.0f);
for (unsigned int x = 0; x < c; x++) {
for (unsigned int y = 0; y < r; y++) {
for (unsigned int z = 0; z < p; z++) {
result.At(z, y) += At(x, y) * other.At(z, x);
}
}
}
return result;
}
Теперь несколько строк тестового кода.
Matrix<4, 4> m1;
// m1 set to
//
// 1 2 3 4
// 5 6 7 8
// 9 10 11 12
// 13 14 15 16
Matrix<1, 4> m2;
// m2 set to
//
// 6
// 3
// 8
// 9
Matrix<1, 4> m3 = m1.MultipliedBy(m2);
Здесь все становится странно. При компиляции (с использованием g++
) без оптимизации (-O0
):
// m3 contains
// 0
// 0
// 0
// 0
При любой оптимизации (-O1
, -O2
или -O3
):
// m3 contains
// 210
// 236
// 262
// 288
Обратите внимание, что даже при оптимизации ответ неверный (проверено внешним калькулятором). Поэтому я сузил его до этого вызова в MultipliedBy
:
Matrix<p, r> result(0.0f);
Если я создаю экземпляр result
, то other
становится недействительным (все значения data_
установлены на 0.0f
). До выделения / инициализации result
, other
все еще действует (6, 3, 8, 9
).
Стоит отметить, что если я умножу две матрицы одинакового (квадратного) размера, я получу полностью корректный и правильный вывод, независимо от уровня оптимизации.
Кто-нибудь знает, что в мире g++
тянет сюда? Я запускаю g++ (GCC) 4.6.1
на mingw
... может это как-то связано с проблемой?