diff options
Diffstat (limited to 'source/blender/blenlib')
-rw-r--r-- | source/blender/blenlib/CMakeLists.txt | 1 | ||||
-rw-r--r-- | source/blender/blenlib/intern/math_solvers.c | 17 | ||||
-rw-r--r-- | source/blender/blenlib/tests/BLI_math_solvers_test.cc | 70 |
3 files changed, 87 insertions, 1 deletions
diff --git a/source/blender/blenlib/CMakeLists.txt b/source/blender/blenlib/CMakeLists.txt index 473cc4d65f3..bcc8c181af4 100644 --- a/source/blender/blenlib/CMakeLists.txt +++ b/source/blender/blenlib/CMakeLists.txt @@ -400,6 +400,7 @@ if(WITH_GTESTS) tests/BLI_math_color_test.cc tests/BLI_math_geom_test.cc tests/BLI_math_matrix_test.cc + tests/BLI_math_solvers_test.cc tests/BLI_math_vector_test.cc tests/BLI_memiter_test.cc tests/BLI_memory_utils_test.cc diff --git a/source/blender/blenlib/intern/math_solvers.c b/source/blender/blenlib/intern/math_solvers.c index cda3d9b66a2..e366d834fc4 100644 --- a/source/blender/blenlib/intern/math_solvers.c +++ b/source/blender/blenlib/intern/math_solvers.c @@ -137,9 +137,24 @@ bool BLI_tridiagonal_solve_cyclic( return false; } + /* Degenerate case not handled correctly by the generic formula. */ + if (count == 1) { + r_x[0] = d[0] / (a[0] + b[0] + c[0]); + + return isfinite(r_x[0]); + } + + /* Degenerate case that works but can be simplified. */ + if (count == 2) { + float a2[2] = {0, a[1] + c[1]}; + float c2[2] = {a[0] + c[0], 0}; + + return BLI_tridiagonal_solve(a2, b, c2, d, r_x, count); + } + + /* If not really cyclic, fall back to the simple solver. */ float a0 = a[0], cN = c[count - 1]; - /* if not really cyclic, fall back to the simple solver */ if (a0 == 0.0f && cN == 0.0f) { return BLI_tridiagonal_solve(a, b, c, d, r_x, count); } diff --git a/source/blender/blenlib/tests/BLI_math_solvers_test.cc b/source/blender/blenlib/tests/BLI_math_solvers_test.cc new file mode 100644 index 00000000000..08cb8571490 --- /dev/null +++ b/source/blender/blenlib/tests/BLI_math_solvers_test.cc @@ -0,0 +1,70 @@ +/* Apache License, Version 2.0 */ + +#include "testing/testing.h" + +#include "BLI_math_solvers.h" + +TEST(math_solvers, Tridiagonal1) +{ + float a[1] = {1}; // ignored + float b[1] = {2}; + float c[1] = {1}; // ignored + float d[1] = {4}; + float x[1]; + + EXPECT_TRUE(BLI_tridiagonal_solve(a, b, c, d, x, 1)); + EXPECT_FLOAT_EQ(x[0], 2); +} + +TEST(math_solvers, Tridiagonal3) +{ + float a[3] = {1, 2, 3}; // 1 ignored + float b[3] = {4, 5, 6}; + float c[3] = {7, 8, 9}; // 9 ignored + float d[3] = {18, 36, 24}; + float x[3]; + + EXPECT_TRUE(BLI_tridiagonal_solve(a, b, c, d, x, 3)); + EXPECT_FLOAT_EQ(x[0], 1); + EXPECT_FLOAT_EQ(x[1], 2); + EXPECT_FLOAT_EQ(x[2], 3); +} + +TEST(math_solvers, CyclicTridiagonal1) +{ + float a[1] = {1}; + float b[1] = {2}; + float c[1] = {1}; + float d[1] = {4}; + float x[1]; + + EXPECT_TRUE(BLI_tridiagonal_solve_cyclic(a, b, c, d, x, 1)); + EXPECT_FLOAT_EQ(x[0], 1); +} + +TEST(math_solvers, CyclicTridiagonal2) +{ + float a[2] = {1, 2}; + float b[2] = {3, 4}; + float c[2] = {5, 6}; + float d[2] = {15, 16}; + float x[2]; + + EXPECT_TRUE(BLI_tridiagonal_solve_cyclic(a, b, c, d, x, 2)); + EXPECT_FLOAT_EQ(x[0], 1); + EXPECT_FLOAT_EQ(x[1], 2); +} + +TEST(math_solvers, CyclicTridiagonal3) +{ + float a[3] = {1, 2, 3}; + float b[3] = {4, 5, 6}; + float c[3] = {7, 8, 9}; + float d[3] = {21, 36, 33}; + float x[3]; + + EXPECT_TRUE(BLI_tridiagonal_solve_cyclic(a, b, c, d, x, 3)); + EXPECT_FLOAT_EQ(x[0], 1); + EXPECT_FLOAT_EQ(x[1], 2); + EXPECT_FLOAT_EQ(x[2], 3); +} |