Optimized matrix inverse and division code (#149)

This commit is contained in:
Christophe Riccio
2014-01-11 16:44:15 +01:00
parent efdfa577ee
commit 90a249b5ff
9 changed files with 255 additions and 207 deletions

View File

@@ -8,6 +8,11 @@
///////////////////////////////////////////////////////////////////////////////////////////////////
#include <glm/matrix.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtc/ulp.hpp>
#include <vector>
#include <ctime>
#include <cstdio>
using namespace glm;
@@ -175,18 +180,71 @@ int test_inverse()
glm::mat2x2 I2x2 = A2x2 * B2x2;
Failed += I2x2 == glm::mat2x2(1) ? 0 : 1;
return Failed;
}
std::size_t const Count(10000000);
template <typename VEC3, typename MAT4>
int test_inverse_perf(std::size_t Instance, char const * Message)
{
std::vector<MAT4> TestInputs;
TestInputs.resize(Count);
std::vector<MAT4> TestOutputs;
TestOutputs.resize(TestInputs.size());
VEC3 Axis(glm::normalize(VEC3(1.0f, 2.0f, 3.0f)));
for(std::size_t i = 0; i < TestInputs.size(); ++i)
{
typename MAT4::value_type f = static_cast<typename MAT4::value_type>(i + Instance) * typename MAT4::value_type(0.1) + typename MAT4::value_type(0.1);
TestInputs[i] = glm::rotate(glm::translate(MAT4(1), Axis * f), f, Axis);
//TestInputs[i] = glm::translate(MAT4(1), Axis * f);
}
std::clock_t StartTime = std::clock();
for(std::size_t i = 0; i < TestInputs.size(); ++i)
TestOutputs[i] = glm::inverse(TestInputs[i]);
std::clock_t EndTime = std::clock();
for(std::size_t i = 0; i < TestInputs.size(); ++i)
TestOutputs[i] = TestOutputs[i] * TestInputs[i];
typename MAT4::value_type Diff(0);
for(std::size_t Entry = 0; Entry < TestOutputs.size(); ++Entry)
{
MAT4 i(1.0);
MAT4 m(TestOutputs[Entry]);
for(glm::length_t y = 0; y < m.length(); ++y)
for(glm::length_t x = 0; x < m[y].length(); ++x)
Diff = glm::max(m[y][x], i[y][x]);
}
//glm::uint Ulp = 0;
//Ulp = glm::max(glm::float_distance(*Dst, *Src), Ulp);
printf("inverse<%s>(%f): %d\n", Message, Diff, EndTime - StartTime);
return 0;
};
int main()
{
int Failed = 0;
Failed += test_matrixCompMult();
Failed += test_outerProduct();
Failed += test_transpose();
Failed += test_determinant();
Failed += test_inverse();
return Failed;
int Error(0);
Error += test_matrixCompMult();
Error += test_outerProduct();
Error += test_transpose();
Error += test_determinant();
Error += test_inverse();
for(std::size_t i = 0; i < 1; ++i)
{
Error += test_inverse_perf<glm::vec3, glm::mat4>(i, "mat4");
Error += test_inverse_perf<glm::dvec3, glm::dmat4>(i, "dmat4");
}
return Error;
}

View File

