Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2007 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 : : (2008) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file TrustRegion.cpp
28 : : * \brief Port Todd Munson's trust region solver to Mesquite
29 : : * \author Jason Kraftcheck (Mesquite Port)
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "TrustRegion.hpp"
34 : : #include "MsqDebug.hpp"
35 : : #include "MsqError.hpp"
36 : : #include "PatchData.hpp"
37 : :
38 : : #define USE_FN_PC1 // Use 1st preconditioner from Todd's code
39 : : // (alternate is whatever is in MsqHessian already)
40 : : #undef DO_STEEP_DESC // Jason's apparently broken hack to fall back to
41 : : // steepest descent search direction
42 : :
43 : : namespace MBMesquite
44 : : {
45 : :
46 : : // Force std::vector to release allocated memory
47 : : template < typename T >
48 : 49 : static inline void free_vector( std::vector< T >& v )
49 : : {
50 [ + - ][ + - ]: 49 : std::vector< T > temp;
51 : 49 : temp.swap( v );
52 : 49 : }
53 : :
54 : 0 : std::string TrustRegion::get_name() const
55 : : {
56 [ # # ]: 0 : return "TrustRegion";
57 : : }
58 : :
59 : 7 : PatchSet* TrustRegion::get_patch_set()
60 : : {
61 : 7 : return PatchSetUser::get_patch_set();
62 : : }
63 : :
64 [ + - ][ + - ]: 7 : TrustRegion::TrustRegion( ObjectiveFunction* of ) : VertexMover( of ), PatchSetUser( true ), mMemento( 0 ) {}
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
65 : :
66 : 14 : TrustRegion::~TrustRegion()
67 : : {
68 [ - + ]: 7 : delete mMemento;
69 : 7 : mMemento = 0;
70 [ - + ]: 7 : }
71 : :
72 : 13 : void TrustRegion::initialize( PatchData& pd, MsqError& err )
73 : : {
74 [ - + ][ # # ]: 13 : mMemento = pd.create_vertices_memento( err );MSQ_CHKERR( err );
75 : 13 : }
76 : :
77 : 7 : void TrustRegion::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
78 : :
79 : 7 : void TrustRegion::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
80 : :
81 : 7 : void TrustRegion::cleanup()
82 : : {
83 : : // release Memento
84 [ + - ]: 7 : delete mMemento;
85 : 7 : mMemento = 0;
86 : : // release temporary array memory
87 : 7 : mHess.clear();
88 : 7 : free_vector( mGrad );
89 : 7 : free_vector( wVect );
90 : 7 : free_vector( zVect );
91 : 7 : free_vector( dVect );
92 : 7 : free_vector( pVect );
93 : 7 : free_vector( rVect );
94 : 7 : free_vector( preCond );
95 : 7 : }
96 : :
97 : 78 : static inline void negate( Vector3D* out, const Vector3D* in, size_t nn )
98 : : {
99 [ + + ]: 252350 : for( size_t i = 0; i < nn; ++i )
100 [ + - ]: 252272 : out[i] = -in[i];
101 : 78 : }
102 : :
103 : : // Do v += s * x, where v and x are arrays of length n
104 : 1043 : static inline void plus_eq_scaled( Vector3D* v, double s, const Vector3D* x, size_t n )
105 : : {
106 : 1043 : Vector3D* end = v + n;
107 [ + + ]: 3470512 : for( ; v != end; ++v, ++x )
108 [ + - ]: 3469469 : *v += s * *x;
109 : 1043 : }
110 : :
111 : : // Do v = s*v - x, where v and x are arrays of length n
112 : 487 : static inline void times_eq_minus( Vector3D* v, double s, const Vector3D* x, size_t n )
113 : : {
114 : 487 : Vector3D* end = v + n;
115 [ + + ]: 1613962 : for( ; v != end; ++v, ++x )
116 : : {
117 : 1613475 : *v *= s;
118 : 1613475 : *v -= *x;
119 : : }
120 : 487 : }
121 : :
122 : 47 : void TrustRegion::compute_preconditioner( MsqError&
123 : : #ifndef USE_FN_PC1
124 : : err
125 : : #endif
126 : : )
127 : : {
128 : : #ifndef USE_FN_PC1
129 : : mHessian.calculate_preconditioner( err );
130 : : #else
131 : : double dia;
132 : 47 : preCond.resize( mHess.size() );
133 [ + + ]: 140389 : for( size_t i = 0; i < mHess.size(); ++i )
134 : : {
135 : 140342 : const Matrix3D& m = *mHess.get_block( i, i );
136 : 140342 : dia = m[0][0] + m[1][1] + m[2][2];
137 [ - + ]: 140342 : preCond[i] = dia < DBL_EPSILON ? 1.0 : 1.0 / dia;
138 : : }
139 : : #endif
140 : 47 : }
141 : :
142 : 565 : void TrustRegion::apply_preconditioner( Vector3D* z, Vector3D* r,
143 : : MsqError&
144 : : #ifndef USE_FN_PC1
145 : : err
146 : : #endif
147 : : )
148 : : {
149 : : #ifndef USE_FN_PC1
150 : : mHessian.apply_preconditioner( z, r, err );
151 : : #else
152 [ + + ]: 1866312 : for( size_t i = 0; i < preCond.size(); ++i )
153 [ + - ]: 1865747 : z[i] = preCond[i] * r[i];
154 : : #endif
155 : 565 : }
156 : :
157 : 7 : void TrustRegion::optimize_vertex_positions( PatchData& pd, MsqError& err )
158 : : {
159 [ + - ]: 7 : TerminationCriterion& term = *get_inner_termination_criterion();
160 [ + - ]: 7 : OFEvaluator& func = get_objective_function_evaluator();
161 : :
162 : 7 : const double cg_tol = 1e-2;
163 : 7 : const double eta_1 = 0.01;
164 : 7 : const double eta_2 = 0.90;
165 : 7 : const double tr_incr = 10;
166 : 7 : const double tr_decr_def = 0.25;
167 : 7 : const double tr_decr_undef = 0.25;
168 : 7 : const double tr_num_tol = 1e-6;
169 : 7 : const int max_cg_iter = 10000;
170 : :
171 : 7 : double radius = 1000; /* delta*delta */
172 : :
173 [ + - ]: 7 : const int nn = pd.num_free_vertices();
174 [ + - ]: 7 : wVect.resize( nn );
175 [ + - ]: 7 : Vector3D* w = arrptr( wVect );
176 [ + - ]: 7 : zVect.resize( nn );
177 [ + - ]: 7 : Vector3D* z = arrptr( zVect );
178 [ + - ]: 7 : dVect.resize( nn );
179 [ + - ]: 7 : Vector3D* d = arrptr( dVect );
180 [ + - ]: 7 : pVect.resize( nn );
181 [ + - ]: 7 : Vector3D* p = arrptr( pVect );
182 [ + - ]: 7 : rVect.resize( nn );
183 [ + - ]: 7 : Vector3D* r = arrptr( rVect );
184 : :
185 : : double norm_r, norm_g;
186 : : double alpha, beta, kappa;
187 : : double rz, rzm1;
188 : : double dMp, norm_d, norm_dp1, norm_p;
189 : : double obj, objn;
190 : :
191 : : int cg_iter;
192 : : bool valid;
193 : :
194 [ + - ]: 7 : mHess.initialize( pd, err ); // hMesh(mesh);
195 [ + - ][ + - ]: 14 : valid = func.update( pd, obj, mGrad, mHess, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
196 [ - + ]: 7 : if( !valid )
197 : : {
198 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "Initial objective function is not valid", MsqError::INVALID_MESH );
199 : 0 : return;
200 : : }
201 [ + - ][ + - ]: 7 : compute_preconditioner( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
202 [ + - ][ + - ]: 7 : pd.recreate_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
203 : :
204 [ + - ][ + + ]: 85 : while( !term.terminate() && ( radius > 1e-20 ) )
[ + - ][ + + ]
205 : : {
206 : :
207 [ + - ][ + - ]: 78 : norm_r = length_squared( arrptr( mGrad ), nn );
208 : 78 : norm_g = sqrt( norm_r );
209 : :
210 : 78 : memset( d, 0, 3 * sizeof( double ) * nn );
211 [ + - ]: 78 : memcpy( r, arrptr( mGrad ),
212 : 156 : nn * sizeof( Vector3D ) ); // memcpy(r, mesh->g, 3*sizeof(double)*nn);
213 : 78 : norm_g *= cg_tol;
214 : :
215 [ + - ][ + - ]: 78 : apply_preconditioner( z, r, err );MSQ_ERRRTN( err ); // prec->apply(z, r, prec, mesh);
[ - + ][ # # ]
[ # # ][ - + ]
216 [ + - ]: 78 : negate( p, z, nn );
217 [ + - ]: 78 : rz = inner( r, z, nn );
218 : :
219 : 78 : dMp = 0;
220 : 78 : norm_p = rz;
221 : 78 : norm_d = 0;
222 : :
223 : 78 : cg_iter = 0;
224 [ + + ][ + - ]: 565 : while( ( sqrt( norm_r ) > norm_g ) &&
[ + + ]
225 : : #ifdef DO_STEEP_DESC
226 : : ( norm_d > tr_num_tol ) &&
227 : : #endif
228 : : ( cg_iter < max_cg_iter ) )
229 : : {
230 : 556 : ++cg_iter;
231 : :
232 : 556 : memset( w, 0, 3 * sizeof( double ) * nn );
233 : : // matmul(w, mHess, p); //matmul(w, mesh, p);
234 [ + - ]: 556 : mHess.product( w, p );
235 : :
236 [ + - ]: 556 : kappa = inner( p, w, nn );
237 [ + + ]: 556 : if( kappa <= 0.0 )
238 : : {
239 : 46 : alpha = ( sqrt( dMp * dMp + norm_p * ( radius - norm_d ) ) - dMp ) / norm_p;
240 [ + - ]: 46 : plus_eq_scaled( d, alpha, p, nn );
241 : 46 : break;
242 : : }
243 : :
244 : 510 : alpha = rz / kappa;
245 : :
246 : 510 : norm_dp1 = norm_d + 2.0 * alpha * dMp + alpha * alpha * norm_p;
247 [ + + ]: 510 : if( norm_dp1 >= radius )
248 : : {
249 : 23 : alpha = ( sqrt( dMp * dMp + norm_p * ( radius - norm_d ) ) - dMp ) / norm_p;
250 [ + - ]: 23 : plus_eq_scaled( d, alpha, p, nn );
251 : 23 : break;
252 : : }
253 : :
254 [ + - ]: 487 : plus_eq_scaled( d, alpha, p, nn );
255 [ + - ]: 487 : plus_eq_scaled( r, alpha, w, nn );
256 [ + - ]: 487 : norm_r = length_squared( r, nn );
257 : :
258 [ + - ][ + - ]: 487 : apply_preconditioner( z, r, err );MSQ_ERRRTN( err ); // prec->apply(z, r, prec, mesh);
[ - + ][ # # ]
[ # # ][ - + ]
259 : :
260 : 487 : rzm1 = rz;
261 [ + - ]: 487 : rz = inner( r, z, nn );
262 : 487 : beta = rz / rzm1;
263 [ + - ]: 487 : times_eq_minus( p, beta, z, nn );
264 : :
265 : 487 : dMp = beta * ( dMp + alpha * norm_p );
266 : 487 : norm_p = rz + beta * beta * norm_p;
267 : 487 : norm_d = norm_dp1;
268 : : }
269 : :
270 : : #ifdef DO_STEEP_DESC
271 : : if( norm_d <= tr_num_tol )
272 : : {
273 : : norm_g = length( arrptr( mGrad ), nn );
274 : : double ll = 1.0;
275 : : if( norm_g < tr_num_tol ) break;
276 : : if( norm_g > radius ) ll = radius / nurm_g;
277 : : for( int i = 0; i < nn; ++i )
278 : : d[i] = ll * mGrad[i];
279 : : }
280 : : #endif
281 : :
282 [ + - ][ + - ]: 78 : alpha = inner( arrptr( mGrad ), d, nn ); // inner(mesh->g, d, nn);
283 : :
284 : 78 : memset( p, 0, 3 * sizeof( double ) * nn );
285 : : // matmul(p, mHess, d); //matmul(p, mesh, d);
286 [ + - ]: 78 : mHess.product( p, d );
287 [ + - ]: 78 : beta = 0.5 * inner( p, d, nn );
288 : 78 : kappa = alpha + beta;
289 : :
290 : : /* Put the new point into the locations */
291 [ + - ][ + - ]: 78 : pd.move_free_vertices_constrained( d, nn, 1.0, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
292 : :
293 [ + - ]: 78 : valid = func.evaluate( pd, objn, err );
294 [ + - ][ + + ]: 78 : if( err.error_code() == err.BARRIER_VIOLATED )
295 [ + - ]: 7 : err.clear(); // barrier violated does not represent an actual error here
296 [ + - ][ - + ]: 78 : MSQ_ERRRTN( err );
[ # # ][ # # ]
[ - + ]
297 : :
298 [ + + ]: 78 : if( !valid )
299 : : {
300 : : /* Function not defined at trial point */
301 : 13 : radius *= tr_decr_undef;
302 [ + - ][ + - ]: 13 : pd.set_to_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
303 : 13 : continue;
304 : : }
305 : :
306 [ + + ][ + - ]: 65 : if( ( fabs( kappa ) <= tr_num_tol ) && ( fabs( objn - obj ) <= tr_num_tol ) ) { kappa = 1; }
307 : : else
308 : : {
309 : 64 : kappa = ( objn - obj ) / kappa;
310 : : }
311 : :
312 [ + + ]: 65 : if( kappa < eta_1 )
313 : : {
314 : : /* Iterate is unacceptable */
315 : 25 : radius *= tr_decr_def;
316 [ + - ][ + - ]: 25 : pd.set_to_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
317 : 25 : continue;
318 : : }
319 : :
320 : : /* Iterate is acceptable */
321 [ + + ]: 40 : if( kappa >= eta_2 )
322 : : {
323 : : /* Iterate is a very good step, increase radius */
324 : 13 : radius *= tr_incr;
325 [ - + ]: 13 : if( radius > 1e20 ) { radius = 1e20; }
326 : : }
327 : :
328 [ + - ]: 40 : func.update( pd, obj, mGrad, mHess, err );
329 [ + - ][ + - ]: 40 : compute_preconditioner( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
330 [ + - ][ + - ]: 40 : pd.recreate_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
331 : :
332 : : // checks stopping criterion
333 [ + - ][ + - ]: 40 : term.accumulate_patch( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
334 [ + - ][ + - ]: 40 : term.accumulate_inner( pd, objn, arrptr( mGrad ), err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
335 : : }
336 : : }
337 : :
338 [ + - ][ + - ]: 12 : } // namespace MBMesquite
|
