math: Mat4f tests: Add test_math_mat4f_01 (simple); test_math_mat4f_02_perf (perfomance)

Native C++ Mat4f multiplications beat JOGL's Matrix4f Java implementation by factor >13 (Java-17/C++ -O3),
utilizing same code logic.

Java:
Summary loops 25000000: I4a   1080 ms total, 0.021600 us/mul, I4a / I2 141.732283%, I4a / I4b 99.447514%
Summary loops 25000000: I4b   1086 ms total, 0.021720 us/mul, I4b / I2 142.519685%, I4b / I4a 100.555556%

C++:
Summary loops 300000000: I4a    981 ms total (981,836 us), 1.636394 ns/mul, I4a / I4b 99.785742%
Summary loops 300000000: I4b    983 ms total (983,944 us), 1.639908 ns/mul, I4b / I4a 100.214718%

Java: 21.600 ns/mul <- 0.021600 us/mul
C++:   1.64 ns/mul

Java / C++: >13

3 files changed, 333 insertions, 0 deletions

diff --git a/include/jau/math/mat4f.hpp b/include/jau/math/mat4f.hpp
index f3c51f1..3e25fe5 100644
--- a/include/jau/math/mat4f.hpp
+++ b/include/jau/math/mat4f.hpp
@@ -693,6 +693,7 @@ class alignas(Value_type) Matrix4 {
     }
 
   private:
