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 QuasiNewton.cpp
28 : : * \brief Port Todd Munson's quasi-Newton solver to Mesquite
29 : : * \author Jason Kraftcheck (Mesquite Port)
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "QuasiNewton.hpp"
34 : : #include "MsqDebug.hpp"
35 : : #include "MsqError.hpp"
36 : : #include "PatchData.hpp"
37 : :
38 : : namespace MBMesquite
39 : : {
40 : :
41 : : // Force std::vector to release allocated memory
42 : : template < typename T >
43 : 0 : static inline void free_vector( std::vector< T >& v )
44 : : {
45 [ # # ][ # # ]: 0 : std::vector< T > temp;
46 : 0 : temp.swap( v );
47 : 0 : }
48 : :
49 : 0 : std::string QuasiNewton::get_name() const
50 : : {
51 [ # # ]: 0 : return "QuasiNewton";
52 : : }
53 : :
54 : 0 : PatchSet* QuasiNewton::get_patch_set()
55 : : {
56 : 0 : return PatchSetUser::get_patch_set();
57 : : }
58 : :
59 [ + - ][ + - ]: 13 : QuasiNewton::QuasiNewton( ObjectiveFunction* of ) : VertexMover( of ), PatchSetUser( true ), mMemento( 0 ) {}
[ + - ][ + + ]
[ + - ][ + + ]
[ + - ][ + - ]
[ # # # #
# # # # #
# # # # #
# # ]
60 : :
61 [ + - ][ + + ]: 14 : QuasiNewton::~QuasiNewton()
[ + - ][ + + ]
62 : : {
63 [ - + ]: 1 : delete mMemento;
64 : 1 : mMemento = 0;
65 [ - + ]: 1 : }
66 : :
67 : 0 : void QuasiNewton::initialize( PatchData& pd, MsqError& err )
68 : : {
69 [ # # ]: 0 : if( !mMemento )
70 : : {
71 [ # # ][ # # ]: 0 : mMemento = pd.create_vertices_memento( err );MSQ_CHKERR( err );
72 : : }
73 : 0 : }
74 : :
75 : 0 : void QuasiNewton::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
76 : :
77 : 0 : void QuasiNewton::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
78 : :
79 : 0 : void QuasiNewton::cleanup()
80 : : {
81 : : // release memento
82 [ # # ]: 0 : delete mMemento;
83 : 0 : mMemento = 0;
84 : :
85 : : // release coordinates
86 [ # # ]: 0 : for( size_t i = 0; i < ( sizeof( w ) / sizeof( w[0] ) ); ++i )
87 : 0 : free_vector( w[i] );
88 : : // release gradients
89 [ # # ]: 0 : for( size_t i = 0; i < ( sizeof( v ) / sizeof( v[0] ) ); ++i )
90 : 0 : free_vector( v[i] );
91 : :
92 : : // release Hessian memory
93 : 0 : free_vector( mHess );
94 : :
95 : : // release temporary array memory
96 : 0 : free_vector( x );
97 : 0 : free_vector( d );
98 : 0 : }
99 : :
100 : : // Do v += s * x, where v and x are arrays of length n
101 : 0 : static inline void plus_eq_scaled( Vector3D* v, double s, const Vector3D* x, size_t n )
102 : : {
103 : 0 : Vector3D* end = v + n;
104 [ # # ]: 0 : for( ; v != end; ++v, ++x )
105 [ # # ]: 0 : *v += s * *x;
106 : 0 : }
107 : :
108 : 0 : void QuasiNewton::solve( Vector3D* z_arr, const Vector3D* v_arr ) const
109 : : {
110 [ # # ]: 0 : SymMatrix3D pd;
111 : :
112 : 0 : const double small = DBL_EPSILON;
113 : 0 : const size_t nn = mHess.size();
114 [ # # ]: 0 : for( size_t i = 0; i < nn; ++i )
115 : : {
116 : :
117 : : // ensure positive definite: perturb a bit if
118 : : // diagonal values are zero.
119 [ # # ]: 0 : SymMatrix3D dd = mHess[i];
120 [ # # ][ # # ]: 0 : while( fabs( dd[0] ) < small || fabs( dd[3] ) < small || fabs( dd[5] ) < small )
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
121 [ # # ][ # # ]: 0 : dd += small;
122 : :
123 : : // factor
124 [ # # ][ # # ]: 0 : pd[0] = 1.0 / dd[0];
125 [ # # ][ # # ]: 0 : pd[1] = dd[1] * pd[0];
[ # # ]
126 [ # # ][ # # ]: 0 : pd[2] = dd[2] * pd[0];
[ # # ]
127 : :
128 [ # # ][ # # ]: 0 : pd[3] = 1.0 / ( dd[3] - dd[1] * pd[1] );
[ # # ][ # # ]
129 [ # # ][ # # ]: 0 : pd[5] = dd[4] - dd[2] * pd[1];
[ # # ][ # # ]
130 [ # # ][ # # ]: 0 : pd[4] = pd[3] * pd[5];
[ # # ]
131 [ # # ][ # # ]: 0 : pd[5] = 1.0 / ( dd[5] - dd[2] * pd[2] - pd[4] * pd[5] );
[ # # ][ # # ]
[ # # ][ # # ]
132 : :
133 [ # # ][ # # ]: 0 : if( pd[0] <= 0.0 || pd[3] <= 0.0 || pd[5] <= 0.0 )
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
134 : : {
135 [ # # ][ # # ]: 0 : if( dd[0] + dd[3] + dd[5] <= 0 )
[ # # ][ # # ]
136 : : {
137 : : // switch to diagonal
138 [ # # ][ # # ]: 0 : pd[0] = 1.0 / fabs( dd[0] );
139 [ # # ]: 0 : pd[1] = 0.0;
140 [ # # ]: 0 : pd[2] = 0.0;
141 [ # # ][ # # ]: 0 : pd[3] = 1.0 / fabs( dd[3] );
142 [ # # ]: 0 : pd[4] = 0.0;
143 [ # # ][ # # ]: 0 : pd[5] = 1.0 / fabs( dd[5] );
144 : : }
145 : : else
146 : : {
147 : : // diagonal preconditioner
148 [ # # ][ # # ]: 0 : pd[0] = pd[3] = pd[5] = 1.0 / ( dd[0] + dd[3] + dd[5] );
[ # # ][ # # ]
[ # # ][ # # ]
149 [ # # ][ # # ]: 0 : pd[1] = pd[2] = pd[4] = 0.0;
[ # # ]
150 : : }
151 : : }
152 : :
153 : : // solve
154 : 0 : const Vector3D& vv = v_arr[i];
155 : 0 : Vector3D& z = z_arr[i];
156 [ # # ][ # # ]: 0 : z[0] = vv[0];
157 [ # # ][ # # ]: 0 : z[1] = vv[1] - pd[1] * z[0];
[ # # ][ # # ]
158 [ # # ][ # # ]: 0 : z[2] = vv[2] - pd[2] * z[0] - pd[4] * z[1];
[ # # ][ # # ]
[ # # ][ # # ]
159 : :
160 [ # # ][ # # ]: 0 : z[0] *= pd[0];
161 [ # # ][ # # ]: 0 : z[1] *= pd[3];
162 [ # # ][ # # ]: 0 : z[2] *= pd[5];
163 : :
164 [ # # ][ # # ]: 0 : z[1] -= pd[4] * z[2];
[ # # ]
165 [ # # ][ # # ]: 0 : z[0] -= pd[1] * z[1] + pd[2] * z[2];
[ # # ][ # # ]
[ # # ]
166 : : }
167 : 0 : }
168 : :
169 : 0 : void QuasiNewton::optimize_vertex_positions( PatchData& pd, MsqError& err )
170 : : {
171 [ # # ]: 0 : TerminationCriterion& term = *get_inner_termination_criterion();
172 [ # # ]: 0 : OFEvaluator& func = get_objective_function_evaluator();
173 : :
174 : 0 : const double sigma = 1e-4;
175 : 0 : const double beta0 = 0.25;
176 : 0 : const double beta1 = 0.80;
177 : 0 : const double tol1 = 1e-8;
178 : 0 : const double epsilon = 1e-10;
179 : :
180 : : // double norm_r; //, norm_g;
181 : : double alpha, beta;
182 : : double obj, objn;
183 : :
184 : : size_t i;
185 : :
186 : : // Initialize stuff
187 [ # # ]: 0 : const size_t nn = pd.num_free_vertices();
188 : : double a[QNVEC], b[QNVEC], r[QNVEC];
189 [ # # ]: 0 : for( i = 0; i < QNVEC; ++i )
190 : 0 : r[i] = 0;
191 [ # # ]: 0 : for( i = 0; i <= QNVEC; ++i )
192 : : {
193 : 0 : v[i].clear();
194 [ # # ][ # # ]: 0 : v[i].resize( nn, Vector3D( 0.0 ) );
195 : 0 : w[i].clear();
196 [ # # ][ # # ]: 0 : w[i].resize( nn, Vector3D( 0.0 ) );
197 : : }
198 [ # # ]: 0 : d.resize( nn );
199 [ # # ]: 0 : mHess.resize( nn ); // hMesh(mesh);
200 : :
201 [ # # ][ # # ]: 0 : bool valid = func.update( pd, obj, v[QNVEC], mHess, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
202 [ # # ]: 0 : if( !valid )
203 : : {
204 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "Initial objective function is not valid", MsqError::INVALID_MESH );
205 : 0 : return;
206 : : }
207 : :
208 [ # # ][ # # ]: 0 : while( !term.terminate() )
209 : : {
210 [ # # ][ # # ]: 0 : pd.recreate_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
211 [ # # ]: 0 : pd.get_free_vertex_coordinates( w[QNVEC] );
212 : :
213 [ # # ]: 0 : x = v[QNVEC];
214 [ # # ]: 0 : for( i = QNVEC; i--; )
215 : : {
216 [ # # ][ # # ]: 0 : a[i] = r[i] * inner( &( w[i][0] ), arrptr( x ), nn );
[ # # ]
217 [ # # ][ # # ]: 0 : plus_eq_scaled( arrptr( x ), -a[i], &v[i][0], nn );
[ # # ]
218 : : }
219 : :
220 [ # # ][ # # ]: 0 : solve( arrptr( d ), arrptr( x ) );
[ # # ]
221 : :
222 [ # # ]: 0 : for( i = QNVEC; i--; )
223 : : {
224 [ # # ][ # # ]: 0 : b[i] = r[i] * inner( &( v[i][0] ), arrptr( d ), nn );
[ # # ]
225 [ # # ][ # # ]: 0 : plus_eq_scaled( arrptr( d ), a[i] - b[i], &( w[i][0] ), nn );
[ # # ]
226 : : }
227 : :
228 [ # # ][ # # ]: 0 : alpha = -inner( &( v[QNVEC][0] ), arrptr( d ), nn ); /* direction is negated */
[ # # ]
229 [ # # ]: 0 : if( alpha > 0.0 )
230 : : {
231 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "No descent.", MsqError::INVALID_MESH );
232 : 0 : return;
233 : : }
234 : :
235 : 0 : alpha *= sigma;
236 : 0 : beta = 1.0;
237 : :
238 [ # # ][ # # ]: 0 : pd.move_free_vertices_constrained( arrptr( d ), nn, -beta, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
239 [ # # ]: 0 : valid = func.evaluate( pd, objn, v[QNVEC], err );
240 [ # # ][ # # ]: 0 : if( err.error_code() == err.BARRIER_VIOLATED )
241 [ # # ]: 0 : err.clear(); // barrier violated does not represent an actual error here
242 [ # # ][ # # ]: 0 : MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ]
243 [ # # ][ # # ]: 0 : if( !valid || ( obj - objn < -alpha * beta - epsilon && length( &( v[QNVEC][0] ), nn ) >= tol1 ) )
[ # # ][ # # ]
[ # # ][ # # ]
244 : : {
245 : :
246 [ # # ]: 0 : if( !valid ) // function not defined at trial point
247 : 0 : beta *= beta0;
248 : : else // unacceptable iterate
249 : 0 : beta *= beta1;
250 : :
251 : : for( ;; )
252 : : {
253 [ # # ]: 0 : if( beta < tol1 )
254 : : {
255 [ # # ][ # # ]: 0 : pd.set_to_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
256 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "Newton step not good", MsqError::INTERNAL_ERROR );
257 : 0 : return;
258 : : }
259 : :
260 [ # # ][ # # ]: 0 : pd.set_free_vertices_constrained( mMemento, arrptr( d ), nn, -beta, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
261 [ # # ]: 0 : valid = func.evaluate( pd, objn, err );
262 [ # # ][ # # ]: 0 : if( err.error_code() == err.BARRIER_VIOLATED )
263 [ # # ]: 0 : err.clear(); // barrier violated does not represent an actual error here
264 [ # # ][ # # ]: 0 : MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ]
265 [ # # ]: 0 : if( !valid ) // function undefined at trial point
266 : 0 : beta *= beta0;
267 [ # # ]: 0 : else if( obj - objn < -alpha * beta - epsilon ) // unacceptlable iterate
268 : 0 : beta *= beta1;
269 : : else
270 : 0 : break;
271 : : }
272 : : }
273 : :
274 [ # # ]: 0 : for( i = 0; i < QNVEC - 1; ++i )
275 : : {
276 : 0 : r[i] = r[i + 1];
277 : 0 : w[i].swap( w[i + 1] );
278 : 0 : v[i].swap( v[i + 1] );
279 : : }
280 : 0 : w[QNVEC - 1].swap( w[0] );
281 : 0 : v[QNVEC - 1].swap( v[0] );
282 : :
283 [ # # ][ # # ]: 0 : func.update( pd, obj, v[QNVEC], mHess, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
284 : : // norm_r = length_squared( &(v[QNVEC][0]), nn );
285 : : // norm_g = sqrt(norm_r);
286 : :
287 : : // checks stopping criterion
288 [ # # ][ # # ]: 0 : term.accumulate_patch( pd, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
289 [ # # ][ # # ]: 0 : term.accumulate_inner( pd, objn, &v[QNVEC][0], err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
290 : : }
291 : : }
292 : :
293 [ + - ][ + - ]: 4 : } // namespace MBMesquite
|
