GPU: Add Math libraries to GPU shaders code

This implement most of the functions provided by the BLI math library.
This is part of the effort to unify GLSL and C++ syntax. Ref T103026.

This also adds some infrastructure to make it possible to run GLSL shader unit
test.

Some code already present in other libs is being copied to the new libs.
This patch does not make use of the new libs outside of the tests.

Note that the test is still crashing when using metal.

source/blender/blenlib/BLI_math_matrix_types.hh

@@ -30,6 +30,8 @@
  * defined outside of the class in the `blender::math` namespace.
  */
 
+#define __BLI_MATH_MATRIX_TYPES_HH__
+
 #include <array>
 #include <cmath>
 #include <iostream>

source/blender/draw/intern/shaders/common_math_lib.glsl

@@ -16,8 +16,11 @@
 #define M_1_PI2 0.101321183642337771443  /* 1/(pi^2) */
 #define M_SQRT2 1.41421356237309504880   /* sqrt(2) */
 #define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */
-#define FLT_MAX 3.402823e+38
-#define FLT_MIN 1.175494e-38
+#ifndef FLT_MAX
+#  define FLT_MAX 3.402823e+38
+#  define FLT_MIN 1.175494e-38
+#  define FLT_EPSILON 1.192092896e-07F
+#endif
 
 vec3 mul(mat3 m, vec3 v)
 {

source/blender/gpu/CMakeLists.txt

@@ -388,6 +388,13 @@ set(GLSL_SRC
   shaders/common/gpu_shader_common_math.glsl
   shaders/common/gpu_shader_common_math_utils.glsl
   shaders/common/gpu_shader_common_mix_rgb.glsl
+  shaders/common/gpu_shader_math_base_lib.glsl
+  shaders/common/gpu_shader_math_fast_lib.glsl
+  shaders/common/gpu_shader_math_matrix_lib.glsl
+  shaders/common/gpu_shader_math_rotation_lib.glsl
+  shaders/common/gpu_shader_math_vector_lib.glsl
+  shaders/common/gpu_shader_test_lib.glsl
+  shaders/common/gpu_shader_utildefines_lib.glsl
 
   shaders/material/gpu_shader_material_add_shader.glsl
   shaders/material/gpu_shader_material_ambient_occlusion.glsl
@@ -484,6 +491,8 @@ set(GLSL_SRC
   shaders/gpu_shader_cfg_world_clip_lib.glsl
   shaders/gpu_shader_colorspace_lib.glsl
 
+  tests/gpu_math_test.glsl
+
   GPU_shader_shared_utils.h
 )
 
@@ -625,6 +634,7 @@ set(SRC_SHADER_CREATE_INFOS
   shaders/infos/gpu_shader_keyframe_shape_info.hh
   shaders/infos/gpu_shader_line_dashed_uniform_color_info.hh
   shaders/infos/gpu_shader_simple_lighting_info.hh
+  shaders/infos/gpu_shader_test_info.hh
   shaders/infos/gpu_shader_text_info.hh
   shaders/infos/gpu_srgb_to_framebuffer_space_info.hh
 )

source/blender/gpu/GPU_shader_shared.h

@@ -7,6 +7,8 @@
 
 #ifndef USE_GPU_SHADER_CREATE_INFO
 #  include "GPU_shader_shared_utils.h"
+
+typedef struct TestOutputRawData TestOutputRawData;
 #endif
 
 struct NodeLinkData {
@@ -66,3 +68,66 @@ struct MultiRectCallData {
   float4 calls_data[MAX_CALLS * 3];
 };
 BLI_STATIC_ASSERT_ALIGN(struct MultiRectCallData, 16)
+
+enum TestStatus {
+  TEST_STATUS_NONE = 0,
+  TEST_STATUS_PASSED = 1,
+  TEST_STATUS_FAILED = 2,
+};
+enum TestType {
+  TEST_TYPE_BOOL = 0,
+  TEST_TYPE_UINT = 1,
+  TEST_TYPE_INT = 2,
+  TEST_TYPE_FLOAT = 3,
+  TEST_TYPE_IVEC2 = 4,
+  TEST_TYPE_IVEC3 = 5,
+  TEST_TYPE_IVEC4 = 6,
+  TEST_TYPE_UVEC2 = 7,
+  TEST_TYPE_UVEC3 = 8,
+  TEST_TYPE_UVEC4 = 9,
+  TEST_TYPE_VEC2 = 10,
+  TEST_TYPE_VEC3 = 11,
+  TEST_TYPE_VEC4 = 12,
+  TEST_TYPE_MAT2X2 = 13,
+  TEST_TYPE_MAT2X3 = 14,
+  TEST_TYPE_MAT2X4 = 15,
+  TEST_TYPE_MAT3X2 = 16,
+  TEST_TYPE_MAT3X3 = 17,
+  TEST_TYPE_MAT3X4 = 18,
+  TEST_TYPE_MAT4X2 = 19,
+  TEST_TYPE_MAT4X3 = 20,
+  TEST_TYPE_MAT4X4 = 21,
+};
+
+/** \note Contains arrays of scalar. To be use only with SSBOs to avoid padding issues. */
+struct TestOutputRawData {
+  uint data[16];
+};
+BLI_STATIC_ASSERT_ALIGN(struct TestOutputRawData, 16)
+
+struct TestOutput {
+  TestOutputRawData expect;
+  TestOutputRawData result;
+  /** TestStatus. */
+  uint status;
+  /** Line error in the glsl file. */
+  int line;
+  /** TestType of expect and result. */
+  uint type;
+  int _pad0;
+};
+BLI_STATIC_ASSERT_ALIGN(struct TestOutput, 16)
+
+#ifdef GPU_SHADER
+TestOutput test_output(
+    TestOutputRawData expect, TestOutputRawData result, bool status, int line, uint type)
+{
+  TestOutput test;
+  test.expect = expect;
+  test.result = result;
+  test.status = status ? TEST_STATUS_PASSED : TEST_STATUS_FAILED;
+  test.line = line;
+  test.type = type;
+  return test;
+}
+#endif

source/blender/gpu/GPU_shader_shared_utils.h

@@ -67,7 +67,9 @@
 #  include "BLI_assert.h"
 
 #  ifdef __cplusplus
-#    include "BLI_float4x4.hh"
+#    ifndef __BLI_MATH_MATRIX_TYPES_HH__
+#      include "BLI_float4x4.hh"
+#    endif
 #    include "BLI_math_vector_types.hh"
 using blender::float2;
 using blender::float3;

source/blender/gpu/shaders/common/gpu_shader_math_base_lib.glsl

@@ -0,0 +1,152 @@
+
+/* WORKAROUND: to guard against double include in EEVEE. */
+#ifndef GPU_SHADER_MATH_BASE_LIB_GLSL
+#define GPU_SHADER_MATH_BASE_LIB_GLSL
+
+#define M_PI 3.14159265358979323846      /* pi */
+#define M_TAU 6.28318530717958647692     /* tau = 2*pi */
+#define M_PI_2 1.57079632679489661923    /* pi/2 */
+#define M_PI_4 0.78539816339744830962    /* pi/4 */
+#define M_SQRT2 1.41421356237309504880   /* sqrt(2) */
+#define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */
+#define M_SQRT3 1.73205080756887729352   /* sqrt(3) */
+#define M_SQRT1_3 0.57735026918962576450 /* 1/sqrt(3) */
+#define M_1_PI 0.318309886183790671538   /* 1/pi */
+#define M_E 2.7182818284590452354        /* e */
+#define M_LOG2E 1.4426950408889634074    /* log_2 e */
+#define M_LOG10E 0.43429448190325182765  /* log_10 e */
+#define M_LN2 0.69314718055994530942     /* log_e 2 */
+#define M_LN10 2.30258509299404568402    /* log_e 10 */
+
+/* `powf` is really slow for raising to integer powers. */
+
+float pow2f(float x)
+{
+  return x * x;
+}
+float pow3f(float x)
+{
+  return x * x * x;
+}
+float pow4f(float x)
+{
+  return pow2f(pow2f(x));
+}
+float pow5f(float x)
+{
+  return pow4f(x) * x;
+}
+float pow6f(float x)
+{
+  return pow2f(pow3f(x));
+}
+float pow7f(float x)
+{
+  return pow6f(x) * x;
+}
+float pow8f(float x)
+{
+  return pow2f(pow4f(x));
+}
+
+int square_i(int v)
+{
+  return v * v;
+}
+uint square_uint(uint v)
+{
+  return v * v;
+}
+float square_f(float v)
+{
+  return v * v;
+}
+
+int cube_i(int v)
+{
+  return v * v * v;
+}
+uint cube_uint(uint v)
+{
+  return v * v * v;
+}
+float cube_f(float v)
+{
+  return v * v * v;
+}
+
+float hypot(float x, float y)
+{
+  return sqrt(x * x + y * y);
+}
+
+float atan2(float y, float x)
+{
+  return atan(y, x);
+}
+
+/**
+ * Safe `a` modulo `b`.
+ * If `b` equal 0 the result will be 0.
+ */
+float safe_mod(float a, float b)
+{
+  return (b != 0.0) ? mod(a, b) : 0.0;
+}
+
+/**
+ * Returns \a a if it is a multiple of \a b or the next multiple or \a b after \b a .
+ * In other words, it is equivalent to `divide_ceil(a, b) * b`.
+ * It is undefined if \a a is negative or \b b is not strictly positive.
+ */
+int ceil_to_multiple(int a, int b)
+{
+  return ((a + b - 1) / b) * b;
+}
+uint ceil_to_multiple(uint a, uint b)
+{
+  return ((a + b - 1u) / b) * b;
+}
+
+/**
+ * Integer division that returns the ceiling, instead of flooring like normal C division.
+ * It is undefined if \a a is negative or \b b is not strictly positive.
+ */
+int divide_ceil(int a, int b)
+{
+  return (a + b - 1) / b;
+}
+uint divide_ceil(uint a, uint b)
+{
+  return (a + b - 1u) / b;
+}
+
+/**
+ * Component wise, use vector to replace min if it is smaller and max if bigger.
+ */
+void min_max(float value, inout float min_v, inout float max_v)
+{
+  min_v = min(value, min_v);
+  max_v = max(value, max_v);
+}
+
+/**
+ * Safe divide `a` by `b`.
+ * If `b` equal 0 the result will be 0.
+ */
+float safe_divide(float a, float b)
+{
+  return (b != 0.0) ? (a / b) : 0.0;
+}
+
+/**
+ * Return true if the difference between`a` and `b` is below the `epsilon` value.
+ */
+bool is_equal(float a, float b, const float epsilon)
+{
+  return abs(a - b) <= epsilon;
+}
+
+/** \} */
+
+#endif /* GPU_SHADER_MATH_BASE_LIB_GLSL */

source/blender/gpu/shaders/common/gpu_shader_math_fast_lib.glsl

@@ -0,0 +1,32 @@
+
+/* WORKAROUND: to guard against double include in EEVEE. */
+#ifndef GPU_SHADER_MATH_FAST_LIB_GLSL
+#define GPU_SHADER_MATH_FAST_LIB_GLSL
+
+/* [Drobot2014a] Low Level Optimizations for GCN. */
+float sqrt_fast(float v)
+{
+  return intBitsToFloat(0x1fbd1df5 + (floatBitsToInt(v) >> 1));
+}
+vec2 sqrt_fast(vec2 v)
+{
+  return intBitsToFloat(0x1fbd1df5 + (floatBitsToInt(v) >> 1));
+}
+
+/* [Eberly2014] GPGPU Programming for Games and Science. */
+float acos_fast(float v)
+{
+  float res = -0.156583 * abs(v) + M_PI_2;
+  res *= sqrt_fast(1.0 - abs(v));
+  return (v >= 0) ? res : M_PI - res;
+}
+vec2 acos_fast(vec2 v)
+{
+  vec2 res = -0.156583 * abs(v) + M_PI_2;
+  res *= sqrt_fast(1.0 - abs(v));
+  v.x = (v.x >= 0) ? res.x : M_PI - res.x;
+  v.y = (v.y >= 0) ? res.y : M_PI - res.y;
+  return v;
+}
+
+#endif /* GPU_SHADER_MATH_FAST_LIB_GLSL */

source/blender/gpu/shaders/common/gpu_shader_math_rotation_lib.glsl

@@ -0,0 +1,186 @@
+
+#pragma BLENDER_REQUIRE(gpu_shader_utildefines_lib.glsl)
+
+/* WORKAROUND: to guard against double include in EEVEE. */
+#ifndef GPU_SHADER_MATH_ROTATION_LIB_GLSL
+#  define GPU_SHADER_MATH_ROTATION_LIB_GLSL
+
+/* -------------------------------------------------------------------- */
+/** \name Rotation Types
+ * \{ */
+
+struct Angle {
+  /* Angle in radian. */
+  float angle;
+
+#  ifdef GPU_METAL
+  Angle() = default;
+  Angle(float angle_) : angle(angle_){};
+#  endif
+};
+
+struct AxisAngle {
+  vec3 axis;
+  float angle;
+
+#  ifdef GPU_METAL
+  AxisAngle() = default;
+  AxisAngle(vec3 axis_, float angle_) : axis(axis_), angle(angle_){};
+#  endif
+};
+
+AxisAngle AxisAngle__identity()
+{
+  return AxisAngle(vec3(0, 1, 0), 0);
+}
+
+struct Quaternion {
+  float x, y, z, w;
+#  ifdef GPU_METAL
+  Quaternion() = default;
+  Quaternion(float x_, float y_, float z_, float w_) : x(x_), y(y_), z(z_), w(w_){};
+#  endif
+};
+
+vec4 as_vec4(Quaternion quat)
+{
+  return vec4(quat.x, quat.y, quat.z, quat.w);
+}
+
+Quaternion Quaternion__identity()
+{
+  return Quaternion(1, 0, 0, 0);
+}
+
+struct EulerXYZ {
+  float x, y, z;
+#  ifdef GPU_METAL
+  EulerXYZ() = default;
+  EulerXYZ(float x_, float y_, float z_) : x(x_), y(y_), z(z_){};
+#  endif
+};
+
+vec3 as_vec3(EulerXYZ eul)
+{
+  return vec3(eul.x, eul.y, eul.z);
+}
+
+EulerXYZ EulerXYZ__identity()
+{
+  return EulerXYZ(0, 0, 0);
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+/** \name Rotation Functions
+ * \{ */
+
+/**
+ * Generic function for implementing slerp
+ * (quaternions and spherical vector coords).
+ *
+ * \param t: factor in [0..1]
+ * \param cosom: dot product from normalized vectors/quats.
+ * \param r_w: calculated weights.
+ */
+vec2 interpolate_dot_slerp(float t, float cosom)
+{
+  vec2 w = vec2(1.0 - t, t);
+  /* Within [-1..1] range, avoid aligned axis. */
+  const float eps = 1e-4;
+  if (abs(cosom) < 1.0 - eps) {
+    float omega = acos(cosom);
+    w = sin(w * omega) / sin(omega);
+  }
+  return w;
+}
+
+Quaternion interpolate(Quaternion a, Quaternion b, float t)
+{
+  vec4 quat = as_vec4(a);
+  float cosom = dot(as_vec4(a), as_vec4(b));
+  /* Rotate around shortest angle. */
+  if (cosom < 0.0) {
+    cosom = -cosom;
+    quat = -quat;
+  }
+  vec2 w = interpolate_dot_slerp(t, cosom);
+  quat = w.x * quat + w.y * as_vec4(b);
+  return Quaternion(UNPACK4(quat));
+}
+
+Quaternion to_quaternion(EulerXYZ eul)
+{
+  float ti = eul.x * 0.5;
+  float tj = eul.y * 0.5;
+  float th = eul.z * 0.5;
+  float ci = cos(ti);
+  float cj = cos(tj);
+  float ch = cos(th);
+  float si = sin(ti);
+  float sj = sin(tj);
+  float sh = sin(th);
+  float cc = ci * ch;
+  float cs = ci * sh;
+  float sc = si * ch;
+  float ss = si * sh;
+
+  Quaternion quat;
+  quat.x = cj * cc + sj * ss;
+  quat.y = cj * sc - sj * cs;
+  quat.z = cj * ss + sj * cc;
+  quat.w = cj * cs - sj * sc;
+  return quat;
+}
+
+Quaternion to_axis_angle(AxisAngle axis_angle)
+{
+  float angle_cos = cos(axis_angle.angle);
+  /** Using half angle identities: sin(angle / 2) = sqrt((1 - angle_cos) / 2) */
+  float sine = sqrt(0.5 - angle_cos * 0.5);
+  float cosine = sqrt(0.5 + angle_cos * 0.5);
+
+  /* TODO(fclem): Optimize. */
+  float angle_sin = sin(axis_angle.angle);
+  if (angle_sin < 0.0) {
+    sine = -sine;
+  }
+
+  Quaternion quat;
+  quat.x = cosine;
+  quat.y = axis_angle.axis.x * sine;
+  quat.z = axis_angle.axis.y * sine;
+  quat.w = axis_angle.axis.z * sine;
+  return quat;
+}
+
+AxisAngle to_axis_angle(Quaternion quat)
+{
+  /* Calculate angle/2, and sin(angle/2). */
+  float ha = acos(quat.x);
+  float si = sin(ha);
+
+  /* From half-angle to angle. */
+  float angle = ha * 2;
+  /* Prevent division by zero for axis conversion. */
+  if (abs(si) < 0.0005) {
+    si = 1.0;
+  }
+
+  vec3 axis = vec3(quat.y, quat.z, quat.w) / si;
+  if (is_zero(axis)) {
+    axis[1] = 1.0;
+  }
+  return AxisAngle(axis, angle);
+}
+
+AxisAngle to_axis_angle(EulerXYZ eul)
+{
+  /* Use quaternions as intermediate representation for now... */
+  return to_axis_angle(to_quaternion(eul));
+}
+
+/** \} */
+
+#endif /* GPU_SHADER_MATH_ROTATION_LIB_GLSL */

source/blender/gpu/shaders/common/gpu_shader_math_vector_lib.glsl

@@ -0,0 +1,535 @@
+
+#pragma BLENDER_REQUIRE(gpu_shader_math_base_lib.glsl)
+
+/* WORKAROUND: to guard against double include in EEVEE. */
+#ifndef GPU_SHADER_MATH_VECTOR_LIB_GLSL
+#  define GPU_SHADER_MATH_VECTOR_LIB_GLSL
+
+/* Metal does not need prototypes. */
+#  ifndef GPU_METAL
+
+/**
+ * Return true if all components is equal to zero.
+ */
+bool is_zero(vec2 vec);
+bool is_zero(vec3 vec);
+bool is_zero(vec4 vec);
+
+/**
+ * Return true if any component is equal to zero.
+ */
+bool is_any_zero(vec2 vec);
+bool is_any_zero(vec3 vec);
+bool is_any_zero(vec4 vec);
+
+/**
+ * Return true if the deference between`a` and `b` is below the `epsilon` value.
+ * Epsilon value is scaled by magnitude of `a` before comparison.
+ */
+bool almost_equal_relative(vec2 a, vec2 b, const float epsilon_factor);
+bool almost_equal_relative(vec3 a, vec3 b, const float epsilon_factor);
+bool almost_equal_relative(vec4 a, vec4 b, const float epsilon_factor);
+
+/**
+ * Safe `a` modulo `b`.
+ * If `b` equal 0 the result will be 0.
+ */
+vec2 safe_mod(vec2 a, vec2 b);
+vec3 safe_mod(vec3 a, vec3 b);
+vec4 safe_mod(vec4 a, vec4 b);
+vec2 safe_mod(vec2 a, float b);
+vec3 safe_mod(vec3 a, float b);
+vec4 safe_mod(vec4 a, float b);
+
+/**
+ * Returns \a a if it is a multiple of \a b or the next multiple or \a b after \b a .
+ * In other words, it is equivalent to `divide_ceil(a, b) * b`.
+ * It is undefined if \a a is negative or \b b is not strictly positive.
+ */
+ivec2 ceil_to_multiple(ivec2 a, ivec2 b);
+ivec3 ceil_to_multiple(ivec3 a, ivec3 b);
+ivec4 ceil_to_multiple(ivec4 a, ivec4 b);
+uvec2 ceil_to_multiple(uvec2 a, uvec2 b);
+uvec3 ceil_to_multiple(uvec3 a, uvec3 b);
+uvec4 ceil_to_multiple(uvec4 a, uvec4 b);
+
+/**
+ * Integer division that returns the ceiling, instead of flooring like normal C division.
+ * It is undefined if \a a is negative or \b b is not strictly positive.
+ */
+ivec2 divide_ceil(ivec2 a, ivec2 b);
+ivec3 divide_ceil(ivec3 a, ivec3 b);
+ivec4 divide_ceil(ivec4 a, ivec4 b);
+uvec2 divide_ceil(uvec2 a, uvec2 b);
+uvec3 divide_ceil(uvec3 a, uvec3 b);
+uvec4 divide_ceil(uvec4 a, uvec4 b);
+
+/**
+ * Component wise, use vector to replace min if it is smaller and max if bigger.
+ */
+void min_max(vec2 vector, inout vec2 min, inout vec2 max);
+void min_max(vec3 vector, inout vec3 min, inout vec3 max);
+void min_max(vec4 vector, inout vec4 min, inout vec4 max);
+
+/**
+ * Safe divide `a` by `b`.
+ * If `b` equal 0 the result will be 0.
+ */
+vec2 safe_divide(vec2 a, vec2 b);
+vec3 safe_divide(vec3 a, vec3 b);
+vec4 safe_divide(vec4 a, vec4 b);
+vec2 safe_divide(vec2 a, float b);
+vec3 safe_divide(vec3 a, float b);
+vec4 safe_divide(vec4 a, float b);
+
+/**
+ * Return the manhattan length of `a`.
+ * This is also the sum of the absolute value of all components.
+ */
+float length_manhattan(vec2 a);
+float length_manhattan(vec3 a);
+float length_manhattan(vec4 a);
+
+/**
+ * Return the length squared of `a`.
+ */
+float length_squared(vec2 a);
+float length_squared(vec3 a);
+float length_squared(vec4 a);
+
+/**
+ * Return the manhattan distance between `a` and `b`.
+ */
+float distance_manhattan(vec2 a, vec2 b);
+float distance_manhattan(vec3 a, vec3 b);
+float distance_manhattan(vec4 a, vec4 b);
+
+/**
+ * Return the squared distance between `a` and `b`.
+ */
+float distance_squared(vec2 a, vec2 b);
+float distance_squared(vec3 a, vec3 b);
+float distance_squared(vec4 a, vec4 b);
+
+/**
+ * Return the projection of `p` onto `v_proj`.
+ */
+vec3 project(vec3 p, vec3 v_proj);
+
+/**
+ * Return normalized version of the `vector` and its length.
+ */
+vec2 normalize_and_get_length(vec2 vector, out float out_length);
+vec3 normalize_and_get_length(vec3 vector, out float out_length);
+vec4 normalize_and_get_length(vec4 vector, out float out_length);
+
+/**
+ * Per component linear interpolation.
+ */
+vec2 interpolate(vec2 a, vec2 b, float t);
+vec3 interpolate(vec3 a, vec3 b, float t);
+vec4 interpolate(vec4 a, vec4 b, float t);
+
+/**
+ * Return half-way point between `a` and  `b`.
+ */
+vec2 midpoint(vec2 a, vec2 b);
+vec3 midpoint(vec3 a, vec3 b);
+vec4 midpoint(vec4 a, vec4 b);
+
+/**
+ * Return `vector` if `incident` and `reference` are pointing in the same direction.
+ */
+// vec2 faceforward(vec2 vector, vec2 incident, vec2 reference); /* Built-in GLSL. */
+
+/**
+ * Return the index of the component with the greatest absolute value.
+ */
+int dominant_axis(vec3 a);
+
+/**
+ * Calculates a perpendicular vector to \a v.
+ * \note Returned vector can be in any perpendicular direction.
+ * \note Returned vector might not the same length as \a v.
+ */
+vec3 orthogonal(vec3 v);
+/**
+ * Calculates a perpendicular vector to \a v.
+ * \note Returned vector is always rotated 90° counter clock wise.
+ */
+vec2 orthogonal(vec2 v);
+
+/**
+ * Return true if the difference between`a` and `b` is below the `epsilon` value.
+ */
+bool is_equal(vec2 a, vec2 b, const float epsilon);
+bool is_equal(vec3 a, vec3 b, const float epsilon);
+bool is_equal(vec4 a, vec4 b, const float epsilon);
+
+#  endif /* GPU_METAL */
+
+/* ---------------------------------------------------------------------- */
+/** \name Implementation
+ * \{ */
+
+#  ifdef GPU_METAL /* Already defined in shader_defines.msl/glsl to move here. */
+bool is_zero(vec2 vec)
+{
+  return all(equal(vec, vec2(0.0)));
+}
+bool is_zero(vec3 vec)
+{
+  return all(equal(vec, vec3(0.0)));
+}
+bool is_zero(vec4 vec)
+{
+  return all(equal(vec, vec4(0.0)));
+}
+#  endif
+
+bool is_any_zero(vec2 vec)
+{
+  return any(equal(vec, vec2(0.0)));
+}
+bool is_any_zero(vec3 vec)
+{
+  return any(equal(vec, vec3(0.0)));
+}
+bool is_any_zero(vec4 vec)
+{
+  return any(equal(vec, vec4(0.0)));
+}
+
+bool almost_equal_relative(vec2 a, vec2 b, const float epsilon_factor)
+{
+  for (int i = 0; i < 2; i++) {
+    if (abs(a[i] - b[i]) > epsilon_factor * abs(a[i])) {
+      return false;
+    }
+  }
+  return true;
+}
+bool almost_equal_relative(vec3 a, vec3 b, const float epsilon_factor)
+{
+  for (int i = 0; i < 3; i++) {
+    if (abs(a[i] - b[i]) > epsilon_factor * abs(a[i])) {
+      return false;
+    }
+  }
+  return true;
+}
+bool almost_equal_relative(vec4 a, vec4 b, const float epsilon_factor)
+{
+  for (int i = 0; i < 4; i++) {
+    if (abs(a[i] - b[i]) > epsilon_factor * abs(a[i])) {
+      return false;
+    }
+  }
+  return true;
+}
+
+vec2 safe_mod(vec2 a, vec2 b)
+{
+  return select(vec2(0), mod(a, b), notEqual(b, vec2(0)));
+}
+vec3 safe_mod(vec3 a, vec3 b)
+{
+  return select(vec3(0), mod(a, b), notEqual(b, vec3(0)));
+}
+vec4 safe_mod(vec4 a, vec4 b)
+{
+  return select(vec4(0), mod(a, b), notEqual(b, vec4(0)));
+}
+
+vec2 safe_mod(vec2 a, float b)
+{
+  return (b != 0.0) ? mod(a, vec2(b)) : vec2(0);
+}
+vec3 safe_mod(vec3 a, float b)
+{
+  return (b != 0.0) ? mod(a, vec3(b)) : vec3(0);
+}
+vec4 safe_mod(vec4 a, float b)
+{
+  return (b != 0.0) ? mod(a, vec4(b)) : vec4(0);
+}
+
+ivec2 ceil_to_multiple(ivec2 a, ivec2 b)
+{
+  return ((a + b - 1) / b) * b;
+}
+ivec3 ceil_to_multiple(ivec3 a, ivec3 b)
+{
+  return ((a + b - 1) / b) * b;
+}
+ivec4 ceil_to_multiple(ivec4 a, ivec4 b)
+{
+  return ((a + b - 1) / b) * b;
+}
+uvec2 ceil_to_multiple(uvec2 a, uvec2 b)
+{
+  return ((a + b - 1u) / b) * b;
+}
+uvec3 ceil_to_multiple(uvec3 a, uvec3 b)
+{
+  return ((a + b - 1u) / b) * b;
+}
+uvec4 ceil_to_multiple(uvec4 a, uvec4 b)
+{
+  return ((a + b - 1u) / b) * b;
+}
+
+ivec2 divide_ceil(ivec2 a, ivec2 b)
+{
+  return (a + b - 1) / b;
+}
+ivec3 divide_ceil(ivec3 a, ivec3 b)
+{
+  return (a + b - 1) / b;
+}
+ivec4 divide_ceil(ivec4 a, ivec4 b)
+{
+  return (a + b - 1) / b;
+}
+uvec2 divide_ceil(uvec2 a, uvec2 b)
+{
+  return (a + b - 1u) / b;
+}
+uvec3 divide_ceil(uvec3 a, uvec3 b)
+{
+  return (a + b - 1u) / b;
+}
+uvec4 divide_ceil(uvec4 a, uvec4 b)
+{
+  return (a + b - 1u) / b;
+}
+
+void min_max(vec2 vector, inout vec2 min_v, inout vec2 max_v)
+{
+  min_v = min(vector, min_v);
+  max_v = max(vector, max_v);
+}
+void min_max(vec3 vector, inout vec3 min_v, inout vec3 max_v)
+{
+  min_v = min(vector, min_v);
+  max_v = max(vector, max_v);
+}
+void min_max(vec4 vector, inout vec4 min_v, inout vec4 max_v)
+{
+  min_v = min(vector, min_v);
+  max_v = max(vector, max_v);
+}
+
+vec2 safe_divide(vec2 a, vec2 b)
+{
+  return select(vec2(0), a / b, notEqual(b, vec2(0)));
+}
+vec3 safe_divide(vec3 a, vec3 b)
+{
+  return select(vec3(0), a / b, notEqual(b, vec3(0)));
+}
+vec4 safe_divide(vec4 a, vec4 b)
+{
+  return select(vec4(0), a / b, notEqual(b, vec4(0)));
+}
+
+vec2 safe_divide(vec2 a, float b)
+{
+  return (b != 0.0) ? (a / b) : vec2(0);
+}
+vec3 safe_divide(vec3 a, float b)
+{
+  return (b != 0.0) ? (a / b) : vec3(0);
+}
+vec4 safe_divide(vec4 a, float b)
+{
+  return (b != 0.0) ? (a / b) : vec4(0);
+}
+
+float length_manhattan(vec2 a)
+{
+  return dot(abs(a), vec2(1));
+}
+float length_manhattan(vec3 a)
+{
+  return dot(abs(a), vec3(1));
+}
+float length_manhattan(vec4 a)
+{
+  return dot(abs(a), vec4(1));
+}
+
+float length_squared(vec2 a)
+{
+  return dot(a, a);
+}
+float length_squared(vec3 a)
+{
+  return dot(a, a);
+}
+float length_squared(vec4 a)
+{
+  return dot(a, a);
+}
+
+float distance_manhattan(vec2 a, vec2 b)
+{
+  return length_manhattan(a - b);
+}
+float distance_manhattan(vec3 a, vec3 b)
+{
+  return length_manhattan(a - b);
+}
+float distance_manhattan(vec4 a, vec4 b)
+{
+  return length_manhattan(a - b);
+}
+
+float distance_squared(vec2 a, vec2 b)
+{
+  return length_squared(a - b);
+}
+float distance_squared(vec3 a, vec3 b)
+{
+  return length_squared(a - b);
+}
+float distance_squared(vec4 a, vec4 b)
+{
+  return length_squared(a - b);
+}
+
+vec3 project(vec3 p, vec3 v_proj)
+{
+  if (is_zero(v_proj)) {
+    return vec3(0);
+  }
+  return v_proj * (dot(p, v_proj) / dot(v_proj, v_proj));
+}
+
+vec2 normalize_and_get_length(vec2 vector, out float out_length)
+{
+  out_length = length_squared(vector);
+  const float threshold = 1e-35f;
+  if (out_length > threshold) {
+    out_length = sqrt(out_length);
+    return vector / out_length;
+  }
+  /* Either the vector is small or one of it's values contained `nan`. */
+  out_length = 0.0;
+  return vec2(0.0);
+}
+vec3 normalize_and_get_length(vec3 vector, out float out_length)
+{
+  out_length = length_squared(vector);
+  const float threshold = 1e-35f;
+  if (out_length > threshold) {
+    out_length = sqrt(out_length);
+    return vector / out_length;
+  }
+  /* Either the vector is small or one of it's values contained `nan`. */
+  out_length = 0.0;
+  return vec3(0.0);
+}
+vec4 normalize_and_get_length(vec4 vector, out float out_length)
+{
+  out_length = length_squared(vector);
+  const float threshold = 1e-35f;
+  if (out_length > threshold) {
+    out_length = sqrt(out_length);
+    return vector / out_length;
+  }
+  /* Either the vector is small or one of it's values contained `nan`. */
+  out_length = 0.0;
+  return vec4(0.0);
+}
+
+vec2 interpolate(vec2 a, vec2 b, float t)
+{
+  return mix(a, b, t);
+}
+vec3 interpolate(vec3 a, vec3 b, float t)
+{
+  return mix(a, b, t);
+}
+vec4 interpolate(vec4 a, vec4 b, float t)
+{
+  return mix(a, b, t);
+}
+
+vec2 midpoint(vec2 a, vec2 b)
+{
+  return (a + b) * 0.5;
+}
+vec3 midpoint(vec3 a, vec3 b)
+{
+  return (a + b) * 0.5;
+}
+vec4 midpoint(vec4 a, vec4 b)
+{
+  return (a + b) * 0.5;
+}
+
+int dominant_axis(vec3 a)
+{
+  vec3 b = abs(a);
+  return ((b.x > b.y) ? ((b.x > b.z) ? 0 : 2) : ((b.y > b.z) ? 1 : 2));
+}
+
+vec3 orthogonal(vec3 v)
+{
+  switch (dominant_axis(v)) {
+    default:
+    case 0:
+      return vec3(-v.y - v.z, v.x, v.x);
+    case 1:
+      return vec3(v.y, -v.x - v.z, v.y);
+    case 2:
+      return vec3(v.z, v.z, -v.x - v.y);
+  }
+}
+
+vec2 orthogonal(vec2 v)
+{
+  return vec2(-v.y, v.x);
+}
+
+bool is_equal(vec2 a, vec2 b, const float epsilon)
+{
+  return all(lessThanEqual(abs(a - b), vec2(epsilon)));
+}
+bool is_equal(vec3 a, vec3 b, const float epsilon)
+{
+  return all(lessThanEqual(abs(a - b), vec3(epsilon)));
+}
+bool is_equal(vec4 a, vec4 b, const float epsilon)
+{
+  return all(lessThanEqual(abs(a - b), vec4(epsilon)));
+}
+
+#  define ASSERT_UNIT_EPSILON 0.0002
+
+/* Checks are flipped so NAN doesn't assert because we're making sure the value was
+ * normalized and in the case we don't want NAN to be raising asserts since there
+ * is nothing to be done in that case. */
+bool is_unit_scale(vec2 v)
+{
+  float test_unit = length_squared(v);
+  return (!(abs(test_unit - 1.0) >= ASSERT_UNIT_EPSILON) ||
+          !(abs(test_unit) >= ASSERT_UNIT_EPSILON));
+}
+bool is_unit_scale(vec3 v)
+{
+  float test_unit = length_squared(v);
+  return (!(abs(test_unit - 1.0) >= ASSERT_UNIT_EPSILON) ||
+          !(abs(test_unit) >= ASSERT_UNIT_EPSILON));
+}
+bool is_unit_scale(vec4 v)
+{
+  float test_unit = length_squared(v);
+  return (!(abs(test_unit - 1.0) >= ASSERT_UNIT_EPSILON) ||
+          !(abs(test_unit) >= ASSERT_UNIT_EPSILON));
+}
+
+/** \} */
+
+#endif /* GPU_SHADER_MATH_VECTOR_LIB_GLSL */

source/blender/gpu/shaders/common/gpu_shader_test_lib.glsl

@@ -0,0 +1,214 @@
+
+// clang-format off
+#ifndef GPU_METAL
+bool is_integer(bool v) { return true; }
+#endif
+bool is_integer(uint v) { return true; }
+bool is_integer(int v) { return true; }
+bool is_integer(float v) { return false; }
+bool is_integer(ivec2 v) { return true; }
+bool is_integer(ivec3 v) { return true; }
+bool is_integer(ivec4 v) { return true; }
+bool is_integer(uvec2 v) { return true; }
+bool is_integer(uvec3 v) { return true; }
+bool is_integer(uvec4 v) { return true; }
+bool is_integer(vec2 v) { return false; }
+bool is_integer(vec3 v) { return false; }
+bool is_integer(vec4 v) { return false; }
+bool is_integer(mat2x2 v) { return false; }
+bool is_integer(mat2x3 v) { return false; }
+bool is_integer(mat2x4 v) { return false; }
+bool is_integer(mat3x2 v) { return false; }
+bool is_integer(mat3x3 v) { return false; }
+bool is_integer(mat3x4 v) { return false; }
+bool is_integer(mat4x2 v) { return false; }
+bool is_integer(mat4x3 v) { return false; }
+bool is_integer(mat4x4 v) { return false; }
+
+int mat_row_len(mat2x2 v) { return 2; }
+int mat_row_len(mat2x3 v) { return 3; }
+int mat_row_len(mat2x4 v) { return 4; }
+int mat_row_len(mat3x2 v) { return 2; }
+int mat_row_len(mat3x3 v) { return 3; }
+int mat_row_len(mat3x4 v) { return 4; }
+int mat_row_len(mat4x2 v) { return 2; }
+int mat_row_len(mat4x3 v) { return 3; }
+int mat_row_len(mat4x4 v) { return 4; }
+
+int mat_col_len(mat2x2 v) { return 2; }
+int mat_col_len(mat2x3 v) { return 2; }
+int mat_col_len(mat2x4 v) { return 2; }
+int mat_col_len(mat3x2 v) { return 3; }
+int mat_col_len(mat3x3 v) { return 3; }
+int mat_col_len(mat3x4 v) { return 3; }
+int mat_col_len(mat4x2 v) { return 4; }
+int mat_col_len(mat4x3 v) { return 4; }
+int mat_col_len(mat4x4 v) { return 4; }
+int mat_col_len(ivec2 v) { return 2; }
+int mat_col_len(ivec3 v) { return 3; }
+int mat_col_len(ivec4 v) { return 4; }
+int mat_col_len(uvec2 v) { return 2; }
+int mat_col_len(uvec3 v) { return 3; }
+int mat_col_len(uvec4 v) { return 4; }
+int mat_col_len(vec2 v) { return 2; }
+int mat_col_len(vec3 v) { return 3; }
+int mat_col_len(vec4 v) { return 4; }
+
+#ifndef GPU_METAL
+uint to_type(bool v) { return TEST_TYPE_BOOL; }
+#endif
+uint to_type(uint v) { return TEST_TYPE_UINT; }
+uint to_type(int v) { return TEST_TYPE_INT; }
+uint to_type(float v) { return TEST_TYPE_FLOAT; }
+uint to_type(ivec2 v) { return TEST_TYPE_IVEC2; }
+uint to_type(ivec3 v) { return TEST_TYPE_IVEC3; }
+uint to_type(ivec4 v) { return TEST_TYPE_IVEC4; }
+uint to_type(uvec2 v) { return TEST_TYPE_UVEC2; }
+uint to_type(uvec3 v) { return TEST_TYPE_UVEC3; }
+uint to_type(uvec4 v) { return TEST_TYPE_UVEC4; }
+uint to_type(vec2 v) { return TEST_TYPE_VEC2; }
+uint to_type(vec3 v) { return TEST_TYPE_VEC3; }
+uint to_type(vec4 v) { return TEST_TYPE_VEC4; }
+uint to_type(mat2x2 v) { return TEST_TYPE_MAT2X2; }
+uint to_type(mat2x3 v) { return TEST_TYPE_MAT2X3; }
+uint to_type(mat2x4 v) { return TEST_TYPE_MAT2X4; }
+uint to_type(mat3x2 v) { return TEST_TYPE_MAT3X2; }
+uint to_type(mat3x3 v) { return TEST_TYPE_MAT3X3; }
+uint to_type(mat3x4 v) { return TEST_TYPE_MAT3X4; }
+uint to_type(mat4x2 v) { return TEST_TYPE_MAT4X2; }
+uint to_type(mat4x3 v) { return TEST_TYPE_MAT4X3; }
+uint to_type(mat4x4 v) { return TEST_TYPE_MAT4X4; }
+// clang-format on
+
+#define WRITE_MATRIX(v) \
+  TestOutputRawData raw; \
+  for (int c = 0; c < mat_col_len(v); c++) { \
+    for (int r = 0; r < mat_row_len(v); r++) { \
+      raw.data[c * mat_row_len(v) + r] = floatBitsToUint(v[c][r]); \
+    } \
+  } \
+  return raw;
+
+#define WRITE_FLOAT_VECTOR(v) \
+  TestOutputRawData raw; \
+  for (int c = 0; c < mat_col_len(v); c++) { \
+    raw.data[c] = floatBitsToUint(v[c]); \
+  } \
+  return raw;
+
+#define WRITE_INT_VECTOR(v) \
+  TestOutputRawData raw; \
+  for (int c = 0; c < mat_col_len(v); c++) { \
+    raw.data[c] = uint(v[c]); \
+  } \
+  return raw;
+
+#define WRITE_FLOAT_SCALAR(v) \
+  TestOutputRawData raw; \
+  raw.data[0] = floatBitsToUint(v); \
+  return raw;
+
+#define WRITE_INT_SCALAR(v) \
+  TestOutputRawData raw; \
+  raw.data[0] = uint(v); \
+  return raw;
+
+// clang-format off
+#ifndef GPU_METAL
+TestOutputRawData as_raw_data(bool v) { WRITE_INT_SCALAR(v); }
+#endif
+TestOutputRawData as_raw_data(uint v) { WRITE_INT_SCALAR(v); }
+TestOutputRawData as_raw_data(int v) { WRITE_INT_SCALAR(v); }
+TestOutputRawData as_raw_data(float v) { WRITE_FLOAT_SCALAR(v); }
+TestOutputRawData as_raw_data(ivec2 v) { WRITE_INT_VECTOR(v); }
+TestOutputRawData as_raw_data(ivec3 v) { WRITE_INT_VECTOR(v); }
+TestOutputRawData as_raw_data(ivec4 v) { WRITE_INT_VECTOR(v); }
+TestOutputRawData as_raw_data(uvec2 v) { WRITE_INT_VECTOR(v); }
+TestOutputRawData as_raw_data(uvec3 v) { WRITE_INT_VECTOR(v); }
+TestOutputRawData as_raw_data(uvec4 v) { WRITE_INT_VECTOR(v); }
+TestOutputRawData as_raw_data(vec2 v) { WRITE_FLOAT_VECTOR(v); }
+TestOutputRawData as_raw_data(vec3 v) { WRITE_FLOAT_VECTOR(v); }
+TestOutputRawData as_raw_data(vec4 v) { WRITE_FLOAT_VECTOR(v); }
+TestOutputRawData as_raw_data(mat2x2 v) { WRITE_MATRIX(v); }
+TestOutputRawData as_raw_data(mat2x3 v) { WRITE_MATRIX(v); }
+TestOutputRawData as_raw_data(mat2x4 v) { WRITE_MATRIX(v); }
+TestOutputRawData as_raw_data(mat3x2 v) { WRITE_MATRIX(v); }
+TestOutputRawData as_raw_data(mat3x3 v) { WRITE_MATRIX(v); }
+TestOutputRawData as_raw_data(mat3x4 v) { WRITE_MATRIX(v); }
+TestOutputRawData as_raw_data(mat4x2 v) { WRITE_MATRIX(v); }
+TestOutputRawData as_raw_data(mat4x3 v) { WRITE_MATRIX(v); }
+TestOutputRawData as_raw_data(mat4x4 v) { WRITE_MATRIX(v); }
+// clang-format on
+
+int g_test_id = 0;
+
+#ifdef GPU_COMPUTE_SHADER
+
+#  define EXPECT_OP(OP, val1, val2) \
+    out_test[g_test_id++] = test_output( \
+        as_raw_data(val1), as_raw_data(val2), bool(OP), int(__LINE__), to_type(val1))
+#else
+
+/** WORKAROUND: Fragment shader variant for older platform. */
+#  define EXPECT_OP(OP, val1, val2) \
+    g_test_id += 1; \
+    if (g_test_id == 1) { \
+      /* Avoid pixels with undefined values. */ \
+      out_test = uvec4(0); \
+    } \
+    if (int(gl_FragCoord.y) == g_test_id - 1) { \
+      TestOutput to = test_output( \
+          as_raw_data(val1), as_raw_data(val2), bool(OP), int(__LINE__), to_type(val1)); \
+      switch (int(gl_FragCoord.x)) { \
+        case 0: \
+          out_test = uvec4( \
+              to.expect.data[0], to.expect.data[1], to.expect.data[2], to.expect.data[3]); \
+          break; \
+        case 1: \
+          out_test = uvec4( \
+              to.expect.data[4], to.expect.data[5], to.expect.data[6], to.expect.data[7]); \
+          break; \
+        case 2: \
+          out_test = uvec4( \
+              to.expect.data[8], to.expect.data[9], to.expect.data[10], to.expect.data[11]); \
+          break; \
+        case 3: \
+          out_test = uvec4( \
+              to.expect.data[12], to.expect.data[13], to.expect.data[14], to.expect.data[15]); \
+          break; \
+        case 4: \
+          out_test = uvec4( \
+              to.result.data[0], to.result.data[1], to.result.data[2], to.result.data[3]); \
+          break; \
+        case 5: \
+          out_test = uvec4( \
+              to.result.data[4], to.result.data[5], to.result.data[6], to.result.data[7]); \
+          break; \
+        case 6: \
+          out_test = uvec4( \
+              to.result.data[8], to.result.data[9], to.result.data[10], to.result.data[11]); \
+          break; \
+        case 7: \
+          out_test = uvec4( \
+              to.result.data[12], to.result.data[13], to.result.data[14], to.result.data[15]); \
+          break; \
+        case 8: \
+          out_test = uvec4(to.status, to.line, to.type, 0); \
+          break; \
+      } \
+    }
+
+#endif
+
+#define EXPECT_EQ(result, expect) EXPECT_OP((result) == (expect), result, expect)
+#define EXPECT_NE(result, expect) EXPECT_OP((result) != (expect), result, expect)
+#define EXPECT_LE(result, expect) EXPECT_OP((result) <= (expect), result, expect)
+#define EXPECT_LT(result, expect) EXPECT_OP((result) < (expect), result, expect)
+#define EXPECT_GE(result, expect) EXPECT_OP((result) >= (expect), result, expect)
+#define EXPECT_GT(result, expect) EXPECT_OP((result) > (expect), result, expect)
+
+#define EXPECT_TRUE(result) EXPECT_OP(result, result, true)
+#define EXPECT_FALSE(result) EXPECT_OP(!result, result, false)
+
+#define EXPECT_NEAR(result, expect, threshold) \
+  EXPECT_OP(is_equal(result, expect, threshold), result, expect)

source/blender/gpu/shaders/common/gpu_shader_utildefines_lib.glsl

@@ -0,0 +1,96 @@
+
+/* WORKAROUND: to guard against double include in EEVEE. */
+#ifndef GPU_SHADER_UTILDEFINES_GLSL
+#define GPU_SHADER_UTILDEFINES_GLSL
+
+#ifndef FLT_MAX
+#  define FLT_MAX uintBitsToFloat(0x7F7FFFFFu)
+#  define FLT_MIN uintBitsToFloat(0x00800000u)
+#  define FLT_EPSILON 1.192092896e-07F
+#  define SHRT_MAX 0x00007FFF
+#  define INT_MAX 0x7FFFFFFF
+#  define USHRT_MAX 0x0000FFFFu
+#  define UINT_MAX 0xFFFFFFFFu
+#endif
+#define NAN_FLT uintBitsToFloat(0x7FC00000u)
+
+#define UNPACK2(a) (a)[0], (a)[1]
+#define UNPACK3(a) (a)[0], (a)[1], (a)[2]
+#define UNPACK4(a) (a)[0], (a)[1], (a)[2], (a)[3]
+
+/**
+ * Clamp input into [0..1] range.
+ */
+#define saturate(a) clamp(a, 0.0, 1.0)
+
+#define isfinite(a) (!isinf(a) && !isnan(a))
+
+/* clang-format off */
+#define in_range_inclusive(val, min_v, max_v) (all(greaterThanEqual(val, min_v)) && all(lessThanEqual(val, max_v)))
+#define in_range_exclusive(val, min_v, max_v) (all(greaterThan(val, min_v)) && all(lessThan(val, max_v)))
+#define in_texture_range(texel, tex) (all(greaterThanEqual(texel, ivec2(0))) && all(lessThan(texel, textureSize(tex, 0).xy)))
+#define in_image_range(texel, tex) (all(greaterThanEqual(texel, ivec2(0))) && all(lessThan(texel, imageSize(tex).xy)))
+/* clang-format on */
+
+bool flag_test(uint flag, uint val)
+{
+  return (flag & val) != 0u;
+}
+bool flag_test(int flag, uint val)
+{
+  return flag_test(uint(flag), val);
+}
+bool flag_test(int flag, int val)
+{
+  return (flag & val) != 0;
+}
+
+void set_flag_from_test(inout uint value, bool test, uint flag)
+{
+  if (test) {
+    value |= flag;
+  }
+  else {
+    value &= ~flag;
+  }
+}
+void set_flag_from_test(inout int value, bool test, int flag)
+{
+  if (test) {
+    value |= flag;
+  }
+  else {
+    value &= ~flag;
+  }
+}
+
+/* Keep define to match C++ implementation. */
+#define SET_FLAG_FROM_TEST(value, test, flag) flag_test(value, test, flag)
+
+/**
+ * Pack two 16-bit uint into one 32-bit uint.
+ */
+uint packUvec2x16(uvec2 data)
+{
+  data = (data & 0xFFFFu) << uvec2(0u, 16u);
+  return data.x | data.y;
+}
+uvec2 unpackUvec2x16(uint data)
+{
+  return (uvec2(data) >> uvec2(0u, 16u)) & uvec2(0xFFFFu);
+}
+
+/**
+ * Pack four 8-bit uint into one 32-bit uint.
+ */
+uint packUvec4x8(uvec4 data)
+{
+  data = (data & 0xFFu) << uvec4(0u, 8u, 16u, 24u);
+  return data.x | data.y | data.z | data.w;
+}
+uvec4 unpackUvec4x8(uint data)
+{
+  return (uvec4(data) >> uvec4(0u, 8u, 16u, 24u)) & uvec4(0xFFu);
+}
+
+#endif /* GPU_SHADER_UTILDEFINES_GLSL */

source/blender/gpu/shaders/infos/gpu_shader_test_info.hh

@@ -0,0 +1,19 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later
+ * Copyright 2022 Blender Foundation. All rights reserved. */
+
+/** \file
+ * \ingroup gpu
+ */
+
+#include "gpu_interface_info.hh"
+#include "gpu_shader_create_info.hh"
+
+GPU_SHADER_CREATE_INFO(gpu_shader_test)
+    .typedef_source("GPU_shader_shared.h")
+    .fragment_out(0, Type::UVEC4, "out_test")
+    .additional_info("draw_fullscreen");
+
+GPU_SHADER_CREATE_INFO(gpu_math_test)
+    .fragment_source("gpu_math_test.glsl")
+    .additional_info("gpu_shader_test")
+    .do_static_compilation(true);

source/blender/gpu/tests/gpu_math_test.glsl

@@ -0,0 +1,290 @@
+/* Directive for resetting the line numbering so the failing tests lines can be printed.
+ * This conflict with the shader compiler error logging scheme.
+ * Comment out for correct compilation error line. */
+#line 5
+
+#pragma BLENDER_REQUIRE(gpu_shader_math_matrix_lib.glsl)
+#pragma BLENDER_REQUIRE(gpu_shader_test_lib.glsl)
+
+#define TEST(a, b) if (true)
+
+void main()
+{
+  TEST(math_matrix, MatrixInverse)
+  {
+    mat3x3 mat = mat3x3__diagonal(2);
+    mat3x3 inv = invert(mat);
+    mat3x3 expect = mat3x3__diagonal(0.5f);
+    EXPECT_NEAR(inv, expect, 1e-5f);
+
+    bool success;
+    mat3x3 m2 = mat3x3__all(1);
+    mat3x3 inv2 = invert(m2, success);
+    mat3x3 expect2 = mat3x3__all(0);
+    EXPECT_NEAR(inv2, expect2, 1e-5f);
+    EXPECT_FALSE(success);
+  }
+
+  TEST(math_matrix, MatrixDeterminant)
+  {
+    mat2x2 m2 = mat2x2(vec2(1, 2), vec2(3, 4));
+    mat3x3 m3 = mat3x3(vec3(1, 2, 3), vec3(-3, 4, -5), vec3(5, -6, 7));
+    mat4x4 m4 = mat4x4(vec4(1, 2, -3, 3), vec4(3, 4, -5, 3), vec4(5, 6, 7, -3), vec4(5, 6, 7, 1));
+    EXPECT_NEAR(determinant(m2), -2.0f, 1e-8f);
+    EXPECT_NEAR(determinant(m3), -16.0f, 1e-8f);
+    EXPECT_NEAR(determinant(m4), -112.0f, 1e-8f);
+  }
+
+  TEST(math_matrix, MatrixAdjoint)
+  {
+    mat2x2 m2 = mat2x2(vec2(1, 2), vec2(3, 4));
+    mat3x3 m3 = mat3x3(vec3(1, 2, 3), vec3(-3, 4, -5), vec3(5, -6, 7));
+    mat4x4 m4 = mat4x4(vec4(1, 2, -3, 3), vec4(3, 4, -5, 3), vec4(5, 6, 7, -3), vec4(5, 6, 7, 1));
+    mat2x2 expect2 = transpose(mat2x2(vec2(4, -3), vec2(-2, 1)));
+    mat3x3 expect3 = transpose(mat3x3(vec3(-2, -4, -2), vec3(-32, -8, 16), vec3(-22, -4, 10)));
+    mat4x4 expect4 = transpose(mat4x4(vec4(232, -184, -8, -0),
+                                      vec4(-128, 88, 16, 0),
+                                      vec4(80, -76, 4, 28),
+                                      vec4(-72, 60, -12, -28)));
+    EXPECT_NEAR(adjoint(m2), expect2, 1e-8f);
+    EXPECT_NEAR(adjoint(m3), expect3, 1e-8f);
+    EXPECT_NEAR(adjoint(m4), expect4, 1e-8f);
+  }
+
+  TEST(math_matrix, MatrixInit)
+  {
+    mat4x4 expect;
+
+    mat4x4 m = from_location(vec3(1, 2, 3));
+    expect = mat4x4(vec4(1, 0, 0, 0), vec4(0, 1, 0, 0), vec4(0, 0, 1, 0), vec4(1, 2, 3, 1));
+    EXPECT_TRUE(is_equal(m, expect, 0.00001));
+
+    expect = transpose(mat4x4(vec4(0.411982, -0.833738, -0.36763, 0),
+                              vec4(-0.0587266, -0.426918, 0.902382, 0),
+                              vec4(-0.909297, -0.350175, -0.224845, 0),
+                              vec4(0, 0, 0, 1)));
+    EulerXYZ euler = EulerXYZ(1, 2, 3);
+    Quaternion quat = to_quaternion(euler);
+    AxisAngle axis_angle = to_axis_angle(euler);
+    m = mat4(from_rotation(euler));
+    EXPECT_NEAR(m, expect, 1e-5);
+    m = mat4(from_rotation(quat));
+    EXPECT_NEAR(m, expect, 1e-5);
+    m = mat4(from_rotation(axis_angle));
+    EXPECT_NEAR(m, expect, 1e-5);
+
+    m = from_scale(vec4(1, 2, 3, 4));
+    expect = mat4x4(vec4(1, 0, 0, 0), vec4(0, 2, 0, 0), vec4(0, 0, 3, 0), vec4(0, 0, 0, 4));
+    EXPECT_TRUE(is_equal(m, expect, 0.00001));
+
+    m = mat4(from_scale(vec3(1, 2, 3)));
+    expect = mat4x4(vec4(1, 0, 0, 0), vec4(0, 2, 0, 0), vec4(0, 0, 3, 0), vec4(0, 0, 0, 1));
+    EXPECT_TRUE(is_equal(m, expect, 0.00001));
+
+    m = mat4(from_scale(vec2(1, 2)));
+    expect = mat4x4(vec4(1, 0, 0, 0), vec4(0, 2, 0, 0), vec4(0, 0, 1, 0), vec4(0, 0, 0, 1));
+    EXPECT_TRUE(is_equal(m, expect, 0.00001));
+
+    m = from_loc_rot(vec3(1, 2, 3), EulerXYZ(1, 2, 3));
+    expect = mat4x4(vec4(0.411982, -0.0587266, -0.909297, 0),
+                    vec4(-0.833738, -0.426918, -0.350175, 0),
+                    vec4(-0.36763, 0.902382, -0.224845, 0),
+                    vec4(1, 2, 3, 1));
+    EXPECT_TRUE(is_equal(m, expect, 0.00001));
+
+    m = from_loc_rot_scale(vec3(1, 2, 3), EulerXYZ(1, 2, 3), vec3(1, 2, 3));
+    expect = mat4x4(vec4(0.411982, -0.0587266, -0.909297, 0),
+                    vec4(-1.66748, -0.853835, -0.700351, 0),
+                    vec4(-1.10289, 2.70714, -0.674535, 0),
+                    vec4(1, 2, 3, 1));
+    EXPECT_TRUE(is_equal(m, expect, 0.00001));
+  }
+
+  TEST(math_matrix, MatrixModify)
+  {
+    const float epsilon = 1e-6;
+    mat4x4 result, expect;
+    mat4x4 m1 = mat4x4(vec4(0, 3, 0, 0), vec4(2, 0, 0, 0), vec4(0, 0, 2, 0), vec4(0, 0, 0, 1));
+
+    expect = mat4x4(vec4(0, 3, 0, 0), vec4(2, 0, 0, 0), vec4(0, 0, 2, 0), vec4(4, 9, 2, 1));
+    result = translate(m1, vec3(3, 2, 1));
+    EXPECT_NEAR(result, expect, epsilon);
+
+    expect = mat4x4(vec4(0, 3, 0, 0), vec4(2, 0, 0, 0), vec4(0, 0, 2, 0), vec4(4, 0, 0, 1));
+    result = translate(m1, vec2(0, 2));
+    EXPECT_NEAR(result, expect, epsilon);
+
+    expect = mat4x4(vec4(0, 0, -2, 0), vec4(2, 0, 0, 0), vec4(0, 3, 0, 0), vec4(0, 0, 0, 1));
+    result = rotate(m1, AxisAngle(vec3(0, 1, 0), M_PI_2));
+    EXPECT_NEAR(result, expect, epsilon);
+
+    expect = mat4x4(vec4(0, 9, 0, 0), vec4(4, 0, 0, 0), vec4(0, 0, 8, 0), vec4(0, 0, 0, 1));
+    result = scale(m1, vec3(3, 2, 4));
+    EXPECT_NEAR(result, expect, epsilon);
+
+    expect = mat4x4(vec4(0, 9, 0, 0), vec4(4, 0, 0, 0), vec4(0, 0, 2, 0), vec4(0, 0, 0, 1));
+    result = scale(m1, vec2(3, 2));
+    EXPECT_NEAR(result, expect, epsilon);
+  }
+
+  TEST(math_matrix, MatrixCompareTest)
+  {
+    mat4x4 m1 = mat4x4(vec4(0, 3, 0, 0), vec4(2, 0, 0, 0), vec4(0, 0, 2, 0), vec4(0, 0, 0, 1));
+    mat4x4 m2 = mat4x4(
+        vec4(0, 3.001, 0, 0), vec4(1.999, 0, 0, 0), vec4(0, 0, 2.001, 0), vec4(0, 0, 0, 1.001));
+    mat4x4 m3 = mat4x4(
+        vec4(0, 3.001, 0, 0), vec4(1, 1, 0, 0), vec4(0, 0, 2.001, 0), vec4(0, 0, 0, 1.001));
+    mat4x4 m4 = mat4x4(vec4(0, 1, 0, 0), vec4(1, 0, 0, 0), vec4(0, 0, 1, 0), vec4(0, 0, 0, 1));
+    mat4x4 m5 = mat4x4(vec4(0, 0, 0, 0), vec4(0, 0, 0, 0), vec4(0, 0, 0, 0), vec4(0, 0, 0, 0));
+    mat4x4 m6 = mat4x4(vec4(1, 0, 0, 0), vec4(0, 1, 0, 0), vec4(0, 0, 1, 0), vec4(0, 0, 0, 1));
+    EXPECT_TRUE(is_equal(m1, m2, 0.01f));
+    EXPECT_FALSE(is_equal(m1, m2, 0.0001f));
+    EXPECT_FALSE(is_equal(m1, m3, 0.01f));
+    EXPECT_TRUE(is_orthogonal(m1));
+    EXPECT_FALSE(is_orthogonal(m3));
+    EXPECT_TRUE(is_orthonormal(m4));
+    EXPECT_FALSE(is_orthonormal(m1));
+    EXPECT_FALSE(is_orthonormal(m3));
+    EXPECT_FALSE(is_uniformly_scaled(m1));
+    EXPECT_TRUE(is_uniformly_scaled(m4));
+    EXPECT_FALSE(is_zero(m4));
+    EXPECT_TRUE(is_zero(m5));
+    EXPECT_TRUE(is_negative(m4));
+    EXPECT_FALSE(is_negative(m5));
+    EXPECT_FALSE(is_negative(m6));
+  }
+
+  TEST(math_matrix, MatrixMethods)
+  {
+    mat4x4 m = mat4x4(vec4(0, 3, 0, 0), vec4(2, 0, 0, 0), vec4(0, 0, 2, 0), vec4(0, 1, 0, 1));
+    EulerXYZ expect_eul = EulerXYZ(0, 0, M_PI_2);
+    Quaternion expect_qt = Quaternion(0, -M_SQRT1_2, M_SQRT1_2, 0);
+    vec3 expect_scale = vec3(3, 2, 2);
+    vec3 expect_location = vec3(0, 1, 0);
+
+    EXPECT_NEAR(as_vec3(to_euler(m)), as_vec3(expect_eul), 0.0002);
+    EXPECT_NEAR(as_vec4(to_quaternion(m)), as_vec4(expect_qt), 0.0002);
+    EXPECT_NEAR(to_scale(m), expect_scale, 0.00001);
+
+    vec4 expect_sz = vec4(3, 2, 2, M_SQRT2);
+    vec4 size;
+    mat4x4 m1 = normalize_and_get_size(m, size);
+    EXPECT_TRUE(is_unit_scale(m1));
+    EXPECT_NEAR(size, expect_sz, 0.0002);
+
+    mat4x4 m2 = normalize(m);
+    EXPECT_TRUE(is_unit_scale(m2));
+
+    EulerXYZ eul;
+    Quaternion qt;
+    vec3 scale;
+    to_rot_scale(mat3x3(m), eul, scale);
+    to_rot_scale(mat3x3(m), qt, scale);
+    EXPECT_NEAR(scale, expect_scale, 0.00001);
+    EXPECT_NEAR(as_vec4(qt), as_vec4(expect_qt), 0.0002);
+    EXPECT_NEAR(as_vec3(eul), as_vec3(expect_eul), 0.0002);
+
+    vec3 loc;
+    to_loc_rot_scale(m, loc, eul, scale);
+    to_loc_rot_scale(m, loc, qt, scale);
+    EXPECT_NEAR(scale, expect_scale, 0.00001);
+    EXPECT_NEAR(loc, expect_location, 0.00001);
+    EXPECT_NEAR(as_vec4(qt), as_vec4(expect_qt), 0.0002);
+    EXPECT_NEAR(as_vec3(eul), as_vec3(expect_eul), 0.0002);
+  }
+
+  TEST(math_matrix, MatrixTranspose)
+  {
+    mat4x4 m = mat4x4(vec4(1, 2, 3, 4), vec4(5, 6, 7, 8), vec4(9, 1, 2, 3), vec4(2, 5, 6, 7));
+    mat4x4 expect = mat4x4(vec4(1, 5, 9, 2), vec4(2, 6, 1, 5), vec4(3, 7, 2, 6), vec4(4, 8, 3, 7));
+    EXPECT_EQ(transpose(m), expect);
+  }
+
+  TEST(math_matrix, MatrixInterpolationRegular)
+  {
+    /* Test 4x4 matrix interpolation without singularity, i.e. without axis flip. */
+
+    /* Transposed matrix, so that the code here is written in the same way as print_m4() outputs.
+     */
+    /* This matrix represents T=(0.1, 0.2, 0.3), R=(40, 50, 60) degrees, S=(0.7, 0.8, 0.9) */
+    mat4x4 m2 = transpose(mat4x4(vec4(0.224976f, -0.333770f, 0.765074f, 0.100000f),
+                                 vec4(0.389669f, 0.647565f, 0.168130f, 0.200000f),
+                                 vec4(-0.536231f, 0.330541f, 0.443163f, 0.300000f),
+                                 vec4(0.000000f, 0.000000f, 0.000000f, 1.000000f)));
+    mat4x4 m1 = mat4x4__identity();
+    mat4x4 result;
+    const float epsilon = 1e-6;
+    result = interpolate_fast(m1, m2, 0.0f);
+    EXPECT_NEAR(result, m1, epsilon);
+    result = interpolate_fast(m1, m2, 1.0f);
+    EXPECT_NEAR(result, m2, epsilon);
+
+    /* This matrix is based on the current implementation of the code, and isn't guaranteed to be
+     * correct. It's just consistent with the current implementation. */
+    mat4x4 expect = transpose(mat4x4(vec4(0.690643f, -0.253244f, 0.484996f, 0.050000f),
+                                     vec4(0.271924f, 0.852623f, 0.012348f, 0.100000f),
+                                     vec4(-0.414209f, 0.137484f, 0.816778f, 0.150000f),
+                                     vec4(0.000000f, 0.000000f, 0.000000f, 1.000000f)));
+    result = interpolate_fast(m1, m2, 0.5f);
+    EXPECT_NEAR(result, expect, epsilon);
+  }
+
+  TEST(math_matrix, MatrixTransform)
+  {
+    vec3 expect, result;
+    const vec3 p = vec3(1, 2, 3);
+    mat4x4 m4 = from_loc_rot(vec3(10, 0, 0), EulerXYZ(M_PI_2, M_PI_2, M_PI_2));
+    mat3x3 m3 = from_rotation(EulerXYZ(M_PI_2, M_PI_2, M_PI_2));
+    mat4x4 pers4 = projection__perspective(-0.1f, 0.1f, -0.1f, 0.1f, -0.1f, -1.0f);
+    mat3x3 pers3 = mat3x3(vec3(1, 0, 0.1f), vec3(0, 1, 0.1f), vec3(0, 0.1f, 1));
+
+    expect = vec3(13, 2, -1);
+    result = transform_point(m4, p);
+    EXPECT_NEAR(result, expect, 1e-2);
+
+    expect = vec3(3, 2, -1);
+    result = transform_point(m3, p);
+    EXPECT_NEAR(result, expect, 1e-5);
+
+    result = transform_direction(m4, p);
+    EXPECT_NEAR(result, expect, 1e-5);
+
+    result = transform_direction(m3, p);
+    EXPECT_NEAR(result, expect, 1e-5);
+
+    expect = vec3(-0.5, -1, -1.7222222);
+    result = project_point(pers4, p);
+    EXPECT_NEAR(result, expect, 1e-5);
+
+    vec2 expect2 = vec2(0.76923, 1.61538);
+    vec2 result2 = project_point(pers3, p.xy);
+    EXPECT_NEAR(result2, expect2, 1e-5);
+  }
+
+  TEST(math_matrix, MatrixProjection)
+  {
+    mat4x4 expect;
+    mat4x4 ortho = projection__orthographic(-0.2f, 0.3f, -0.2f, 0.4f, -0.2f, -0.5f);
+    mat4x4 pers1 = projection__perspective(-0.2f, 0.3f, -0.2f, 0.4f, -0.2f, -0.5f);
+    mat4x4 pers2 = projection__perspective_fov(
+        atan(-0.2f), atan(0.3f), atan(-0.2f), atan(0.4f), -0.2f, -0.5f);
+
+    expect = transpose(mat4x4(vec4(4.0f, 0.0f, 0.0f, -0.2f),
+                              vec4(0.0f, 3.33333f, 0.0f, -0.333333f),
+                              vec4(0.0f, 0.0f, 6.66667f, -2.33333f),
+                              vec4(0.0f, 0.0f, 0.0f, 1.0f)));
+    EXPECT_NEAR(ortho, expect, 1e-5);
+
+    expect = transpose(mat4x4(vec4(-0.8f, 0.0f, 0.2f, 0.0f),
+                              vec4(0.0f, -0.666667f, 0.333333f, 0.0f),
+                              vec4(0.0f, 0.0f, -2.33333f, 0.666667f),
+                              vec4(0.0f, 0.0f, -1.0f, 1.0f)));
+    EXPECT_NEAR(pers1, expect, 1e-4);
+
+    expect = transpose(mat4x4(vec4(4.0f, 0.0f, 0.2f, 0.0f),
+                              vec4(0.0f, 3.33333f, 0.333333f, 0.0f),
+                              vec4(0.0f, 0.0f, -2.33333f, 0.666667f),
+                              vec4(0.0f, 0.0f, -1.0f, 1.0f)));
+    EXPECT_NEAR(pers2, expect, 1e-4);
+  }
+}

source/blender/gpu/tests/gpu_shader_test.cc

@@ -2,16 +2,24 @@
 
 #include "testing/testing.h"
 
+#include "BLI_math_matrix_types.hh"
+#include "GPU_batch.h"
+#include "GPU_batch_presets.h"
 #include "GPU_capabilities.h"
 #include "GPU_compute.h"
+#include "GPU_context.h"
+#include "GPU_framebuffer.h"
 #include "GPU_index_buffer.h"
 #include "GPU_shader.h"
+#include "GPU_shader_shared.h"
 #include "GPU_texture.h"
 #include "GPU_vertex_buffer.h"
 #include "GPU_vertex_format.h"
 
 #include "MEM_guardedalloc.h"
 
+#include "gpu_shader_create_info.hh"
+#include "gpu_shader_dependency_private.h"
 #include "gpu_testing.hh"
 
 namespace blender::gpu::tests {
@@ -302,4 +310,206 @@ void main() {
 }
 GPU_TEST(gpu_shader_ssbo_binding)
 
+static void test_gpu_texture_read()
+{
+  GPU_render_begin();
+
+  eGPUTextureUsage usage = GPU_TEXTURE_USAGE_ATTACHMENT | GPU_TEXTURE_USAGE_HOST_READ;
+  GPUTexture *rgba32u = GPU_texture_create_2d_ex("rgba32u", 1, 1, 1, GPU_RGBA32UI, usage, nullptr);
+  GPUTexture *rgba16u = GPU_texture_create_2d_ex("rgba16u", 1, 1, 1, GPU_RGBA16UI, usage, nullptr);
+  GPUTexture *rgba32f = GPU_texture_create_2d_ex("rgba32f", 1, 1, 1, GPU_RGBA32F, usage, nullptr);
+
+  const float4 fcol = {0.0f, 1.3f, -231.0f, 1000.0f};
+  const uint4 ucol = {0, 1, 2, 12223};
+  GPU_texture_clear(rgba32u, GPU_DATA_UINT, ucol);
+  GPU_texture_clear(rgba16u, GPU_DATA_UINT, ucol);
+  GPU_texture_clear(rgba32f, GPU_DATA_FLOAT, fcol);
+
+  GPU_finish();
+
+  uint4 *rgba32u_data = (uint4 *)GPU_texture_read(rgba32u, GPU_DATA_UINT, 0);
+  uint4 *rgba16u_data = (uint4 *)GPU_texture_read(rgba16u, GPU_DATA_UINT, 0);
+  float4 *rgba32f_data = (float4 *)GPU_texture_read(rgba32f, GPU_DATA_FLOAT, 0);
+
+  EXPECT_EQ(ucol, *rgba32u_data);
+  EXPECT_EQ(ucol, *rgba16u_data);
+  EXPECT_EQ(fcol, *rgba32f_data);
+
+  MEM_freeN(rgba32u_data);
+  MEM_freeN(rgba16u_data);
+  MEM_freeN(rgba32f_data);
+
+  GPU_texture_free(rgba32u);
+  GPU_texture_free(rgba16u);
+  GPU_texture_free(rgba32f);
+
+  GPU_render_end();
+}
+GPU_TEST(gpu_texture_read)
+
+static std::string print_test_data(const TestOutputRawData &raw, TestType type)
+{
+  std::stringstream ss;
+  switch (type) {
+    case TEST_TYPE_BOOL:
+    case TEST_TYPE_UINT:
+      ss << *reinterpret_cast<const uint *>(&raw);
+      break;
+    case TEST_TYPE_INT:
+      ss << *reinterpret_cast<const int *>(&raw);
+      break;
+    case TEST_TYPE_FLOAT:
+      ss << *reinterpret_cast<const float *>(&raw);
+      break;
+    case TEST_TYPE_IVEC2:
+      ss << *reinterpret_cast<const int2 *>(&raw);
+      break;
+    case TEST_TYPE_IVEC3:
+      ss << *reinterpret_cast<const int3 *>(&raw);
+      break;
+    case TEST_TYPE_IVEC4:
+      ss << *reinterpret_cast<const int4 *>(&raw);
+      break;
+    case TEST_TYPE_UVEC2:
+      ss << *reinterpret_cast<const uint2 *>(&raw);
+      break;
+    case TEST_TYPE_UVEC3:
+      ss << *reinterpret_cast<const uint3 *>(&raw);
+      break;
+    case TEST_TYPE_UVEC4:
+      ss << *reinterpret_cast<const uint4 *>(&raw);
+      break;
+    case TEST_TYPE_VEC2:
+      ss << *reinterpret_cast<const float2 *>(&raw);
+      break;
+    case TEST_TYPE_VEC3:
+      ss << *reinterpret_cast<const float3 *>(&raw);
+      break;
+    case TEST_TYPE_VEC4:
+      ss << *reinterpret_cast<const float4 *>(&raw);
+      break;
+    case TEST_TYPE_MAT2X2:
+      ss << *reinterpret_cast<const float2x2 *>(&raw);
+      break;
+    case TEST_TYPE_MAT2X3:
+      ss << *reinterpret_cast<const float2x3 *>(&raw);
+      break;
+    case TEST_TYPE_MAT2X4:
+      ss << *reinterpret_cast<const float2x4 *>(&raw);
+      break;
+    case TEST_TYPE_MAT3X2:
+      ss << *reinterpret_cast<const float3x2 *>(&raw);
+      break;
+    case TEST_TYPE_MAT3X3:
+      ss << *reinterpret_cast<const float3x3 *>(&raw);
+      break;
+    case TEST_TYPE_MAT3X4:
+      ss << *reinterpret_cast<const float3x4 *>(&raw);
+      break;
+    case TEST_TYPE_MAT4X2:
+      ss << *reinterpret_cast<const float4x2 *>(&raw);
+      break;
+    case TEST_TYPE_MAT4X3:
+      ss << *reinterpret_cast<const float4x3 *>(&raw);
+      break;
+    case TEST_TYPE_MAT4X4:
+      ss << *reinterpret_cast<const float4x4 *>(&raw);
+      break;
+    default:
+      ss << *reinterpret_cast<const MatBase<uint, 4, 4> *>(&raw);
+      break;
+  }
+  return ss.str();
+}
+
+static StringRef print_test_line(StringRefNull test_src, int64_t test_line)
+{
+  /* Start at line one like the line report scheme. */
+  int64_t line = 1;
+  int64_t last_pos = 0;
+  int64_t pos = 0;
+  while ((pos = test_src.find('\n', pos)) != std::string::npos) {
+    if (line == test_line) {
+      return test_src.substr(last_pos, pos - last_pos);
+    }
+    pos += 1; /* Skip newline */
+    last_pos = pos;
+    line++;
+  }
+  return "";
+}
+
+static void gpu_shader_lib_test(const char *test_src_name)
+{
+  using namespace shader;
+
+  GPU_render_begin();
+
+  ShaderCreateInfo create_info(test_src_name);
+  create_info.fragment_source(test_src_name);
+  create_info.additional_info("gpu_shader_test");
+
+  StringRefNull test_src = gpu_shader_dependency_get_source(test_src_name);
+
+  GPUShader *shader = GPU_shader_create_from_info(
+      reinterpret_cast<GPUShaderCreateInfo *>(&create_info));
+  GPU_shader_bind(shader);
+
+  int test_count = 0;
+  /* Count tests. */
+  int64_t pos = 0;
+  StringRefNull target = "EXPECT_";
+  while ((pos = test_src.find(target, pos)) != std::string::npos) {
+    test_count++;
+    pos += sizeof("EXPECT_");
+  }
+
+  int test_output_px_len = divide_ceil_u(sizeof(TestOutput), 4 * 4);
+
+  eGPUTextureUsage usage = GPU_TEXTURE_USAGE_ATTACHMENT | GPU_TEXTURE_USAGE_HOST_READ;
+  GPUTexture *tex = GPU_texture_create_2d_ex(
+      "tx", test_output_px_len, test_count, 1, GPU_RGBA32UI, usage, nullptr);
+  GPUFrameBuffer *fb = GPU_framebuffer_create("test_fb");
+  GPU_framebuffer_ensure_config(&fb, {GPU_ATTACHMENT_NONE, GPU_ATTACHMENT_TEXTURE(tex)});
+  GPU_framebuffer_bind(fb);
+
+  GPU_batch_draw_advanced(GPU_batch_preset_quad(), 0, 3, 0, 1);
+
+  GPU_finish();
+
+  TestOutput *test_data = (TestOutput *)GPU_texture_read(tex, GPU_DATA_UINT, 0);
+  Span<TestOutput> tests(test_data, test_count);
+
+  for (const TestOutput &test : tests) {
+    if (ELEM(test.status, TEST_STATUS_NONE, TEST_STATUS_PASSED)) {
+      continue;
+    }
+    else if (test.status == TEST_STATUS_FAILED) {
+      ADD_FAILURE_AT(test_src_name, test.line)
+          << "Value of: " << print_test_line(test_src, test.line) << "\n"
+          << "  Actual: " << print_test_data(test.expect, TestType(test.type)) << "\n"
+          << "Expected: " << print_test_data(test.result, TestType(test.type)) << "\n";
+    }
+    else {
+      BLI_assert_unreachable();
+    }
+  }
+
+  MEM_freeN(test_data);
+
+  /* Cleanup. */
+  GPU_shader_unbind();
+  GPU_shader_free(shader);
+  GPU_framebuffer_free(fb);
+  GPU_texture_free(tex);
+
+  GPU_render_end();
+}
+
+static void test_gpu_math_lib()
+{
+  gpu_shader_lib_test("gpu_math_test.glsl");
+}
+GPU_TEST(gpu_math_lib)
+
 }  // namespace blender::gpu::tests