Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2009 Sandia National Laboratories. Developed at the
5 : : University of Wisconsin--Madison under SNL contract number
6 : : 624796. The U.S. Government and the University of Wisconsin
7 : : retain certain rights to this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : (2009) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file TMPDerivs.hpp
28 : : * \brief Utility methods for calculating derivatives of TMP metrics
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #ifndef MSQ_TMP_DERIVS_HPP
33 : : #define MSQ_TMP_DERIVS_HPP
34 : :
35 : : #include "Mesquite.hpp"
36 : : #include "MsqMatrix.hpp"
37 : :
38 : : /** If defined, then \f$ (A \otimes B)_{i,j} = row(A,i)^t \otimes row(B,j) \f$
39 : : otherwise \f$ (A \otimes B)_{i,j} = column(A,i) \otimes column(B,j)^t \f$
40 : : */
41 : : #define MSQ_ROW_BASED_OUTER_PRODUCT
42 : :
43 : : namespace MBMesquite
44 : : {
45 : :
46 : : /**\brief \f$ R *= s \f$ */
47 : : template < unsigned D >
48 : : inline void hess_scale_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha );
49 : :
50 : 0 : inline void hess_scale( MsqMatrix< 3, 3 > R[6], double alpha )
51 : : {
52 : 0 : hess_scale_t< 3 >( R, alpha );
53 : 0 : }
54 : :
55 : 0 : inline void hess_scale( MsqMatrix< 2, 2 > R[3], double alpha )
56 : : {
57 : 0 : hess_scale_t< 2 >( R, alpha );
58 : 0 : }
59 : :
60 : : /**\brief \f$ R = \alpha I_9 \f$
61 : : *
62 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
63 : : * symmetric matrix.
64 : : */
65 : : inline void set_scaled_I( MsqMatrix< 3, 3 > R[6], double alpha );
66 : :
67 : : /**\brief \f$ R += \alpha I_9 \f$
68 : : *
69 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
70 : : * symmetric matrix.
71 : : */
72 : : inline void pluseq_scaled_I( MsqMatrix< 3, 3 > R[6], double alpha );
73 : :
74 : : /**\brief \f$ R += \alpha I_9 \f$
75 : : *
76 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
77 : : * symmetric matrix.
78 : : */
79 : : inline void pluseq_scaled_I( MsqMatrix< 3, 3 >& R, double alpha );
80 : :
81 : : /**\brief \f$ R = \alpha I_4 \f$
82 : : *
83 : : *\param R The 3 blocks of the upper triangular portion of a 4x4
84 : : * symmetric matrix.
85 : : */
86 : : inline void set_scaled_I( MsqMatrix< 2, 2 > R[3], double alpha );
87 : :
88 : : /**\brief \f$ R += \alpha I_4 \f$
89 : : *
90 : : *\param R The 3 blocks of the upper triangular portion of a 4x4
91 : : * symmetric matrix.
92 : : */
93 : : inline void pluseq_scaled_I( MsqMatrix< 2, 2 > R[3], double alpha );
94 : :
95 : : /**\brief \f$ R += \alpha I_4 \f$
96 : : *
97 : : *\param R The 3 blocks of the upper triangular portion of a 4x4
98 : : * symmetric matrix.
99 : : */
100 : : inline void pluseq_scaled_I( MsqMatrix< 2, 2 >& R, double alpha );
101 : :
102 : : /**\brief \f$ R += \alpha \frac{\partial}{\partial T}det(T) \f$
103 : : *
104 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
105 : : * symmetric matrix.
106 : : */
107 : : inline void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T );
108 : :
109 : : /**\brief \f$ R += \alpha \frac{\partial}{\partial T}det(T) \f$
110 : : *
111 : : *\param R The 3 blocks of the upper triangular portion of a 4x4
112 : : * symmetric matrix.
113 : : */
114 : : inline void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha );
115 : 0 : inline void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& )
116 : : {
117 : 0 : pluseq_scaled_2nd_deriv_of_det( R, alpha );
118 : 0 : }
119 : :
120 : : /**\brief \f$ R = \alpha \frac{\partial}{\partial T}det(T) \f$
121 : : *
122 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
123 : : * symmetric matrix.
124 : : */
125 : : inline void set_scaled_2nd_deriv_of_det( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T );
126 : :
127 : : /**\brief \f$ R = \alpha \frac{\partial}{\partial T}det(T) \f$
128 : : *
129 : : *\param R The 3 blocks of the upper triangular portion of a 4x4
130 : : * symmetric matrix.
131 : : */
132 : : inline void set_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha );
133 : : inline void set_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& )
134 : : {
135 : : set_scaled_2nd_deriv_of_det( R, alpha );
136 : : }
137 : :
138 : : /**\brief \f$ R += \alpha \frac{\partial^2}{\partial T^2}tr(adj T) \f$
139 : : *
140 : : * Note: 2nd deriv is a constant, so T is not passed as an argument.
141 : : */
142 : : inline void pluseq_scaled_2nd_deriv_tr_adj( MsqMatrix< 2, 2 > R[3], double alpha );
143 : :
144 : : /**\brief \f$ R += \alpha \frac{\partial^2}{\partial T^2}tr(adj T) \f$
145 : : *
146 : : * Note: 2nd deriv is a constant, so T is not passed as an argument.
147 : : */
148 : : inline void pluseq_scaled_2nd_deriv_tr_adj( MsqMatrix< 3, 3 > R[6], double alpha );
149 : :
150 : : /**\brief \f$ R += \alpha \frac{\partial^2}{\partial T^2}|adj T|^2 \f$
151 : : */
152 : : inline void set_scaled_2nd_deriv_norm_sqr_adj( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& T );
153 : :
154 : : /**\brief \f$ R += \alpha \frac{\partial^2}{\partial T^2}|adj T|^2 \f$
155 : : */
156 : : inline void set_scaled_2nd_deriv_norm_sqr_adj( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T );
157 : :
158 : : /**\brief \f$ R += \alpha \left( M \otimes M \right) \f$
159 : : *
160 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
161 : : * symmetric matrix.
162 : : */
163 : : template < unsigned D >
164 : : inline void pluseq_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
165 : : const MsqMatrix< D, D >& M );
166 : 118180 : inline void pluseq_scaled_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& M )
167 : : {
168 : 118180 : pluseq_scaled_outer_product_t< 3 >( R, alpha, M );
169 : 118180 : }
170 : :
171 : 0 : inline void pluseq_scaled_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& M )
172 : : {
173 : 0 : pluseq_scaled_outer_product_t< 2 >( R, alpha, M );
174 : 0 : }
175 : :
176 : : /**\brief \f$ R = \alpha \left( M \otimes M \right) \f$
177 : : *
178 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
179 : : * symmetric matrix.
180 : : */
181 : : template < unsigned D >
182 : : inline void set_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
183 : : const MsqMatrix< D, D >& M );
184 : :
185 : 236360 : inline void set_scaled_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& M )
186 : : {
187 : 236360 : set_scaled_outer_product_t< 3 >( R, alpha, M );
188 : 236360 : }
189 : :
190 : 0 : inline void set_scaled_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& M )
191 : : {
192 : 0 : set_scaled_outer_product_t< 2 >( R, alpha, M );
193 : 0 : }
194 : :
195 : : /**\brief \f$ R += \alpha \left( A \otimes B + B \otimes A \right) \f$
196 : : *
197 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
198 : : * symmetric matrix.
199 : : */
200 : : inline void pluseq_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A,
201 : : const MsqMatrix< 3, 3 >& B );
202 : :
203 : : /**\brief \f$ R += \alpha \left( A \otimes B + B \otimes A \right) \f$
204 : : *
205 : : *\param R The 3 blocks of the upper triangular portion of a 4x4
206 : : * symmetric matrix.
207 : : */
208 : : inline void pluseq_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A,
209 : : const MsqMatrix< 2, 2 >& B );
210 : :
211 : : /**\brief \f$ R = \alpha \left( A \otimes B + B \otimes A \right) \f$
212 : : *
213 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
214 : : * symmetric matrix.
215 : : */
216 : : inline void set_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A,
217 : : const MsqMatrix< 3, 3 >& B );
218 : :
219 : : /**\brief \f$ R = \alpha \left( A \otimes B + B \otimes A \right) \f$
220 : : *
221 : : *\param R The 3 blocks of the upper triangular portion of a 4x4
222 : : * symmetric matrix.
223 : : */
224 : : inline void set_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A,
225 : : const MsqMatrix< 2, 2 >& B );
226 : :
227 : : /**\brief \f$ R += \alpha (I \otimes I) \f$
228 : : *
229 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
230 : : * symmetric matrix.
231 : : */
232 : : inline void pluseq_scaled_outer_product_I_I( MsqMatrix< 3, 3 > R[6], double alpha );
233 : :
234 : : /**\brief \f$ R += \alpha (I \otimes I) \f$
235 : : *
236 : : *\param R The 3 blocks of the upper triangular portion of a 4x4
237 : : * symmetric matrix.
238 : : */
239 : : inline void pluseq_scaled_outer_product_I_I( MsqMatrix< 2, 2 > R[3], double alpha );
240 : :
241 : : /**\brief \f$ R += I \otimes A \f$
242 : : *
243 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
244 : : * symmetric matrix.
245 : : */
246 : : inline void pluseq_I_outer_product( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A );
247 : :
248 : : /**\brief \f$ R += A \otimes I \f$
249 : : *
250 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
251 : : * symmetric matrix.
252 : : */
253 : : inline void pluseq_outer_product_I( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A );
254 : :
255 : : /**\brief \f$ R += \alpha \left( I \otimes A + A \otimes I \right) \f$
256 : : *
257 : : *\param R The 6 blocks of the upper triangular portion of a 9x9
258 : : * symmetric matrix.
259 : : */
260 : : template < unsigned D >
261 : : inline void pluseq_scaled_sum_outer_product_I_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
262 : : const MsqMatrix< D, D >& A );
263 : 0 : inline void pluseq_scaled_sum_outer_product_I( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A_in )
264 : : {
265 : 0 : pluseq_scaled_sum_outer_product_I_t< 3 >( R, alpha, A_in );
266 : 0 : }
267 : :
268 : 0 : inline void pluseq_scaled_sum_outer_product_I( MsqMatrix< 2, 2 > R[2], double alpha, const MsqMatrix< 2, 2 >& A_in )
269 : : {
270 : 0 : pluseq_scaled_sum_outer_product_I_t< 2 >( R, alpha, A_in );
271 : 0 : }
272 : :
273 : : /**\brief \f$ \frac{\partial^2 f}{\partial (AZ)^2} \Rightarrow \frac{\partial^2 f}{\partial A^2} \f$
274 : : *
275 : : * Given the second derivatives of a function with respect to
276 : : * a matrix product, calculate the derivatives of the function
277 : : * with respect to the first matrix in the product.
278 : : */
279 : : template < unsigned D >
280 : : inline void second_deriv_wrt_product_factor_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], const MsqMatrix< D, D >& Z );
281 : :
282 : 236360 : inline void second_deriv_wrt_product_factor( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& Z )
283 : : {
284 : 236360 : second_deriv_wrt_product_factor_t< 3 >( R, Z );
285 : 236360 : }
286 : :
287 : : inline void second_deriv_wrt_product_factor( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& Z )
288 : : {
289 : : second_deriv_wrt_product_factor_t< 2 >( R, Z );
290 : : }
291 : :
292 : : /**\brief \f$ R = \alpha * \frac{\partial^2 \psi(T)}{\partial T^2} \f$
293 : : *
294 : : * \f$ \psi(T) = \sqrt{ |T|^2 + 2 \tau } \f$
295 : : */
296 : : inline void set_scaled_2nd_deriv_wrt_psi( MsqMatrix< 2, 2 > R[3], const double alpha, const double psi,
297 : : const MsqMatrix< 2, 2 >& T );
298 : :
299 : : /**\brief \f$ R = R + \alpha * Z \f$
300 : : */
301 : : template < unsigned D >
302 : : inline void pluseq_scaled_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
303 : : const MsqMatrix< D, D > Z[D * ( D + 1 ) / 2] );
304 : :
305 : 0 : void set_scaled_I( MsqMatrix< 3, 3 > R[6], double alpha )
306 : : {
307 : 0 : R[0] = R[3] = R[5] = MsqMatrix< 3, 3 >( alpha );
308 : 0 : R[1] = R[2] = R[4] = MsqMatrix< 3, 3 >( 0.0 );
309 : 0 : }
310 : :
311 : 1046099 : void pluseq_scaled_I( MsqMatrix< 3, 3 >& R, double alpha )
312 : : {
313 : 1046099 : R( 0, 0 ) += alpha;
314 : 1046099 : R( 1, 1 ) += alpha;
315 : 1046099 : R( 2, 2 ) += alpha;
316 : 1046099 : }
317 : :
318 : 562262 : void pluseq_scaled_I( MsqMatrix< 2, 2 >& R, double alpha )
319 : : {
320 : 562262 : R( 0, 0 ) += alpha;
321 : 562262 : R( 1, 1 ) += alpha;
322 : 562262 : }
323 : :
324 : 236360 : void pluseq_scaled_I( MsqMatrix< 3, 3 > R[6], double alpha )
325 : : {
326 : 236360 : pluseq_scaled_I( R[0], alpha );
327 : 236360 : pluseq_scaled_I( R[3], alpha );
328 : 236360 : pluseq_scaled_I( R[5], alpha );
329 : 236360 : }
330 : :
331 : 0 : void set_scaled_I( MsqMatrix< 2, 2 > R[3], double alpha )
332 : : {
333 : 0 : R[0] = R[2] = MsqMatrix< 2, 2 >( alpha );
334 : 0 : R[1] = MsqMatrix< 2, 2 >( 0.0 );
335 : 0 : }
336 : :
337 : 0 : void pluseq_scaled_I( MsqMatrix< 2, 2 > R[3], double alpha )
338 : : {
339 : 0 : pluseq_scaled_I( R[0], alpha );
340 : 0 : pluseq_scaled_I( R[2], alpha );
341 : 0 : }
342 : :
343 : 236360 : void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T )
344 : : {
345 : 236360 : MsqMatrix< 3, 3 > A( T );
346 [ + - ]: 236360 : A *= alpha;
347 : :
348 [ + - ][ + - ]: 236360 : R[1]( 0, 1 ) += A( 2, 2 );
349 [ + - ][ + - ]: 236360 : R[1]( 1, 0 ) -= A( 2, 2 );
350 [ + - ][ + - ]: 236360 : R[1]( 0, 2 ) -= A( 2, 1 );
351 [ + - ][ + - ]: 236360 : R[1]( 2, 0 ) += A( 2, 1 );
352 [ + - ][ + - ]: 236360 : R[1]( 1, 2 ) += A( 2, 0 );
353 [ + - ][ + - ]: 236360 : R[1]( 2, 1 ) -= A( 2, 0 );
354 : :
355 [ + - ][ + - ]: 236360 : R[2]( 0, 1 ) -= A( 1, 2 );
356 [ + - ][ + - ]: 236360 : R[2]( 1, 0 ) += A( 1, 2 );
357 [ + - ][ + - ]: 236360 : R[2]( 0, 2 ) += A( 1, 1 );
358 [ + - ][ + - ]: 236360 : R[2]( 2, 0 ) -= A( 1, 1 );
359 [ + - ][ + - ]: 236360 : R[2]( 1, 2 ) -= A( 1, 0 );
360 [ + - ][ + - ]: 236360 : R[2]( 2, 1 ) += A( 1, 0 );
361 : :
362 [ + - ][ + - ]: 236360 : R[4]( 0, 1 ) += A( 0, 2 );
363 [ + - ][ + - ]: 236360 : R[4]( 1, 0 ) -= A( 0, 2 );
364 [ + - ][ + - ]: 236360 : R[4]( 0, 2 ) -= A( 0, 1 );
365 [ + - ][ + - ]: 236360 : R[4]( 2, 0 ) += A( 0, 1 );
366 [ + - ][ + - ]: 236360 : R[4]( 1, 2 ) += A( 0, 0 );
367 [ + - ][ + - ]: 236360 : R[4]( 2, 1 ) -= A( 0, 0 );
368 : 236360 : }
369 : :
370 : 0 : void set_scaled_2nd_deriv_of_det( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T )
371 : : {
372 : 0 : MsqMatrix< 3, 3 > A( T );
373 [ # # ]: 0 : A *= alpha;
374 : :
375 [ # # ]: 0 : R[0] = R[3] = R[5] = MsqMatrix< 3, 3 >( 0.0 );
376 : :
377 [ # # ][ # # ]: 0 : R[1]( 0, 0 ) = R[1]( 1, 1 ) = R[1]( 2, 2 ) = 0;
[ # # ]
378 [ # # ][ # # ]: 0 : R[1]( 0, 1 ) = A( 2, 2 );
379 [ # # ][ # # ]: 0 : R[1]( 1, 0 ) = -A( 2, 2 );
380 [ # # ][ # # ]: 0 : R[1]( 0, 2 ) = -A( 2, 1 );
381 [ # # ][ # # ]: 0 : R[1]( 2, 0 ) = A( 2, 1 );
382 [ # # ][ # # ]: 0 : R[1]( 1, 2 ) = A( 2, 0 );
383 [ # # ][ # # ]: 0 : R[1]( 2, 1 ) = -A( 2, 0 );
384 : :
385 [ # # ][ # # ]: 0 : R[2]( 0, 0 ) = R[2]( 1, 1 ) = R[2]( 2, 2 ) = 0;
[ # # ]
386 [ # # ][ # # ]: 0 : R[2]( 0, 1 ) = -A( 1, 2 );
387 [ # # ][ # # ]: 0 : R[2]( 1, 0 ) = A( 1, 2 );
388 [ # # ][ # # ]: 0 : R[2]( 0, 2 ) = A( 1, 1 );
389 [ # # ][ # # ]: 0 : R[2]( 2, 0 ) = -A( 1, 1 );
390 [ # # ][ # # ]: 0 : R[2]( 1, 2 ) = -A( 1, 0 );
391 [ # # ][ # # ]: 0 : R[2]( 2, 1 ) = A( 1, 0 );
392 : :
393 [ # # ][ # # ]: 0 : R[4]( 0, 0 ) = R[4]( 1, 1 ) = R[4]( 2, 2 ) = 0;
[ # # ]
394 [ # # ][ # # ]: 0 : R[4]( 0, 1 ) = A( 0, 2 );
395 [ # # ][ # # ]: 0 : R[4]( 1, 0 ) = -A( 0, 2 );
396 [ # # ][ # # ]: 0 : R[4]( 0, 2 ) = -A( 0, 1 );
397 [ # # ][ # # ]: 0 : R[4]( 2, 0 ) = A( 0, 1 );
398 [ # # ][ # # ]: 0 : R[4]( 1, 2 ) = A( 0, 0 );
399 [ # # ][ # # ]: 0 : R[4]( 2, 1 ) = -A( 0, 0 );
400 : 0 : }
401 : :
402 : 0 : void pluseq_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha )
403 : : {
404 : 0 : R[1]( 0, 1 ) += alpha;
405 : 0 : R[1]( 1, 0 ) -= alpha;
406 : 0 : }
407 : :
408 : : void set_scaled_2nd_deriv_of_det( MsqMatrix< 2, 2 > R[3], double alpha )
409 : : {
410 : : R[0] = R[2] = MsqMatrix< 2, 2 >( 0.0 );
411 : : R[1]( 0, 0 ) = 0.0;
412 : : R[1]( 0, 1 ) = alpha;
413 : : R[1]( 1, 0 ) = -alpha;
414 : : R[1]( 1, 1 ) = 0.0;
415 : : }
416 : :
417 : : void pluseq_scaled_2nd_deriv_tr_adj( MsqMatrix< 2, 2 >*, double )
418 : : {
419 : : // 2nd derivative is zero
420 : : }
421 : :
422 : 0 : void pluseq_scaled_2nd_deriv_tr_adj( MsqMatrix< 3, 3 > R[6], double alpha )
423 : : {
424 : 0 : R[1]( 0, 1 ) += alpha;
425 : 0 : R[1]( 1, 0 ) -= alpha;
426 : 0 : R[2]( 0, 2 ) += alpha;
427 : 0 : R[2]( 2, 0 ) -= alpha;
428 : 0 : R[4]( 1, 2 ) += alpha;
429 : 0 : R[4]( 2, 1 ) -= alpha;
430 : 0 : }
431 : :
432 : : void set_scaled_2nd_deriv_norm_sqr_adj( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& )
433 : : {
434 : : set_scaled_I( R, 2 * alpha );
435 : : }
436 : :
437 : 0 : void set_scaled_2nd_deriv_norm_sqr_adj( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& T )
438 : : {
439 [ # # ]: 0 : set_scaled_outer_product( R, 1, T );
440 : : double tmp01, tmp02, tmp12;
441 : : // R[1] = 2*R[1] - transpose(R[1]);
442 [ # # ]: 0 : tmp01 = R[1]( 0, 1 );
443 [ # # ]: 0 : tmp02 = R[1]( 0, 2 );
444 [ # # ]: 0 : tmp12 = R[1]( 1, 2 );
445 [ # # ][ # # ]: 0 : R[1]( 0, 1 ) = 2 * R[1]( 0, 1 ) - R[1]( 1, 0 );
[ # # ]
446 [ # # ][ # # ]: 0 : R[1]( 0, 2 ) = 2 * R[1]( 0, 2 ) - R[1]( 2, 0 );
[ # # ]
447 [ # # ][ # # ]: 0 : R[1]( 1, 2 ) = 2 * R[1]( 1, 2 ) - R[1]( 2, 1 );
[ # # ]
448 [ # # ][ # # ]: 0 : R[1]( 1, 0 ) = 2 * R[1]( 1, 0 ) - tmp01;
449 [ # # ][ # # ]: 0 : R[1]( 2, 0 ) = 2 * R[1]( 2, 0 ) - tmp02;
450 [ # # ][ # # ]: 0 : R[1]( 2, 1 ) = 2 * R[1]( 2, 1 ) - tmp12;
451 : : // R[2] = 2*R[2] - transpose(R[1]);
452 [ # # ]: 0 : tmp01 = R[2]( 0, 1 );
453 [ # # ]: 0 : tmp02 = R[2]( 0, 2 );
454 [ # # ]: 0 : tmp12 = R[2]( 1, 2 );
455 [ # # ][ # # ]: 0 : R[2]( 0, 1 ) = 2 * R[2]( 0, 1 ) - R[2]( 1, 0 );
[ # # ]
456 [ # # ][ # # ]: 0 : R[2]( 0, 2 ) = 2 * R[2]( 0, 2 ) - R[2]( 2, 0 );
[ # # ]
457 [ # # ][ # # ]: 0 : R[2]( 1, 2 ) = 2 * R[2]( 1, 2 ) - R[2]( 2, 1 );
[ # # ]
458 [ # # ][ # # ]: 0 : R[2]( 1, 0 ) = 2 * R[2]( 1, 0 ) - tmp01;
459 [ # # ][ # # ]: 0 : R[2]( 2, 0 ) = 2 * R[2]( 2, 0 ) - tmp02;
460 [ # # ][ # # ]: 0 : R[2]( 2, 1 ) = 2 * R[2]( 2, 1 ) - tmp12;
461 : : // R[4] = 2*R[4] - transpose(R[1]);
462 [ # # ]: 0 : tmp01 = R[4]( 0, 1 );
463 [ # # ]: 0 : tmp02 = R[4]( 0, 2 );
464 [ # # ]: 0 : tmp12 = R[4]( 1, 2 );
465 [ # # ][ # # ]: 0 : R[4]( 0, 1 ) = 2 * R[4]( 0, 1 ) - R[4]( 1, 0 );
[ # # ]
466 [ # # ][ # # ]: 0 : R[4]( 0, 2 ) = 2 * R[4]( 0, 2 ) - R[4]( 2, 0 );
[ # # ]
467 [ # # ][ # # ]: 0 : R[4]( 1, 2 ) = 2 * R[4]( 1, 2 ) - R[4]( 2, 1 );
[ # # ]
468 [ # # ][ # # ]: 0 : R[4]( 1, 0 ) = 2 * R[4]( 1, 0 ) - tmp01;
469 [ # # ][ # # ]: 0 : R[4]( 2, 0 ) = 2 * R[4]( 2, 0 ) - tmp02;
470 [ # # ][ # # ]: 0 : R[4]( 2, 1 ) = 2 * R[4]( 2, 1 ) - tmp12;
471 : :
472 [ # # ][ # # ]: 0 : const MsqMatrix< 3, 3 > TpT = transpose( T ) * T;
473 [ # # ]: 0 : R[0] -= TpT;
474 [ # # ]: 0 : R[3] -= TpT;
475 [ # # ]: 0 : R[5] -= TpT;
476 : :
477 [ # # ][ # # ]: 0 : const MsqMatrix< 3, 3 > TTp = T * transpose( T );
478 [ # # ]: 0 : const double ns = sqr_Frobenius( T );
479 [ # # ][ # # ]: 0 : R[0]( 0, 0 ) += ns - TTp( 0, 0 );
480 [ # # ][ # # ]: 0 : R[0]( 1, 1 ) += ns - TTp( 0, 0 );
481 [ # # ][ # # ]: 0 : R[0]( 2, 2 ) += ns - TTp( 0, 0 );
482 [ # # ][ # # ]: 0 : R[1]( 0, 0 ) += -TTp( 0, 1 );
483 [ # # ][ # # ]: 0 : R[1]( 1, 1 ) += -TTp( 0, 1 );
484 [ # # ][ # # ]: 0 : R[1]( 2, 2 ) += -TTp( 0, 1 );
485 [ # # ][ # # ]: 0 : R[2]( 0, 0 ) += -TTp( 0, 2 );
486 [ # # ][ # # ]: 0 : R[2]( 1, 1 ) += -TTp( 0, 2 );
487 [ # # ][ # # ]: 0 : R[2]( 2, 2 ) += -TTp( 0, 2 );
488 [ # # ][ # # ]: 0 : R[3]( 0, 0 ) += ns - TTp( 1, 1 );
489 [ # # ][ # # ]: 0 : R[3]( 1, 1 ) += ns - TTp( 1, 1 );
490 [ # # ][ # # ]: 0 : R[3]( 2, 2 ) += ns - TTp( 1, 1 );
491 [ # # ][ # # ]: 0 : R[4]( 0, 0 ) += -TTp( 1, 2 );
492 [ # # ][ # # ]: 0 : R[4]( 1, 1 ) += -TTp( 1, 2 );
493 [ # # ][ # # ]: 0 : R[4]( 2, 2 ) += -TTp( 1, 2 );
494 [ # # ][ # # ]: 0 : R[5]( 0, 0 ) += ns - TTp( 2, 2 );
495 [ # # ][ # # ]: 0 : R[5]( 1, 1 ) += ns - TTp( 2, 2 );
496 [ # # ][ # # ]: 0 : R[5]( 2, 2 ) += ns - TTp( 2, 2 );
497 : :
498 : 0 : alpha *= 2;
499 [ # # ]: 0 : R[0] *= alpha;
500 [ # # ]: 0 : R[1] *= alpha;
501 [ # # ]: 0 : R[2] *= alpha;
502 [ # # ]: 0 : R[3] *= alpha;
503 [ # # ]: 0 : R[4] *= alpha;
504 [ # # ]: 0 : R[5] *= alpha;
505 : 0 : }
506 : :
507 : : #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
508 : : template < unsigned D >
509 : 118180 : void pluseq_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
510 : : {
511 : 118180 : MsqMatrix< D, D > aM( M );
512 [ + - ][ # # ]: 118180 : aM *= alpha;
513 : 118180 : unsigned h = 0;
514 [ + + ][ # # ]: 472720 : for( unsigned i = 0; i < D; ++i )
515 [ + + ][ # # ]: 1063620 : for( unsigned j = i; j < D; ++j )
516 [ + - ][ + - ]: 709080 : R[h++] += transpose( aM.row( i ) ) * M.row( j );
[ + - ][ + - ]
[ + - ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
517 : 118180 : }
518 : : #else
519 : : template < unsigned D >
520 : : void pluseq_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
521 : : {
522 : : MsqMatrix< D, D > aM( transpose( M ) );
523 : : aM *= alpha;
524 : : unsigned h = 0;
525 : : for( unsigned i = 0; i < D; ++i )
526 : : for( unsigned j = i; j < D; ++j )
527 : : R[h++] += M.column( i ) * aM.row( 0 );
528 : : }
529 : : #endif
530 : :
531 : : #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
532 : : template < unsigned D >
533 : 236360 : void set_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
534 : : {
535 : 236360 : MsqMatrix< D, D > aM( M );
536 [ # # ][ + - ]: 236360 : aM *= alpha;
537 : 236360 : unsigned h = 0;
538 [ # # ][ + + ]: 945440 : for( unsigned i = 0; i < D; ++i )
539 [ # # ][ + + ]: 2127240 : for( unsigned j = i; j < D; ++j )
540 [ # # ][ # # ]: 1418160 : R[h++] = transpose( aM.row( i ) ) * M.row( j );
[ # # ][ # # ]
[ + - ][ + - ]
[ + - ][ + - ]
541 : 236360 : }
542 : : #else
543 : : template < unsigned D >
544 : : void set_scaled_outer_product_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha, const MsqMatrix< D, D >& M )
545 : : {
546 : : MsqMatrix< D, D > aM( transpose( M ) );
547 : : aM *= alpha;
548 : : unsigned h = 0;
549 : : for( unsigned i = 0; i < D; ++i )
550 : : for( unsigned j = i; j < D; ++j )
551 : : R[h++] = M.column( i ) * aM.row( j );
552 : : }
553 : : #endif
554 : :
555 : : #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
556 : 236360 : void pluseq_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A_in,
557 : : const MsqMatrix< 3, 3 >& B )
558 : : {
559 : : // apply scalar first
560 [ + - ]: 236360 : MsqMatrix< 3, 3 > A( A_in ), tmp;
561 [ + - ]: 236360 : A *= alpha;
562 : :
563 : : // block 0,0
564 [ + - ][ + - ]: 236360 : tmp = transpose( A.row( 0 ) ) * B.row( 0 );
[ + - ][ + - ]
565 [ + - ]: 236360 : R[0] += tmp;
566 [ + - ][ + - ]: 236360 : R[0] += transpose( tmp );
567 : :
568 : : // block 1,1
569 [ + - ][ + - ]: 236360 : tmp = transpose( A.row( 1 ) ) * B.row( 1 );
[ + - ][ + - ]
570 [ + - ]: 236360 : R[3] += tmp;
571 [ + - ][ + - ]: 236360 : R[3] += transpose( tmp );
572 : :
573 : : // block 2,2
574 [ + - ][ + - ]: 236360 : tmp = transpose( A.row( 2 ) ) * B.row( 2 );
[ + - ][ + - ]
575 [ + - ]: 236360 : R[5] += tmp;
576 [ + - ][ + - ]: 236360 : R[5] += transpose( tmp );
577 : :
578 : : // block 0,1
579 [ + - ][ + - ]: 236360 : R[1] += transpose( A.row( 0 ) ) * B.row( 1 ) + transpose( B.row( 0 ) ) * A.row( 1 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
580 : :
581 : : // block 0,2
582 [ + - ][ + - ]: 236360 : R[2] += transpose( A.row( 0 ) ) * B.row( 2 ) + transpose( B.row( 0 ) ) * A.row( 2 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
583 : :
584 : : // block 1,2
585 [ + - ][ + - ]: 236360 : R[4] += transpose( A.row( 1 ) ) * B.row( 2 ) + transpose( B.row( 1 ) ) * A.row( 2 );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
586 : 236360 : }
587 : : void set_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A_in,
588 : : const MsqMatrix< 3, 3 >& B )
589 : : {
590 : : // apply scalar first
591 : : MsqMatrix< 3, 3 > A( A_in );
592 : : A *= alpha;
593 : :
594 : : // block 0,0
595 : : R[0] = transpose( A.row( 0 ) ) * B.row( 0 );
596 : : R[0] += transpose( R[0] );
597 : :
598 : : // block 1,1
599 : : R[3] = transpose( A.row( 1 ) ) * B.row( 1 );
600 : : R[3] += transpose( R[3] );
601 : :
602 : : // block 2,2
603 : : R[5] = transpose( A.row( 2 ) ) * B.row( 2 );
604 : : R[5] += transpose( R[5] );
605 : :
606 : : // block 0,1
607 : : R[1] = transpose( A.row( 0 ) ) * B.row( 1 );
608 : : R[1] += transpose( B.row( 0 ) ) * A.row( 1 );
609 : :
610 : : // block 0,2
611 : : R[2] = transpose( A.row( 0 ) ) * B.row( 2 );
612 : : R[2] += transpose( B.row( 0 ) ) * A.row( 2 );
613 : :
614 : : // block 1,2
615 : : R[4] = transpose( A.row( 1 ) ) * B.row( 2 );
616 : : R[4] += transpose( B.row( 1 ) ) * A.row( 2 );
617 : : }
618 : : #else
619 : : void pluseq_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A_in,
620 : : const MsqMatrix< 3, 3 >& B )
621 : : {
622 : : // apply scalar first
623 : : MsqMatrix< 3, 3 > A( A_in ), tmp;
624 : : A *= alpha;
625 : :
626 : : // block 0,0
627 : : tmp = A.column( 0 ) * transpose( B.column( 0 ) );
628 : : R[0] += tmp;
629 : : R[0] += transpose( tmp );
630 : :
631 : : // block 1,1
632 : : tmp = A.column( 1 ) * transpose( B.column( 1 ) );
633 : : R[3] += tmp;
634 : : R[3] += transpose( tmp );
635 : :
636 : : // block 2,2
637 : : tmp = A.column( 2 ) * transpose( B.column( 2 ) );
638 : : R[5] += tmp;
639 : : R[5] += transpose( tmp );
640 : :
641 : : // block 0,1
642 : : R[1] += A.column( 0 ) * transpose( B.column( 1 ) ) + B.column( 0 ) * transpose( A.column( 1 ) );
643 : :
644 : : // block 0,2
645 : : R[2] += A.column( 0 ) * transpose( B.column( 2 ) ) + B.column( 0 ) * transpose( A.column( 2 ) );
646 : :
647 : : // block 1,2
648 : : R[4] += A.column( 1 ) * transpose( B.column( 2 ) ) + B.column( 1 ) * transpose( A.column( 2 ) );
649 : : }
650 : : void set_scaled_sum_outer_product( MsqMatrix< 3, 3 > R[6], double alpha, const MsqMatrix< 3, 3 >& A_in,
651 : : const MsqMatrix< 3, 3 >& B )
652 : : {
653 : : // apply scalar first
654 : : MsqMatrix< 3, 3 > A( A_in );
655 : : A *= alpha;
656 : :
657 : : // block 0,0
658 : : R[0] = A.column( 0 ) * transpose( B.column( 0 ) );
659 : : R[0] += transpose( R[0] );
660 : :
661 : : // block 1,1
662 : : R[3] = A.column( 1 ) * transpose( B.column( 1 ) );
663 : : R[3] += transpose( R[3] );
664 : :
665 : : // block 2,2
666 : : R[5] = A.column( 2 ) * transpose( B.column( 2 ) );
667 : : R[5] += transpose( R[5] );
668 : :
669 : : // block 0,1
670 : : R[1] = A.column( 0 ) * transpose( B.column( 1 ) );
671 : : R[1] += B.column( 0 ) * transpose( A.column( 1 ) );
672 : :
673 : : // block 0,2
674 : : R[2] = A.column( 0 ) * transpose( B.column( 2 ) );
675 : : R[2] += B.column( 0 ) * transpose( A.column( 2 ) );
676 : :
677 : : // block 1,2
678 : : R[4] = A.column( 1 ) * transpose( B.column( 2 ) );
679 : : R[4] += B.column( 1 ) * transpose( A.column( 2 ) );
680 : : }
681 : : #endif
682 : :
683 : : #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
684 : 0 : void pluseq_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A_in,
685 : : const MsqMatrix< 2, 2 >& B )
686 : : {
687 : : // apply scalar first
688 [ # # ]: 0 : MsqMatrix< 2, 2 > A( A_in ), tmp;
689 [ # # ]: 0 : A *= alpha;
690 : :
691 : : // block 0,0
692 [ # # ][ # # ]: 0 : tmp = transpose( A.row( 0 ) ) * B.row( 0 );
[ # # ][ # # ]
693 [ # # ]: 0 : R[0] += tmp;
694 [ # # ][ # # ]: 0 : R[0] += transpose( tmp );
695 : :
696 : : // block 1,1
697 [ # # ][ # # ]: 0 : tmp = transpose( A.row( 1 ) ) * B.row( 1 );
[ # # ][ # # ]
698 [ # # ]: 0 : R[2] += tmp;
699 [ # # ][ # # ]: 0 : R[2] += transpose( tmp );
700 : :
701 : : // block 0,1
702 [ # # ][ # # ]: 0 : R[1] += transpose( A.row( 0 ) ) * B.row( 1 ) + transpose( B.row( 0 ) ) * A.row( 1 );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
703 : 0 : }
704 : 0 : void set_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A_in,
705 : : const MsqMatrix< 2, 2 >& B )
706 : : {
707 : : // apply scalar first
708 : 0 : MsqMatrix< 2, 2 > A( A_in );
709 [ # # ]: 0 : A *= alpha;
710 : :
711 : : // block 0,0
712 [ # # ][ # # ]: 0 : R[0] = transpose( A.row( 0 ) ) * B.row( 0 );
[ # # ][ # # ]
713 [ # # ][ # # ]: 0 : R[0] += transpose( R[0] );
714 : :
715 : : // block 1,1
716 [ # # ][ # # ]: 0 : R[2] = transpose( A.row( 1 ) ) * B.row( 1 );
[ # # ][ # # ]
717 [ # # ][ # # ]: 0 : R[2] += transpose( R[2] );
718 : :
719 : : // block 0,1
720 [ # # ][ # # ]: 0 : R[1] = transpose( A.row( 0 ) ) * B.row( 1 );
[ # # ][ # # ]
721 [ # # ][ # # ]: 0 : R[1] += transpose( B.row( 0 ) ) * A.row( 1 );
[ # # ][ # # ]
[ # # ]
722 : 0 : }
723 : : #else
724 : : void pluseq_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A_in,
725 : : const MsqMatrix< 2, 2 >& B )
726 : : {
727 : : // apply scalar first
728 : : MsqMatrix< 2, 2 > A( A_in ), tmp;
729 : : A *= alpha;
730 : :
731 : : // block 0,0
732 : : tmp = A.column( 0 ) * transpose( B.column( 0 ) );
733 : : R[0] += tmp;
734 : : R[0] += transpose( tmp );
735 : :
736 : : // block 1,1
737 : : tmp = A.column( 1 ) * transpose( B.column( 1 ) );
738 : : R[2] += tmp;
739 : : R[2] += transpose( tmp );
740 : :
741 : : // block 0,1
742 : : R[1] += A.column( 0 ) * transpose( B.column( 1 ) ) + B.column( 0 ) * transpose( A.column( 1 ) );
743 : : }
744 : : void set_scaled_sum_outer_product( MsqMatrix< 2, 2 > R[3], double alpha, const MsqMatrix< 2, 2 >& A_in,
745 : : const MsqMatrix< 2, 2 >& B )
746 : : {
747 : : // apply scalar first
748 : : MsqMatrix< 2, 2 > A( A_in );
749 : : A *= alpha;
750 : :
751 : : // block 0,0
752 : : R[0] = A.column( 0 ) * transpose( B.column( 0 ) );
753 : : R[0] += transpose( R[0] );
754 : :
755 : : // block 1,1
756 : : R[2] = A.column( 1 ) * transpose( B.column( 1 ) );
757 : : R[2] += transpose( R[2] );
758 : :
759 : : // block 0,1
760 : : R[1] = A.column( 0 ) * transpose( B.column( 1 ) );
761 : : R[1] += B.column( 0 ) * transpose( A.column( 1 ) );
762 : : }
763 : : #endif
764 : :
765 : 0 : void pluseq_scaled_outer_product_I_I( MsqMatrix< 3, 3 > R[6], double alpha )
766 : : {
767 : 0 : R[0]( 0, 0 ) += alpha;
768 : 0 : R[1]( 0, 1 ) += alpha;
769 : 0 : R[2]( 0, 2 ) += alpha;
770 : 0 : R[3]( 1, 1 ) += alpha;
771 : 0 : R[4]( 1, 2 ) += alpha;
772 : 0 : R[5]( 2, 2 ) += alpha;
773 : 0 : }
774 : :
775 : 0 : void pluseq_scaled_outer_product_I_I( MsqMatrix< 2, 2 > R[3], double alpha )
776 : : {
777 : 0 : R[0]( 0, 0 ) += alpha;
778 : 0 : R[1]( 0, 1 ) += alpha;
779 : 0 : R[2]( 1, 1 ) += alpha;
780 : 0 : }
781 : :
782 : : #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
783 : 0 : inline void pluseq_I_outer_product( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
784 : : {
785 [ # # ]: 0 : R[0].add_row( 0, A.row( 0 ) );
786 [ # # ]: 0 : R[1].add_row( 0, A.row( 1 ) );
787 [ # # ]: 0 : R[2].add_row( 0, A.row( 2 ) );
788 [ # # ]: 0 : R[3].add_row( 1, A.row( 1 ) );
789 [ # # ]: 0 : R[4].add_row( 1, A.row( 2 ) );
790 [ # # ]: 0 : R[5].add_row( 2, A.row( 2 ) );
791 : 0 : }
792 : 0 : inline void pluseq_I_outer_product( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
793 : : {
794 [ # # ]: 0 : R[0].add_row( 0, A.row( 0 ) );
795 [ # # ]: 0 : R[1].add_row( 0, A.row( 1 ) );
796 [ # # ]: 0 : R[2].add_row( 1, A.row( 1 ) );
797 : 0 : }
798 : : #else
799 : : inline void pluseq_I_outer_product( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
800 : : {
801 : : R[0].add_row( 0, transpose( A.column( 0 ) ) );
802 : : R[1].add_row( 0, transpose( A.column( 1 ) ) );
803 : : R[2].add_row( 0, transpose( A.column( 2 ) ) );
804 : : R[3].add_row( 1, transpose( A.column( 1 ) ) );
805 : : R[4].add_row( 1, transpose( A.column( 2 ) ) );
806 : : R[5].add_row( 2, transpose( A.column( 2 ) ) );
807 : : }
808 : : inline void pluseq_I_outer_product( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
809 : : {
810 : : R[0].add_row( 0, transpose( A.column( 0 ) ) );
811 : : R[1].add_row( 0, transpose( A.column( 1 ) ) );
812 : : R[2].add_row( 1, transpose( A.column( 1 ) ) );
813 : : }
814 : : #endif
815 : :
816 : : #ifdef MSQ_ROW_BASED_OUTER_PRODUCT
817 : 0 : inline void pluseq_outer_product_I( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
818 : : {
819 [ # # ][ # # ]: 0 : R[0].add_column( 0, transpose( A.row( 0 ) ) );
820 [ # # ][ # # ]: 0 : R[1].add_column( 1, transpose( A.row( 0 ) ) );
821 [ # # ][ # # ]: 0 : R[2].add_column( 2, transpose( A.row( 0 ) ) );
822 [ # # ][ # # ]: 0 : R[3].add_column( 1, transpose( A.row( 1 ) ) );
823 [ # # ][ # # ]: 0 : R[4].add_column( 2, transpose( A.row( 1 ) ) );
824 [ # # ][ # # ]: 0 : R[5].add_column( 2, transpose( A.row( 2 ) ) );
825 : 0 : }
826 : 0 : inline void pluseq_outer_product_I( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
827 : : {
828 [ # # ][ # # ]: 0 : R[0].add_column( 0, transpose( A.row( 0 ) ) );
829 [ # # ][ # # ]: 0 : R[1].add_column( 1, transpose( A.row( 0 ) ) );
830 [ # # ][ # # ]: 0 : R[2].add_column( 1, transpose( A.row( 1 ) ) );
831 : 0 : }
832 : : #else
833 : : inline void pluseq_outer_product_I( MsqMatrix< 3, 3 > R[6], const MsqMatrix< 3, 3 >& A )
834 : : {
835 : : R[0].add_column( 0, A.column( 0 ) );
836 : : R[1].add_column( 1, A.column( 0 ) );
837 : : R[2].add_column( 2, A.column( 0 ) );
838 : : R[3].add_column( 1, A.column( 1 ) );
839 : : R[4].add_column( 2, A.column( 1 ) );
840 : : R[5].add_column( 2, A.column( 2 ) );
841 : : }
842 : : inline void pluseq_outer_product_I( MsqMatrix< 2, 2 > R[3], const MsqMatrix< 2, 2 >& A )
843 : : {
844 : : R[0].add_column( 0, A.column( 0 ) );
845 : : R[1].add_column( 1, A.column( 0 ) );
846 : : R[2].add_column( 1, A.column( 1 ) );
847 : : }
848 : : #endif
849 : :
850 : : template < unsigned D >
851 : 0 : void pluseq_scaled_sum_outer_product_I_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
852 : : const MsqMatrix< D, D >& A_in )
853 : : {
854 : : // apply scalar first
855 : 0 : MsqMatrix< D, D > A( A_in );
856 [ # # ][ # # ]: 0 : A *= alpha;
857 [ # # ][ # # ]: 0 : pluseq_I_outer_product( R, A );
858 [ # # ][ # # ]: 0 : pluseq_outer_product_I( R, A );
859 : 0 : }
860 : :
861 : : template < unsigned D >
862 : 236360 : void second_deriv_wrt_product_factor_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], const MsqMatrix< D, D >& Z )
863 : : {
864 [ + - ]: 236360 : const MsqMatrix< D, D > Zt = transpose( Z );
865 [ + + ]: 1654520 : for( unsigned i = 0; i < D * ( D + 1 ) / 2; ++i )
866 [ + - ][ + - ]: 1418160 : R[i] = ( Z * R[i] ) * Zt;
867 : 236360 : }
868 : :
869 : 0 : void set_scaled_2nd_deriv_wrt_psi( MsqMatrix< 2, 2 > R[3], const double alpha, const double psi,
870 : : const MsqMatrix< 2, 2 >& T )
871 : : {
872 : 0 : const double t = trace( T );
873 : 0 : const double f = alpha / ( psi * psi * psi );
874 : 0 : const double s = T( 0, 1 ) - T( 1, 0 );
875 : 0 : R[0]( 0, 0 ) = R[1]( 0, 1 ) = R[2]( 1, 1 ) = f * s * s;
876 : 0 : R[0]( 0, 1 ) = R[0]( 1, 0 ) = R[1]( 1, 1 ) = -f * s * t;
877 : 0 : R[1]( 0, 0 ) = R[2]( 0, 1 ) = R[2]( 1, 0 ) = f * s * t;
878 : 0 : R[0]( 1, 1 ) = R[2]( 0, 0 ) = -( R[1]( 1, 0 ) = -f * t * t );
879 : 0 : }
880 : :
881 : : template < unsigned D >
882 : : inline void pluseq_scaled( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha,
883 : : const MsqMatrix< D, D > Z[D * ( D + 1 ) / 2] )
884 : : {
885 : : for( unsigned i = 0; i < D * ( D + 1 ) / 2; ++i )
886 : : R[i] += alpha * Z[i];
887 : : }
888 : :
889 : : /**\brief \f$ R *= s */
890 : : template < unsigned D >
891 : 0 : inline void hess_scale_t( MsqMatrix< D, D > R[D * ( D + 1 ) / 2], double alpha )
892 : : {
893 [ # # ][ # # ]: 0 : for( unsigned i = 0; i < D * ( D + 1 ) / 2; ++i )
894 : 0 : R[i] *= alpha;
895 : 0 : }
896 : :
897 : : } // namespace MBMesquite
898 : :
899 : : #endif
|
