Fix T53410: 3D Text always recalculated
[blender.git] / intern / cycles / util / util_math_matrix.h
1 /*
2  * Copyright 2011-2017 Blender Foundation
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16
17 #ifndef __UTIL_MATH_MATRIX_H__
18 #define __UTIL_MATH_MATRIX_H__
19
20 CCL_NAMESPACE_BEGIN
21
22 #define MAT(A, size, row, col) A[(row)*(size)+(col)]
23
24 /* Variants that use a constant stride on GPUS. */
25 #ifdef __KERNEL_GPU__
26 #  define MATS(A, n, r, c, s) A[((r)*(n)+(c))*(s)]
27 /* Element access when only the lower-triangular elements are stored. */
28 #  define MATHS(A, r, c, s) A[((r)*((r)+1)/2+(c))*(s)]
29 #  define VECS(V, i, s) V[(i)*(s)]
30 #else
31 #  define MATS(A, n, r, c, s) MAT(A, n, r, c)
32 #  define MATHS(A, r, c, s) A[(r)*((r)+1)/2+(c)]
33 #  define VECS(V, i, s) V[i]
34 #endif
35
36 /* Zeroing helpers. */
37
38 ccl_device_inline void math_vector_zero(float *v, int n)
39 {
40         for(int i = 0; i < n; i++) {
41                 v[i] = 0.0f;
42         }
43 }
44
45 ccl_device_inline void math_matrix_zero(float *A, int n)
46 {
47         for(int row = 0; row < n; row++) {
48                 for(int col = 0; col <= row; col++) {
49                         MAT(A, n, row, col) = 0.0f;
50                 }
51         }
52 }
53
54 /* Elementary vector operations. */
55
56 ccl_device_inline void math_vector_add(float *a, const float *ccl_restrict b, int n)
57 {
58         for(int i = 0; i < n; i++) {
59                 a[i] += b[i];
60         }
61 }
62
63 ccl_device_inline void math_vector_mul(float *a, const float *ccl_restrict b, int n)
64 {
65         for(int i = 0; i < n; i++) {
66                 a[i] *= b[i];
67         }
68 }
69
70 ccl_device_inline void math_vector_mul_strided(ccl_global float *a, const float *ccl_restrict b, int astride, int n)
71 {
72         for(int i = 0; i < n; i++) {
73                 a[i*astride] *= b[i];
74         }
75 }
76
77 ccl_device_inline void math_vector_scale(float *a, float b, int n)
78 {
79         for(int i = 0; i < n; i++) {
80                 a[i] *= b;
81         }
82 }
83
84 ccl_device_inline void math_vector_max(float *a, const float *ccl_restrict b, int n)
85 {
86         for(int i = 0; i < n; i++) {
87                 a[i] = max(a[i], b[i]);
88         }
89 }
90
91 ccl_device_inline void math_vec3_add(float3 *v, int n, float *x, float3 w)
92 {
93         for(int i = 0; i < n; i++) {
94                 v[i] += w*x[i];
95         }
96 }
97
98 ccl_device_inline void math_vec3_add_strided(ccl_global float3 *v, int n, float *x, float3 w, int stride)
99 {
100         for(int i = 0; i < n; i++) {
101                 v[i*stride] += w*x[i];
102         }
103 }
104
105 /* Elementary matrix operations.
106  * Note: TriMatrix refers to a square matrix that is symmetric, and therefore its upper-triangular part isn't stored. */
107
108 ccl_device_inline void math_trimatrix_add_diagonal(ccl_global float *A, int n, float val, int stride)
109 {
110         for(int row = 0; row < n; row++) {
111                 MATHS(A, row, row, stride) += val;
112         }
113 }
114
115 /* Add Gramian matrix of v to A.
116  * The Gramian matrix of v is vt*v, so element (i,j) is v[i]*v[j]. */
117 ccl_device_inline void math_matrix_add_gramian(float *A,
118                                                   int n,
119                                                   const float *ccl_restrict v,
120                                                   float weight)
121 {
122         for(int row = 0; row < n; row++) {
123                 for(int col = 0; col <= row; col++) {
124                         MAT(A, n, row, col) += v[row]*v[col]*weight;
125                 }
126         }
127 }
128
129 /* Add Gramian matrix of v to A.
130  * The Gramian matrix of v is vt*v, so element (i,j) is v[i]*v[j]. */
131 ccl_device_inline void math_trimatrix_add_gramian_strided(ccl_global float *A,
132                                                           int n,
133                                                           const float *ccl_restrict v,
134                                                           float weight,
135                                                           int stride)
136 {
137         for(int row = 0; row < n; row++) {
138                 for(int col = 0; col <= row; col++) {
139                         MATHS(A, row, col, stride) += v[row]*v[col]*weight;
140                 }
141         }
142 }
143
144 /* Transpose matrix A inplace. */
145 ccl_device_inline void math_matrix_transpose(ccl_global float *A, int n, int stride)
146 {
147         for(int i = 0; i < n; i++) {
148                 for(int j = 0; j < i; j++) {
149                         float temp = MATS(A, n, i, j, stride);
150                         MATS(A, n, i, j, stride) = MATS(A, n, j, i, stride);
151                         MATS(A, n, j, i, stride) = temp;
152                 }
153         }
154 }
155
156 /* Solvers for matrix problems */
157
158 /* In-place Cholesky-Banachiewicz decomposition of the square, positive-definite matrix A
159  * into a lower triangular matrix L so that A = L*L^T. A is being overwritten by L.
160  * Also, only the lower triangular part of A is ever accessed. */
161 ccl_device void math_trimatrix_cholesky(ccl_global float *A, int n, int stride)
162 {
163         for(int row = 0; row < n; row++) {
164                 for(int col = 0; col <= row; col++) {
165                         float sum_col = MATHS(A, row, col, stride);
166                         for(int k = 0; k < col; k++) {
167                                 sum_col -= MATHS(A, row, k, stride) * MATHS(A, col, k, stride);
168                         }
169                         if(row == col) {
170                                 sum_col = sqrtf(max(sum_col, 0.0f));
171                         }
172                         else {
173                                 sum_col /= MATHS(A, col, col, stride);
174                         }
175                         MATHS(A, row, col, stride) = sum_col;
176                 }
177         }
178 }
179
180 /* Solve A*S=y for S given A and y, where A is symmetrical positive-semidefinite and both inputs are destroyed in the process.
181  *
182  * We can apply Cholesky decomposition to find a lower triangular L so that L*Lt = A.
183  * With that we get (L*Lt)*S = L*(Lt*S) = L*b = y, defining b as Lt*S.
184  * Since L is lower triangular, finding b is relatively easy since y is known.
185  * Then, the remaining problem is Lt*S = b, which again can be solved easily.
186  *
187  * This is useful for solving the normal equation S=inv(Xt*W*X)*Xt*W*y, since Xt*W*X is
188  * symmetrical positive-semidefinite by construction, so we can just use this function with A=Xt*W*X and y=Xt*W*y. */
189 ccl_device_inline void math_trimatrix_vec3_solve(ccl_global float *A, ccl_global float3 *y, int n, int stride)
190 {
191         /* Since the first entry of the design row is always 1, the upper-left element of XtWX is a good
192          * heuristic for the amount of pixels considered (with weighting), therefore the amount of correction
193          * is scaled based on it. */
194         math_trimatrix_add_diagonal(A, n, 3e-7f*A[0], stride); /* Improve the numerical stability. */
195         math_trimatrix_cholesky(A, n, stride); /* Replace A with L so that L*Lt = A. */
196
197         /* Use forward substitution to solve L*b = y, replacing y by b. */
198         for(int row = 0; row < n; row++) {
199                 float3 sum = VECS(y, row, stride);
200                 for(int col = 0; col < row; col++)
201                         sum -= MATHS(A, row, col, stride) * VECS(y, col, stride);
202                 VECS(y, row, stride) = sum / MATHS(A, row, row, stride);
203         }
204
205         /* Use backward substitution to solve Lt*S = b, replacing b by S. */
206         for(int row = n-1; row >= 0; row--) {
207                 float3 sum = VECS(y, row, stride);
208                 for(int col = row+1; col < n; col++)
209                         sum -= MATHS(A, col, row, stride) * VECS(y, col, stride);
210                 VECS(y, row, stride) = sum / MATHS(A, row, row, stride);
211         }
212 }
213
214 /* Perform the Jacobi Eigenvalue Methon on matrix A.
215  * A is assumed to be a symmetrical matrix, therefore only the lower-triangular part is ever accessed.
216  * The algorithm overwrites the contents of A.
217  *
218  * After returning, A will be overwritten with D, which is (almost) diagonal,
219  * and V will contain the eigenvectors of the original A in its rows (!),
220  * so that A = V^T*D*V. Therefore, the diagonal elements of D are the (sorted) eigenvalues of A.
221  */
222 ccl_device void math_matrix_jacobi_eigendecomposition(float *A, ccl_global float *V, int n, int v_stride)
223 {
224         const float singular_epsilon = 1e-9f;
225
226         for(int row = 0; row < n; row++) {
227                 for(int col = 0; col < n; col++) {
228                         MATS(V, n, row, col, v_stride) = (col == row) ? 1.0f : 0.0f;
229                 }
230         }
231
232         for(int sweep = 0; sweep < 8; sweep++) {
233                 float off_diagonal = 0.0f;
234                 for(int row = 1; row < n; row++) {
235                         for(int col = 0; col < row; col++) {
236                                 off_diagonal += fabsf(MAT(A, n, row, col));
237                         }
238                 }
239                 if(off_diagonal < 1e-7f) {
240                         /* The matrix has nearly reached diagonal form.
241                          * Since the eigenvalues are only used to determine truncation, their exact values aren't required - a relative error of a few ULPs won't matter at all. */
242                         break;
243                 }
244
245                 /* Set the threshold for the small element rotation skip in the first sweep:
246                  * Skip all elements that are less than a tenth of the average off-diagonal element. */
247                 float threshold = 0.2f*off_diagonal / (n*n);
248
249                 for(int row = 1; row < n; row++) {
250                         for(int col = 0; col < row; col++) {
251                                 /* Perform a Jacobi rotation on this element that reduces it to zero. */
252                                 float element = MAT(A, n, row, col);
253                                 float abs_element = fabsf(element);
254
255                                 /* If we're in a later sweep and the element already is very small, just set it to zero and skip the rotation. */
256                                 if(sweep > 3 && abs_element <= singular_epsilon*fabsf(MAT(A, n, row, row)) && abs_element <= singular_epsilon*fabsf(MAT(A, n, col, col))) {
257                                         MAT(A, n, row, col) = 0.0f;
258                                         continue;
259                                 }
260
261                                 if(element == 0.0f) {
262                                         continue;
263                                 }
264
265                                 /* If we're in one of the first sweeps and the element is smaller than the threshold, skip it. */
266                                 if(sweep < 3 && (abs_element < threshold)) {
267                                         continue;
268                                 }
269
270                                 /* Determine rotation: The rotation is characterized by its angle phi - or, in the actual implementation, sin(phi) and cos(phi).
271                                  * To find those, we first compute their ratio - that might be unstable if the angle approaches 90°, so there's a fallback for that case.
272                                  * Then, we compute sin(phi) and cos(phi) themselves. */
273                                 float singular_diff = MAT(A, n, row, row) - MAT(A, n, col, col);
274                                 float ratio;
275                                 if(abs_element > singular_epsilon*fabsf(singular_diff)) {
276                                         float cot_2phi = 0.5f*singular_diff / element;
277                                         ratio = 1.0f / (fabsf(cot_2phi) + sqrtf(1.0f + cot_2phi*cot_2phi));
278                                         if(cot_2phi < 0.0f) ratio = -ratio; /* Copy sign. */
279                                 }
280                                 else {
281                                         ratio = element / singular_diff;
282                                 }
283
284                                 float c = 1.0f / sqrtf(1.0f + ratio*ratio);
285                                 float s = ratio*c;
286                                 /* To improve numerical stability by avoiding cancellation, the update equations are reformulized to use sin(phi) and tan(phi/2) instead. */
287                                 float tan_phi_2 = s / (1.0f + c);
288
289                                 /* Update the singular values in the diagonal. */
290                                 float singular_delta = ratio*element;
291                                 MAT(A, n, row, row) += singular_delta;
292                                 MAT(A, n, col, col) -= singular_delta;
293
294                                 /* Set the element itself to zero. */
295                                 MAT(A, n, row, col) = 0.0f;
296
297                                 /* Perform the actual rotations on the matrices. */
298 #define ROT(M, r1, c1, r2, c2, stride)                                   \
299                                 {                                                        \
300                                         float M1 = MATS(M, n, r1, c1, stride);               \
301                                         float M2 = MATS(M, n, r2, c2, stride);               \
302                                         MATS(M, n, r1, c1, stride) -= s*(M2 + tan_phi_2*M1); \
303                                         MATS(M, n, r2, c2, stride) += s*(M1 - tan_phi_2*M2); \
304                                 }
305
306                                 /* Split into three parts to ensure correct accesses since we only store the lower-triangular part of A. */
307                                 for(int i = 0    ; i < col; i++) ROT(A, col, i, row, i, 1);
308                                 for(int i = col+1; i < row; i++) ROT(A, i, col, row, i, 1);
309                                 for(int i = row+1; i < n  ; i++) ROT(A, i, col, i, row, 1);
310
311                                 for(int i = 0    ; i < n  ; i++) ROT(V, col, i, row, i, v_stride);
312 #undef ROT
313                         }
314                 }
315         }
316
317         /* Sort eigenvalues and the associated eigenvectors. */
318         for(int i = 0; i < n - 1; i++) {
319                 float v = MAT(A, n, i, i);
320                 int k = i;
321                 for(int j = i; j < n; j++) {
322                         if(MAT(A, n, j, j) >= v) {
323                                 v = MAT(A, n, j, j);
324                                 k = j;
325                         }
326                 }
327                 if(k != i) {
328                         /* Swap eigenvalues. */
329                         MAT(A, n, k, k) = MAT(A, n, i, i);
330                         MAT(A, n, i, i) = v;
331                         /* Swap eigenvectors. */
332                         for(int j = 0; j < n; j++) {
333                                 float v = MATS(V, n, i, j, v_stride);
334                                 MATS(V, n, i, j, v_stride) = MATS(V, n, k, j, v_stride);
335                                 MATS(V, n, k, j, v_stride) = v;
336                         }
337                 }
338         }
339 }
340
341 #ifdef __KERNEL_SSE3__
342 ccl_device_inline void math_vector_zero_sse(float4 *A, int n)
343 {
344         for(int i = 0; i < n; i++) {
345                 A[i] = make_float4(0.0f);
346         }
347 }
348
349 ccl_device_inline void math_matrix_zero_sse(float4 *A, int n)
350 {
351         for(int row = 0; row < n; row++) {
352                 for(int col = 0; col <= row; col++) {
353                         MAT(A, n, row, col) = make_float4(0.0f);
354                 }
355         }
356 }
357
358 /* Add Gramian matrix of v to A.
359  * The Gramian matrix of v is v^T*v, so element (i,j) is v[i]*v[j]. */
360 ccl_device_inline void math_matrix_add_gramian_sse(float4 *A, int n, const float4 *ccl_restrict v, float4 weight)
361 {
362         for(int row = 0; row < n; row++) {
363                 for(int col = 0; col <= row; col++) {
364                         MAT(A, n, row, col) = MAT(A, n, row, col) + v[row] * v[col] * weight;
365                 }
366         }
367 }
368
369 ccl_device_inline void math_vector_add_sse(float4 *V, int n, const float4 *ccl_restrict a)
370 {
371         for(int i = 0; i < n; i++) {
372                 V[i] += a[i];
373         }
374 }
375
376 ccl_device_inline void math_vector_mul_sse(float4 *V, int n, const float4 *ccl_restrict a)
377 {
378         for(int i = 0; i < n; i++) {
379                 V[i] *= a[i];
380         }
381 }
382
383 ccl_device_inline void math_vector_max_sse(float4 *a, const float4 *ccl_restrict b, int n)
384 {
385         for(int i = 0; i < n; i++) {
386                 a[i] = max(a[i], b[i]);
387         }
388 }
389
390 ccl_device_inline void math_matrix_hsum(float *A, int n, const float4 *ccl_restrict B)
391 {
392         for(int row = 0; row < n; row++) {
393                 for(int col = 0; col <= row; col++) {
394                         MAT(A, n, row, col) = reduce_add(MAT(B, n, row, col))[0];
395                 }
396         }
397 }
398 #endif
399
400 #undef MAT
401
402 CCL_NAMESPACE_END
403
404 #endif  /* __UTIL_MATH_MATRIX_H__ */