20 float r_eigen_values[3],
21 float r_eigen_vectors[3][3])
25 if ((m3[0][1] != m3[1][0]) || (m3[0][2] != m3[2][0]) || (m3[1][2] != m3[2][1])) {
31 3, (
const float *)m3, r_eigen_values, (
float *)r_eigen_vectors);
34 void BLI_svd_m3(
const float m3[3][3],
float r_U[3][3],
float r_S[3],
float r_V[3][3])
42 const float *
a,
const float *
b,
const float *
c,
const float *d,
float *r_x,
const int count)
49 double *c1 = (
double *)
MEM_mallocN(bytes * 2,
"tridiagonal_c1d1");
50 double *d1 = c1 +
count;
57 double c_prev, d_prev, x_prev;
61 c1[0] = c_prev = ((
double)
c[0]) /
b[0];
62 d1[0] = d_prev = ((
double)d[0]) /
b[0];
64 for (i = 1; i <
count; i++) {
65 double denum =
b[i] -
a[i] * c_prev;
67 c1[i] = c_prev =
c[i] / denum;
68 d1[i] = d_prev = (d[i] -
a[i] * d_prev) / denum;
74 r_x[--i] = ((
float)x_prev);
77 x_prev = d1[i] - c1[i] * x_prev;
78 r_x[i] = ((
float)x_prev);
87 const float *
a,
const float *
b,
const float *
c,
const float *d,
float *r_x,
const int count)
95 r_x[0] = d[0] / (
a[0] +
b[0] +
c[0]);
102 const float a2[2] = {0,
a[1] +
c[1]};
103 const float c2[2] = {
a[0] +
c[0], 0};
109 float a0 =
a[0], cN =
c[
count - 1];
111 if (a0 == 0.0f && cN == 0.0f) {
115 size_t bytes =
sizeof(
float) * (
unsigned)
count;
116 float *tmp = (
float *)
MEM_mallocN(bytes * 2,
"tridiagonal_ex");
117 float *b2 = tmp +
count;
124 memcpy(b2,
b, bytes);
128 memset(tmp, 0, bytes);
138 float coeff = (r_x[0] + r_x[
count - 1]) / (1.0f + tmp[0] + tmp[
count - 1]);
140 for (
int i = 0; i <
count; i++) {
141 r_x[i] -= coeff * tmp[i];
157 const float x_init[3],
160 float fdelta[3], fdeltav, next_fdeltav;
161 float jacobian[3][3], step[3],
x[3], x_next[3];
167 func_delta(userdata,
x, fdelta);
171 printf(
"START (%g, %g, %g) %g %g\n",
x[0],
x[1],
x[2], fdeltav,
epsilon);
174 for (
int i = 0; i == 0 || (i < max_iterations && fdeltav >
epsilon); i++) {
176 func_jacobian(userdata,
x, jacobian);
186 if (func_correction) {
188 printf(
"%3d * (%g, %g, %g)\n", i, x_next[0], x_next[1], x_next[2]);
191 if (!func_correction(userdata,
x, step, x_next)) {
196 func_delta(userdata, x_next, fdelta);
200 printf(
"%3d ? (%g, %g, %g) %g\n", i, x_next[0], x_next[1], x_next[2], next_fdeltav);
204 while (next_fdeltav > fdeltav && next_fdeltav >
epsilon) {
205 float g0 =
sqrtf(fdeltav), g1 =
sqrtf(next_fdeltav);
206 float g01 = -g0 /
len_v3(step);
207 float det = 2.0f * (g1 - g0 - g01);
208 float l = (det == 0.0f) ? 0.1f : -g01 / det;
214 func_delta(userdata, x_next, fdelta);
218 printf(
"%3d . (%g, %g, %g) %g\n", i, x_next[0], x_next[1], x_next[2], next_fdeltav);
223 fdeltav = next_fdeltav;
226 bool success = (fdeltav <=
epsilon);
229 printf(
"%s (%g, %g, %g) %g\n", success ?
"OK " :
"FAIL",
x[0],
x[1],
x[2], fdeltav);
typedef float(TangentPoint)[2]
bool invert_m3(float R[3][3])
void mul_v3_m3v3(float r[3], const float M[3][3], const float a[3])
bool(* Newton3D_CorrectionFunc)(void *userdata, const float x[3], float step[3], float x_next[3])
void(* Newton3D_JacobianFunc)(void *userdata, const float x[3], float r_jacobian[3][3])
void(* Newton3D_DeltaFunc)(void *userdata, const float x[3], float r_delta[3])
MINLINE float len_squared_v3(const float v[3]) ATTR_WARN_UNUSED_RESULT
MINLINE void sub_v3_v3v3(float r[3], const float a[3], const float b[3])
MINLINE void mul_v3_fl(float r[3], float f)
MINLINE void copy_v3_v3(float r[3], const float a[3])
MINLINE float len_v3(const float a[3]) ATTR_WARN_UNUSED_RESULT
Strict compiler flags for areas of code we want to ensure don't do conversions without us knowing abo...
typedef double(DMatrix)[4][4]
Read Guarded memory(de)allocation.
ATTR_WARN_UNUSED_RESULT const BMLoop * l
bool EIG_self_adjoint_eigen_solve(const int size, const float *matrix, float *r_eigen_values, float *r_eigen_vectors)
void(* MEM_freeN)(void *vmemh)
void *(* MEM_mallocN)(size_t len, const char *str)
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])
Solve a generic f(x) = 0 equation using Newton's method.
void BLI_svd_m3(const float m3[3][3], float r_U[3][3], float r_S[3], float r_V[3][3])
Compute the SVD (Singular Values Decomposition) of given 3D matrix (m3 = USV*).
bool BLI_tridiagonal_solve(const float *a, const float *b, const float *c, const float *d, float *r_x, const int count)
Solve a tridiagonal system of equations:
bool BLI_tridiagonal_solve_cyclic(const float *a, const float *b, const float *c, const float *d, float *r_x, const int count)
Solve a possibly cyclic tridiagonal system using the Sherman-Morrison formula.
bool BLI_eigen_solve_selfadjoint_m3(const float m3[3][3], float r_eigen_values[3], float r_eigen_vectors[3][3])
Compute the eigen values and/or vectors of given 3D symmetric (aka adjoint) matrix.
static const pxr::TfToken b("b", pxr::TfToken::Immortal)
void EIG_svd_square_matrix(const int size, const float *matrix, float *r_U, float *r_S, float *r_V)