@@ -16,19 +16,19 @@
void print(glm::dmat4 const & Mat0)
{
printf("mat4(\n");
printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[0][0], Mat0[0][1], Mat0[0][2], Mat0[0][3]);
printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[1][0], Mat0[1][1], Mat0[1][2], Mat0[1][3]);
printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[2][0], Mat0[2][1], Mat0[2][2], Mat0[2][3]);
printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f))\n\n", Mat0[3][0], Mat0[3][1], Mat0[3][2], Mat0[3][3]);
printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[0][0], Mat0[0][1], Mat0[0][2], Mat0[0][3]);
printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[1][0], Mat0[1][1], Mat0[1][2], Mat0[1][3]);
printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[2][0], Mat0[2][1], Mat0[2][2], Mat0[2][3]);
printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f))\n\n", Mat0[3][0], Mat0[3][1], Mat0[3][2], Mat0[3][3]);
}
void print(glm::mat4 const & Mat0)
{
printf("mat4(\n");
printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[0][0], Mat0[0][1], Mat0[0][2], Mat0[0][3]);
printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[1][0], Mat0[1][1], Mat0[1][2], Mat0[1][3]);
printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f)\n", Mat0[2][0], Mat0[2][1], Mat0[2][2], Mat0[2][3]);
printf("\tvec4(%2.3f, %2.3f, %2.3f, %2.3f))\n\n", Mat0[3][0], Mat0[3][1], Mat0[3][2], Mat0[3][3]);
printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[0][0], Mat0[0][1], Mat0[0][2], Mat0[0][3]);
printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[1][0], Mat0[1][1], Mat0[1][2], Mat0[1][3]);
printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f)\n", Mat0[2][0], Mat0[2][1], Mat0[2][2], Mat0[2][3]);
printf("\tvec4(%2.9f, %2.9f, %2.9f, %2.9f))\n\n", Mat0[3][0], Mat0[3][1], Mat0[3][2], Mat0[3][3]);
}
int test_inverse_mat4x4()
@@ -107,6 +107,66 @@ int test_inverse()
Error += glm::all(glm::epsilonEqual(Identity[3], glm::vec4(0.0f, 0.0f, 0.0f, 1.0f), glm::vec4(0.01f))) ? 0 : 1;
}
{
glm::highp_mat4 const Matrix(
glm::highp_vec4(0.6f, 0.2f, 0.3f, 0.4f),
glm::highp_vec4(0.2f, 0.7f, 0.5f, 0.3f),
glm::highp_vec4(0.3f, 0.5f, 0.7f, 0.2f),
glm::highp_vec4(0.4f, 0.3f, 0.2f, 0.6f));
glm::highp_mat4 const Inverse = glm::inverse(Matrix);
glm::highp_mat4 const Identity = Matrix * Inverse;
printf("highp_mat4 inverse\n");
print(Matrix);
print(Inverse);
print(Identity);
Error += glm::all(glm::epsilonEqual(Identity[0], glm::highp_vec4(1.0f, 0.0f, 0.0f, 0.0f), glm::highp_vec4(0.01f))) ? 0 : 1;
Error += glm::all(glm::epsilonEqual(Identity[1], glm::highp_vec4(0.0f, 1.0f, 0.0f, 0.0f), glm::highp_vec4(0.01f))) ? 0 : 1;
Error += glm::all(glm::epsilonEqual(Identity[2], glm::highp_vec4(0.0f, 0.0f, 1.0f, 0.0f), glm::highp_vec4(0.01f))) ? 0 : 1;
Error += glm::all(glm::epsilonEqual(Identity[3], glm::highp_vec4(0.0f, 0.0f, 0.0f, 1.0f), glm::highp_vec4(0.01f))) ? 0 : 1;
}
{
glm::mediump_mat4 const Matrix(
glm::mediump_vec4(0.6f, 0.2f, 0.3f, 0.4f),
glm::mediump_vec4(0.2f, 0.7f, 0.5f, 0.3f),
glm::mediump_vec4(0.3f, 0.5f, 0.7f, 0.2f),
glm::mediump_vec4(0.4f, 0.3f, 0.2f, 0.6f));
glm::mediump_mat4 const Inverse = glm::inverse(Matrix);
glm::mediump_mat4 const Identity = Matrix * Inverse;
printf("mediump_mat4 inverse\n");
print(Matrix);
print(Inverse);
print(Identity);
Error += glm::all(glm::epsilonEqual(Identity[0], glm::mediump_vec4(1.0f, 0.0f, 0.0f, 0.0f), glm::mediump_vec4(0.01f))) ? 0 : 1;
Error += glm::all(glm::epsilonEqual(Identity[1], glm::mediump_vec4(0.0f, 1.0f, 0.0f, 0.0f), glm::mediump_vec4(0.01f))) ? 0 : 1;
Error += glm::all(glm::epsilonEqual(Identity[2], glm::mediump_vec4(0.0f, 0.0f, 1.0f, 0.0f), glm::mediump_vec4(0.01f))) ? 0 : 1;
Error += glm::all(glm::epsilonEqual(Identity[3], glm::mediump_vec4(0.0f, 0.0f, 0.0f, 1.0f), glm::mediump_vec4(0.01f))) ? 0 : 1;
}
{
glm::lowp_mat4 const Matrix(
glm::lowp_vec4(0.6f, 0.2f, 0.3f, 0.4f),
glm::lowp_vec4(0.2f, 0.7f, 0.5f, 0.3f),
glm::lowp_vec4(0.3f, 0.5f, 0.7f, 0.2f),
glm::lowp_vec4(0.4f, 0.3f, 0.2f, 0.6f));
glm::lowp_mat4 const Inverse = glm::inverse(Matrix);
glm::lowp_mat4 const Identity = Matrix * Inverse;
printf("lowp_mat4 inverse\n");
print(Matrix);
print(Inverse);
print(Identity);
Error += glm::all(glm::epsilonEqual(Identity[0], glm::lowp_vec4(1.0f, 0.0f, 0.0f, 0.0f), glm::lowp_vec4(0.01f))) ? 0 : 1;
Error += glm::all(glm::epsilonEqual(Identity[1], glm::lowp_vec4(0.0f, 1.0f, 0.0f, 0.0f), glm::lowp_vec4(0.01f))) ? 0 : 1;
Error += glm::all(glm::epsilonEqual(Identity[2], glm::lowp_vec4(0.0f, 0.0f, 1.0f, 0.0f), glm::lowp_vec4(0.01f))) ? 0 : 1;
Error += glm::all(glm::epsilonEqual(Identity[3], glm::lowp_vec4(0.0f, 0.0f, 0.0f, 1.0f), glm::lowp_vec4(0.01f))) ? 0 : 1;
}
{
glm::mat4 const Matrix(
glm::vec4(0.6f, 0.2f, 0.3f, 0.4f),