diff options
author | Campbell Barton <ideasman42@gmail.com> | 2019-04-17 07:17:24 +0300 |
---|---|---|
committer | Campbell Barton <ideasman42@gmail.com> | 2019-04-17 07:21:24 +0300 |
commit | e12c08e8d170b7ca40f204a5b0423c23a9fbc2c1 (patch) | |
tree | 8cf3453d12edb177a218ef8009357518ec6cab6a /source/blender/blenlib/intern/math_solvers.c | |
parent | b3dabc200a4b0399ec6b81f2ff2730d07b44fcaa (diff) |
ClangFormat: apply to source, most of intern
Apply clang format as proposed in T53211.
For details on usage and instructions for migrating branches
without conflicts, see:
https://wiki.blender.org/wiki/Tools/ClangFormat
Diffstat (limited to 'source/blender/blenlib/intern/math_solvers.c')
-rw-r--r-- | source/blender/blenlib/intern/math_solvers.c | 271 |
1 files changed, 139 insertions, 132 deletions
diff --git a/source/blender/blenlib/intern/math_solvers.c b/source/blender/blenlib/intern/math_solvers.c index 2669ac89b25..a6331aaadd8 100644 --- a/source/blender/blenlib/intern/math_solvers.c +++ b/source/blender/blenlib/intern/math_solvers.c @@ -39,19 +39,19 @@ * \return r_eigen_values the computed eigen values (NULL if not needed). * \return r_eigen_vectors the computed eigen vectors (NULL if not needed). */ -bool BLI_eigen_solve_selfadjoint_m3(const float m3[3][3], float r_eigen_values[3], float r_eigen_vectors[3][3]) +bool BLI_eigen_solve_selfadjoint_m3(const float m3[3][3], + float r_eigen_values[3], + float r_eigen_vectors[3][3]) { #ifndef NDEBUG - /* We must assert given matrix is self-adjoint (i.e. symmetric) */ - if ((m3[0][1] != m3[1][0]) || - (m3[0][2] != m3[2][0]) || - (m3[1][2] != m3[2][1])) - { - BLI_assert(0); - } + /* We must assert given matrix is self-adjoint (i.e. symmetric) */ + if ((m3[0][1] != m3[1][0]) || (m3[0][2] != m3[2][0]) || (m3[1][2] != m3[2][1])) { + BLI_assert(0); + } #endif - return EIG_self_adjoint_eigen_solve(3, (const float *)m3, r_eigen_values, (float *)r_eigen_vectors); + return EIG_self_adjoint_eigen_solve( + 3, (const float *)m3, r_eigen_values, (float *)r_eigen_vectors); } /** @@ -64,7 +64,7 @@ bool BLI_eigen_solve_selfadjoint_m3(const float m3[3][3], float r_eigen_values[3 */ void BLI_svd_m3(const float m3[3][3], float r_U[3][3], float r_S[3], float r_V[3][3]) { - EIG_svd_square_matrix(3, (const float *)m3, (float *)r_U, (float *)r_S, (float *)r_V); + EIG_svd_square_matrix(3, (const float *)m3, (float *)r_U, (float *)r_S, (float *)r_V); } /***************************** Simple Solvers ************************************/ @@ -79,48 +79,49 @@ void BLI_svd_m3(const float m3[3][3], float r_U[3][3], float r_S[3], float r_V[3 * \param r_x: output vector, may be shared with any of the input ones * \return true if success */ -bool BLI_tridiagonal_solve(const float *a, const float *b, const float *c, const float *d, float *r_x, const int count) +bool BLI_tridiagonal_solve( + const float *a, const float *b, const float *c, const float *d, float *r_x, const int count) { - if (count < 1) { - return false; - } + if (count < 1) { + return false; + } - size_t bytes = sizeof(double) * (unsigned)count; - double *c1 = (double *)MEM_mallocN(bytes * 2, "tridiagonal_c1d1"); - double *d1 = c1 + count; + size_t bytes = sizeof(double) * (unsigned)count; + double *c1 = (double *)MEM_mallocN(bytes * 2, "tridiagonal_c1d1"); + double *d1 = c1 + count; - if (!c1) { - return false; - } + if (!c1) { + return false; + } - int i; - double c_prev, d_prev, x_prev; + int i; + double c_prev, d_prev, x_prev; - /* forward pass */ + /* forward pass */ - c1[0] = c_prev = ((double)c[0]) / b[0]; - d1[0] = d_prev = ((double)d[0]) / b[0]; + c1[0] = c_prev = ((double)c[0]) / b[0]; + d1[0] = d_prev = ((double)d[0]) / b[0]; - for (i = 1; i < count; i++) { - double denum = b[i] - a[i] * c_prev; + for (i = 1; i < count; i++) { + double denum = b[i] - a[i] * c_prev; - c1[i] = c_prev = c[i] / denum; - d1[i] = d_prev = (d[i] - a[i] * d_prev) / denum; - } + c1[i] = c_prev = c[i] / denum; + d1[i] = d_prev = (d[i] - a[i] * d_prev) / denum; + } - /* back pass */ + /* back pass */ - x_prev = d_prev; - r_x[--i] = ((float)x_prev); + x_prev = d_prev; + r_x[--i] = ((float)x_prev); - while (--i >= 0) { - x_prev = d1[i] - c1[i] * x_prev; - r_x[i] = ((float)x_prev); - } + while (--i >= 0) { + x_prev = d1[i] - c1[i] * x_prev; + r_x[i] = ((float)x_prev); + } - MEM_freeN(c1); + MEM_freeN(c1); - return isfinite(x_prev); + return isfinite(x_prev); } /** @@ -129,53 +130,53 @@ bool BLI_tridiagonal_solve(const float *a, const float *b, const float *c, const * \param r_x: output vector, may be shared with any of the input ones * \return true if success */ -bool BLI_tridiagonal_solve_cyclic(const float *a, const float *b, const float *c, const float *d, float *r_x, const int count) +bool BLI_tridiagonal_solve_cyclic( + const float *a, const float *b, const float *c, const float *d, float *r_x, const int count) { - if (count < 1) { - return false; - } + if (count < 1) { + return false; + } - float a0 = a[0], cN = c[count - 1]; + 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); - } + /* 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); + } - size_t bytes = sizeof(float) * (unsigned)count; - float *tmp = (float *)MEM_mallocN(bytes * 2, "tridiagonal_ex"); - float *b2 = tmp + count; + size_t bytes = sizeof(float) * (unsigned)count; + float *tmp = (float *)MEM_mallocN(bytes * 2, "tridiagonal_ex"); + float *b2 = tmp + count; - if (!tmp) { - return false; - } + if (!tmp) { + return false; + } - /* prepare the noncyclic system; relies on tridiagonal_solve ignoring values */ - memcpy(b2, b, bytes); - b2[0] -= a0; - b2[count - 1] -= cN; + /* prepare the noncyclic system; relies on tridiagonal_solve ignoring values */ + memcpy(b2, b, bytes); + b2[0] -= a0; + b2[count - 1] -= cN; - memset(tmp, 0, bytes); - tmp[0] = a0; - tmp[count - 1] = cN; + memset(tmp, 0, bytes); + tmp[0] = a0; + tmp[count - 1] = cN; - /* solve for partial solution and adjustment vector */ - bool success = - BLI_tridiagonal_solve(a, b2, c, tmp, tmp, count) && - BLI_tridiagonal_solve(a, b2, c, d, r_x, count); + /* solve for partial solution and adjustment vector */ + bool success = BLI_tridiagonal_solve(a, b2, c, tmp, tmp, count) && + BLI_tridiagonal_solve(a, b2, c, d, r_x, count); - /* apply adjustment */ - if (success) { - float coeff = (r_x[0] + r_x[count - 1]) / (1.0f + tmp[0] + tmp[count - 1]); + /* apply adjustment */ + if (success) { + float coeff = (r_x[0] + r_x[count - 1]) / (1.0f + tmp[0] + tmp[count - 1]); - for (int i = 0; i < count; i++) { - r_x[i] -= coeff * tmp[i]; - } - } + for (int i = 0; i < count; i++) { + r_x[i] -= coeff * tmp[i]; + } + } - MEM_freeN(tmp); + MEM_freeN(tmp); - return success; + return success; } /** @@ -192,82 +193,88 @@ bool BLI_tridiagonal_solve_cyclic(const float *a, const float *b, const float *c * \param result: Final result. * \return true if success */ -bool BLI_newton3d_solve( - Newton3D_DeltaFunc func_delta, Newton3D_JacobianFunc func_jacobian, Newton3D_CorrectionFunc func_correction, void *userdata, - float epsilon, int max_iterations, bool trace, const float x_init[3], float result[3]) +bool BLI_newton3d_solve(Newton3D_DeltaFunc func_delta, + Newton3D_JacobianFunc func_jacobian, + Newton3D_CorrectionFunc func_correction, + void *userdata, + float epsilon, + int max_iterations, + bool trace, + const float x_init[3], + float result[3]) { - float fdelta[3], fdeltav, next_fdeltav; - float jacobian[3][3], step[3], x[3], x_next[3]; + float fdelta[3], fdeltav, next_fdeltav; + float jacobian[3][3], step[3], x[3], x_next[3]; - epsilon *= epsilon; + epsilon *= epsilon; - copy_v3_v3(x, x_init); + copy_v3_v3(x, x_init); - func_delta(userdata, x, fdelta); - fdeltav = len_squared_v3(fdelta); + func_delta(userdata, x, fdelta); + fdeltav = len_squared_v3(fdelta); - if (trace) { - printf("START (%g, %g, %g) %g\n", x[0], x[1], x[2], fdeltav); - } + if (trace) { + printf("START (%g, %g, %g) %g\n", x[0], x[1], x[2], fdeltav); + } - for (int i = 0; i < max_iterations && fdeltav > epsilon; i++) { - /* Newton's method step. */ - func_jacobian(userdata, x, jacobian); + for (int i = 0; i < max_iterations && fdeltav > epsilon; i++) { + /* Newton's method step. */ + func_jacobian(userdata, x, jacobian); - if (!invert_m3(jacobian)) { - return false; - } + if (!invert_m3(jacobian)) { + return false; + } - mul_v3_m3v3(step, jacobian, fdelta); - sub_v3_v3v3(x_next, x, step); + mul_v3_m3v3(step, jacobian, fdelta); + sub_v3_v3v3(x_next, x, step); - /* Custom out-of-bounds value correction. */ - if (func_correction) { - if (trace) { - printf("%3d * (%g, %g, %g)\n", i, x_next[0], x_next[1], x_next[2]); - } + /* Custom out-of-bounds value correction. */ + if (func_correction) { + if (trace) { + printf("%3d * (%g, %g, %g)\n", i, x_next[0], x_next[1], x_next[2]); + } - if (!func_correction(userdata, x, step, x_next)) { - return false; - } - } + if (!func_correction(userdata, x, step, x_next)) { + return false; + } + } - func_delta(userdata, x_next, fdelta); - next_fdeltav = len_squared_v3(fdelta); + func_delta(userdata, x_next, fdelta); + next_fdeltav = len_squared_v3(fdelta); - if (trace) { - printf("%3d ? (%g, %g, %g) %g\n", i, x_next[0], x_next[1], x_next[2], next_fdeltav); - } + if (trace) { + printf("%3d ? (%g, %g, %g) %g\n", i, x_next[0], x_next[1], x_next[2], next_fdeltav); + } - /* Line search correction. */ - while (next_fdeltav > fdeltav) { - float g0 = sqrtf(fdeltav), g1 = sqrtf(next_fdeltav); - float g01 = -g0 / len_v3(step); - float det = 2.0f * (g1 - g0 - g01); - float l = (det == 0.0f) ? 0.1f : -g01 / det; - CLAMP_MIN(l, 0.1f); + /* Line search correction. */ + while (next_fdeltav > fdeltav) { + float g0 = sqrtf(fdeltav), g1 = sqrtf(next_fdeltav); + float g01 = -g0 / len_v3(step); + float det = 2.0f * (g1 - g0 - g01); + float l = (det == 0.0f) ? 0.1f : -g01 / det; + CLAMP_MIN(l, 0.1f); - mul_v3_fl(step, l); - sub_v3_v3v3(x_next, x, step); + mul_v3_fl(step, l); + sub_v3_v3v3(x_next, x, step); - func_delta(userdata, x_next, fdelta); - next_fdeltav = len_squared_v3(fdelta); + func_delta(userdata, x_next, fdelta); + next_fdeltav = len_squared_v3(fdelta); - if (trace) { - printf("%3d . (%g, %g, %g) %g\n", i, x_next[0], x_next[1], x_next[2], next_fdeltav); - } - } + if (trace) { + printf("%3d . (%g, %g, %g) %g\n", i, x_next[0], x_next[1], x_next[2], next_fdeltav); + } + } - copy_v3_v3(x, x_next); - fdeltav = next_fdeltav; - } + copy_v3_v3(x, x_next); + fdeltav = next_fdeltav; + } - bool success = (fdeltav <= epsilon); + bool success = (fdeltav <= epsilon); - if (trace) { - printf("%s (%g, %g, %g) %g\n", success ? "OK " : "FAIL", x[0], x[1], x[2], fdeltav); - } + if (trace) { + printf("%s (%g, %g, %g) %g\n", success ? "OK " : "FAIL", x[0], x[1], x[2], fdeltav); + } - copy_v3_v3(result, x); - return success; + copy_v3_v3(result, x); + return success; } |