+    /** Returns the maximum abs(mxy) field */
     value_type absMax() const noexcept {
         value_type max = std::abs(m00);
         max = std::max(max, std::abs(m01));
diff --git a/test/test_math_mat4f_01.cpp b/test/test_math_mat4f_01.cpp
new file mode 100644
index 0000000..74dec57
--- /dev/null
+++ b/test/test_math_mat4f_01.cpp
@@ -0,0 +1,162 @@
+/*
+ * Author: Sven Gothel <[email protected]>
+ * Copyright (c) 2024 Gothel Software e.K.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be
+ * included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+#include <thread>
+#include <cassert>
+#include <cinttypes>
+#include <cstring>
+
+#include <jau/test/catch2_ext.hpp>
+
+#include <jau/int_math.hpp>
+#include <jau/float_math.hpp>
+#include <jau/math/vec2f.hpp>
+#include <jau/math/vec2i.hpp>
+#include <jau/math/vec3f.hpp>
+#include <jau/math/vec4f.hpp>
+#include <jau/math/mat4f.hpp>
+#include <jau/math/quaternion.hpp>
+#include <jau/math/aabbox2f.hpp>
+#include <jau/math/aabbox3f.hpp>
+#include <jau/math/mat4f.hpp>
+#include <jau/math/recti.hpp>
+#include <jau/math/math_error.hpp>
+
+using namespace jau;
+using namespace jau::math;
+
+static const float EPSILON = std::numeric_limits<float>::epsilon();
+
+static const float mI_0[] = {    1,    0,    0,    0,
+                                 0,    1,    0,    0,
+                                 0,    0,    1,    0,
+                                 0,    0,    0,    1 };
+static const Mat4f mI(mI_0);
+
+static const float m1_0[] = {    1,    3,    4,    0,
+                                 6,    7,    8,    5,
+                                98,    7,    6,    9,
+                                54,    3,    2,    5 };
+static const Mat4f m1(m1_0);
+
+static const float m1T_0[] = {   1,    6,   98,   54,
+                                 3,    7,    7,    3,
+                                 4,    8,    6,    2,
+                                 0,    5,    9,    5 };
+static const Mat4f m1T(m1T_0);
+
+static const float m2_0[] = {    1,    6,   98,   54,
+                                 3,    7,    7,    3,
+                                 4,    8,    6,    2,
+                                 0,    5,    9,    5 };
+static const Mat4f m2(m2_0);
+
+static const float m2xm1_0[] = {   26,   59,  143,   71,
+                                   59,  174,  730,  386,
+                                  143,  730, 9770, 5370,
+                                   71,  386, 5370, 2954 };
+static const Mat4f m2xm1(m2xm1_0);
+
+static const float m1xm2_0[] = {12557,  893,  748, 1182,
+                                  893,  116,  116,  113,
+                                  748,  116,  120,  104,
+                                 1182,  113,  104,  131 };
+static const Mat4f m1xm2(m1xm2_0);
+
+TEST_CASE( "Test 00 Load Get", "[mat4f][linear_algebra][math]" ) {
+    {
+        Mat4f m;
+        REQUIRE(mI == m);
+    }
+    {
+        float f16[16];
+        m1.get(f16);
+        COMPARE_NARRAYS_EPS(m1_0, f16, 16, EPSILON);
+
+        Mat4f m;
+        m.load(f16);
+        REQUIRE(m1 == m);
+    }
+}
+
+TEST_CASE( "Test 01 Mul", "[mat4f][linear_algebra][math]" ) {
+    {
+        REQUIRE(m1xm2 == m1 * m2);
+        Mat4f m; m.mul(m1, m2);
+        REQUIRE(m1xm2 == m);
+    }
+    {
+        REQUIRE(m2xm1 == m2 * m1);
+        Mat4f m; m.mul(m2, m1);
+        REQUIRE(m2xm1 == m);
+    }
+}
+
+TEST_CASE( "Test 02 Transpose", "[mat4f][linear_algebra][math]" ) {
+    REQUIRE(m1T == Mat4f(m1).transpose());
+    REQUIRE(m1T == Mat4f().transpose(m1));
+}
+
+TEST_CASE( "Test 80 LookAtNegZ", "[mat4f][linear_algebra][math]" ) {
+    Mat4f tmp;
+    Mat4f m;
+    // Look towards -z
+    m.setToLookAt(
+            Vec3f(0, 0,  0),  // eye
+            Vec3f(0, 0, -1),  // center
+            Vec3f(0, 1,  0),  // up
+            tmp);
+
+    /**
+     * The 3 rows of the matrix (= the 3 columns of the array/buffer) should be: side, up, -forward.
+     */
+    Mat4f exp( { 1, 0, 0, 0,
+                 0, 1, 0, 0,
+                 0, 0, 1, 0,
+                 0, 0, 0, 1  } );
+
+    REQUIRE(exp == m);
+}
+
+TEST_CASE( "Test 81 LookAtPosY", "[mat4f][linear_algebra][math]" ) {
+    Mat4f tmp;
+    Mat4f m;
+    // Look towards -z
+    m.setToLookAt(
+            Vec3f(0, 0, 0),  // eye
+            Vec3f(0, 1, 0),  // center
+            Vec3f(0, 0, 1),  // up
+            tmp);
+
+    /**
+     * The 3 rows of the matrix (= the 3 columns of the array/buffer) should be: side, up, -forward.
+     */
+    Mat4f exp( { 1, 0,  0, 0,
+                 0, 0, -1, 0,
+                 0, 1,  0, 0,
+                 0, 0,  0, 1
+            } );
+
+    REQUIRE(exp == m);
+}
+
diff --git a/test/test_math_mat4f_02_perf.cpp b/test/test_math_mat4f_02_perf.cpp
new file mode 100644
index 0000000..dfe6ab9
--- /dev/null
+++ b/test/test_math_mat4f_02_perf.cpp
@@ -0,0 +1,170 @@
+/*
+ * Author: Sven Gothel <[email protected]>
+ * Copyright (c) 2024 Gothel Software e.K.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be
+ * included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+#include <thread>
+#include <cinttypes>
+#include <cstring>
+
+#include <jau/test/catch2_ext.hpp>
+
+#include <jau/int_math.hpp>
+#include <jau/float_math.hpp>
+#include <jau/math/vec2f.hpp>
+#include <jau/math/vec2i.hpp>
+#include <jau/math/vec3f.hpp>
+#include <jau/math/vec4f.hpp>
+#include <jau/math/mat4f.hpp>
+#include <jau/math/quaternion.hpp>
+#include <jau/math/aabbox2f.hpp>
+#include <jau/math/aabbox3f.hpp>
+#include <jau/math/mat4f.hpp>
+#include <jau/math/recti.hpp>
+#include <jau/math/math_error.hpp>
+
+using namespace jau;
+using namespace jau::math;
+
+static const float m1_0[] = {    1,    3,    4,    0,
+                                 6,    7,    8,    5,
+                                98,    7,    6,    9,
+                                54,    3,    2,    5 };
+static const Mat4f m1(m1_0);
+
+static const float m2_0[] = {    1,    6,   98,   54,
+                                 3,    7,    7,    3,
+                                 4,    8,    6,    2,
+                                 0,    5,    9,    5 };
+static const Mat4f m2(m2_0);
+
+TEST_CASE( "Test 05 Perf01", "[mat4f][linear_algebra][math]" ) {
+    Mat4f res_m;
+
+    const size_t warmups = 1000_u64;
+    const size_t loops = 300_u64*1000000_u64;
+    jau::fraction_i64 tI4a = fractions_i64::zero;
+    jau::fraction_i64 tI4b = fractions_i64::zero;
+
+    const uint64_t tI5Max = 1000; // 1s
+    size_t loops5a = 0;
+    jau::fraction_i64 tI5a = fractions_i64::zero;
+    size_t loops5b = 0;
+    jau::fraction_i64 tI5b = fractions_i64::zero;
+
+    // avoid optimizing out unused computation results by simply adding up determinat
+    double dr = 1;
+
+    //
+    // Mat4f
+    //
+
+    // warm-up
+    for(size_t i=0; i<warmups; i++) {
+        res_m = m1 * m2;
+        dr += res_m.determinant();
+        res_m = m2 * m1;
+        dr += res_m.determinant();
+    }
+
+    jau::fraction_timespec t_0 = jau::getMonotonicTime();
+    for(size_t i=0; i<loops; i++) {
+        res_m = m1 * m2;
+        dr += res_m.determinant();
+        res_m = m2 * m1;
+        dr += res_m.determinant();
+    }
+    tI4a = (getMonotonicTime() - t_0).to_fraction_i64();
+    REQUIRE( dr > 0 );
+
+    // warm-up
+    for(size_t i=0; i<warmups; i++) {
+        res_m.load(m1);
+        res_m.mul(m2);
+        dr += res_m.determinant();
+        res_m.load(m2);
+        res_m.mul(m1);
+        dr += res_m.determinant();
+    }
+
+    t_0 = jau::getMonotonicTime();
+    for(size_t i=0; i<loops; i++) {
+        res_m.load(m1);
+        res_m.mul(m2);
+        dr += res_m.determinant();
+        res_m.load(m2);
+        res_m.mul(m1);
+        dr += res_m.determinant();
+    }
+    tI4b = (getMonotonicTime() - t_0).to_fraction_i64();
+    REQUIRE( dr > 0 );
+
+    tI5a = fractions_i64::zero;
+    t_0 = jau::getMonotonicTime();
+    uint64_t t_5 = jau::getCurrentMilliseconds();
+    uint64_t td_5=0;
+    while( td_5 < tI5Max ) {
+        res_m = m1 * m2;
+        dr += res_m.determinant();
+        res_m = m2 * m1;
+        dr += res_m.determinant();
+        ++loops5a;
+        // if( 0 == loops5a % 1000000 ) {
+            td_5 = jau::getCurrentMilliseconds() - t_5;
+        // }
+    }
+    tI5a = (getMonotonicTime() - t_0).to_fraction_i64();
+    REQUIRE( dr > 0 );
+
+    tI5b = fractions_i64::zero;
+    t_0 = jau::getMonotonicTime();
+    t_5 = jau::getCurrentMilliseconds();
+    td_5=0;
+    while( td_5 < tI5Max ) {
+        res_m.load(m1);
+        res_m.mul(m2);
+        dr += res_m.determinant();
+        res_m.load(m2);
+        res_m.mul(m1);
+        dr += res_m.determinant();
+        ++loops5b;
+        // if( 0 == loops5b % 1000000 ) {
+            td_5 = jau::getCurrentMilliseconds() - t_5;
+        // }
+    }
+    tI5b = (getMonotonicTime() - t_0).to_fraction_i64();
+    REQUIRE( dr > 0 );
+
+    printf("Checkmark %f\n", dr);
+    printf("Summary loops %6zu: I4a %6s ms total (%s us), %f ns/mul, I4a / I4b %f%%\n", loops,
+            jau::to_decstring(tI4a.to_ms()).c_str(), jau::to_decstring(tI4a.to_us()).c_str(),
+            (double)tI4a.to_ns()/2.0/(double)loops, tI4a.to_double()/tI4b.to_double()*100.0);
+    printf("Summary loops %6zu: I4b %6s ms total (%s us), %f ns/mul, I4b / I4a %f%%\n", loops,
+            jau::to_decstring(tI4b.to_ms()).c_str(), jau::to_decstring(tI4b.to_us()).c_str(),
+            (double)tI4b.to_ns()/2.0/(double)loops, tI4b.to_double()/tI4a.to_double()*100.0);
+
+    printf("Summary loops %6zu: I5a %6s ms total, %f ns/mul, I5a / I5b %f%%\n", loops5a,
+            jau::to_decstring(tI5a.to_ms()).c_str(),
+            (double)tI5a.to_ns()/2.0/(double)loops5a, tI5a.to_double()/tI5b.to_double()*100.0);
+    printf("Summary loops %6zu: I5b %6s ms total, %f ns/mul, I5b / I5a %f%%\n", loops5b,
+            jau::to_decstring(tI5b.to_ms()).c_str(),
+            (double)tI5b.to_ns()/2.0/(double)loops5b, tI5b.to_double()/tI5a.to_double()*100.0);
+}