Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2004 Sandia Corporation and Argonne National
5 : : Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 : : with Sandia Corporation, the U.S. Government retains certain
7 : : rights in 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 : : [email protected], [email protected], [email protected],
24 : : [email protected], [email protected], [email protected]
25 : :
26 : : ***************************************************************** */
27 : : /*!
28 : : \file ConjugateGradient.cpp
29 : : \brief
30 : :
31 : : The Conjugate Gradient class is a concrete vertex mover
32 : : which performs conjugate gradient minimizaiton.
33 : :
34 : : \author Michael Brewer
35 : : \date 2002-06-19
36 : : */
37 : :
38 : : #include "ConjugateGradient.hpp"
39 : : #include <math.h>
40 : : #include "MsqDebug.hpp"
41 : : #include "MsqTimer.hpp"
42 : : //#include "MsqFreeVertexIndexIterator.hpp"
43 : :
44 : : namespace MBMesquite
45 : : {
46 : :
47 : : extern int get_parallel_rank();
48 : : extern int get_parallel_size();
49 : :
50 : 0 : std::string ConjugateGradient::get_name() const
51 : : {
52 [ # # ]: 0 : return "ConjugateGradient";
53 : : }
54 : :
55 : 19 : PatchSet* ConjugateGradient::get_patch_set()
56 : : {
57 : 19 : return PatchSetUser::get_patch_set();
58 : : }
59 : :
60 : 49 : ConjugateGradient::ConjugateGradient( ObjectiveFunction* objective )
61 [ + - ][ + - ]: 49 : : VertexMover( objective ), PatchSetUser( true ), pMemento( NULL ), conjGradDebug( 0 )
[ + - ][ + - ]
62 : : {
63 : 49 : }
64 : :
65 : 4 : ConjugateGradient::ConjugateGradient( ObjectiveFunction* objective, MsqError& err )
66 [ + - ][ + - ]: 4 : : VertexMover( objective ), PatchSetUser( true ), pMemento( NULL ), conjGradDebug( 0 )
[ + - ][ + - ]
67 : : {
68 : : // Michael:: default to global?
69 [ + - ]: 4 : set_debugging_level( 0 );
70 : : // set the default inner termination criterion
71 [ + - ]: 4 : TerminationCriterion* default_crit = get_inner_termination_criterion();
72 [ - + ]: 4 : if( default_crit == NULL )
73 : : {
74 : : MSQ_SETERR( err )
75 : : ( "QualityImprover did not create a default inner "
76 : : "termination criterion.",
77 [ # # ][ # # ]: 0 : MsqError::INVALID_STATE );
78 : 0 : return;
79 : : }
80 : : else
81 : : {
82 [ + - ][ + - ]: 4 : default_crit->add_iteration_limit( 5 );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
83 : : }
84 : : }
85 : :
86 : 107 : ConjugateGradient::~ConjugateGradient()
87 : : {
88 : : // Checks that cleanup() has been called.
89 [ - + ]: 53 : assert( pMemento == NULL );
90 [ - + ]: 54 : }
91 : :
92 : 22 : void ConjugateGradient::initialize( PatchData& pd, MsqError& err )
93 : : {
94 : 22 : if( get_parallel_size() )
95 : : {
96 : : MSQ_DBGOUT( 2 ) << "\nP[" << get_parallel_rank() << "] "
97 : : << "o Performing Conjugate Gradient optimization.\n";
98 : : }
99 : : else
100 : : {
101 : : MSQ_DBGOUT( 2 ) << "\no Performing Conjugate Gradient optimization.\n";
102 : : }
103 : 22 : pMemento = pd.create_vertices_memento( err );
104 : 22 : }
105 : :
106 : 345 : void ConjugateGradient::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
107 : :
108 : : /*!Performs Conjugate gradient minimization on the PatchData, pd.*/
109 : 343 : void ConjugateGradient::optimize_vertex_positions( PatchData& pd, MsqError& err )
110 : : {
111 : : // pd.reorder();
112 : :
113 : : MSQ_FUNCTION_TIMER( "ConjugateGradient::optimize_vertex_positions" );
114 : :
115 [ + - ]: 343 : Timer c_timer;
116 [ + - ]: 343 : size_t num_vert = pd.num_free_vertices();
117 [ - + ]: 343 : if( num_vert < 1 )
118 : : {
119 : : MSQ_DBGOUT( 1 ) << "\nEmpty free vertex list in ConjugateGradient\n";
120 : 343 : return;
121 : : }
122 : : /*
123 : : //zero out arrays
124 : : int zero_loop=0;
125 : : while(zero_loop<arraySize){
126 : : fGrad[zero_loop].set(0,0,0);
127 : : pGrad[zero_loop].set(0,0,0);
128 : : fNewGrad[zero_loop].set(0,0,0);
129 : : ++zero_loop;
130 : : }
131 : : */
132 : :
133 : : // get OF evaluator
134 [ + - ]: 343 : OFEvaluator& objFunc = get_objective_function_evaluator();
135 : :
136 : : size_t ind;
137 : : // Michael cull list: possibly set soft_fixed flags here
138 : :
139 : : // MsqFreeVertexIndexIterator free_iter(pd, err); MSQ_ERRRTN(err);
140 : :
141 : 343 : double f = 0;
142 : : // Michael, this isn't equivalent to CUBIT because we only want to check
143 : : // the objective function value of the 'bad' elements
144 : : // if invalid initial patch set an error.
145 [ + - ]: 343 : bool temp_bool = objFunc.update( pd, f, fGrad, err );
146 [ - + ]: 343 : assert( fGrad.size() == num_vert );
147 [ + - ][ - + ]: 343 : if( MSQ_CHKERR( err ) ) return;
[ # # ][ # # ]
[ - + ]
148 [ - + ]: 343 : if( !temp_bool )
149 : : {
150 : : MSQ_SETERR( err )
151 : : ( "Conjugate Gradient not able to get valid gradient "
152 : : "and function values on intial patch.",
153 [ # # ][ # # ]: 0 : MsqError::INVALID_MESH );
154 : 0 : return;
155 : : }
156 : 343 : double grad_norm = MSQ_MAX_CAP;
157 : :
158 [ - + ]: 343 : if( conjGradDebug > 0 )
159 : : {
160 : : MSQ_PRINT( 2 )( "\nCG's DEGUB LEVEL = %i \n", conjGradDebug );
161 [ # # ][ # # ]: 0 : grad_norm = Linf( arrptr( fGrad ), fGrad.size() );
162 : : MSQ_PRINT( 2 )( "\nCG's FIRST VALUE = %f,grad_norm = %f", f, grad_norm );
163 : : MSQ_PRINT( 2 )( "\n TIME %f", c_timer.since_birth() );
164 : 0 : grad_norm = MSQ_MAX_CAP;
165 : : }
166 : :
167 : : // Initializing pGrad (search direction).
168 [ + - ]: 343 : pGrad.resize( fGrad.size() );
169 [ + + ]: 3775 : for( ind = 0; ind < num_vert; ++ind )
170 [ + - ][ + - ]: 3432 : pGrad[ind] = ( -fGrad[ind] );
[ + - ][ + - ]
171 : :
172 : 343 : int j = 0; // total nb of step size changes ... not used much
173 : 343 : int i = 0; // iteration counter
174 : 343 : unsigned m = 0; //
175 : 343 : double alp = MSQ_MAX_CAP; // alp: scale factor of search direction
176 : : // we know inner_criterion is false because it was checked in
177 : : // loop_over_mesh before being sent here.
178 [ + - ]: 343 : TerminationCriterion* term_crit = get_inner_termination_criterion();
179 : :
180 : : // while ((i<maxIteration && alp>stepBound && grad_norm>normGradientBound)
181 : : // && !inner_criterion){
182 [ + - ][ + + ]: 1673 : while( !term_crit->terminate() )
183 : : {
184 : 1330 : ++i;
185 : : // std::cout<<"\Michael delete i = "<<i;
186 : 1330 : int k = 0;
187 [ + - ]: 1330 : alp = get_step( pd, f, k, err );
188 : 1330 : j += k;
189 : 1330 : if( conjGradDebug > 2 ) { MSQ_PRINT( 2 )( "\n Alp initial, alp = %20.18f", alp ); }
190 : :
191 : : // if alp == 0, revert to steepest descent search direction
192 [ + + ]: 1330 : if( alp == 0 )
193 : : {
194 [ + + ]: 2963 : for( m = 0; m < num_vert; ++m )
195 : : {
196 [ + - ][ + - ]: 2142 : pGrad[m] = ( -fGrad[m] );
[ + - ][ + - ]
197 : : }
198 [ + - ]: 821 : alp = get_step( pd, f, k, err );
199 : 821 : j += k;
200 [ - + ]: 821 : if( conjGradDebug > 1 )
201 : : {
202 : : MSQ_PRINT( 2 )( "\n CG's search direction reset." );
203 : 0 : if( conjGradDebug > 2 ) MSQ_PRINT( 2 )( "\n Alp was zero, alp = %20.18f", alp );
204 : : }
205 : : }
206 [ + + ]: 1330 : if( alp != 0 )
207 : : {
208 [ + - ][ + - ]: 510 : pd.move_free_vertices_constrained( arrptr( pGrad ), num_vert, alp, err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
209 : :
210 [ + - ][ - + ]: 510 : if( !objFunc.update( pd, f, fNewGrad, err ) )
211 : : {
212 : : MSQ_SETERR( err )
213 : : ( "Error inside Conjugate Gradient, vertices moved "
214 : : "making function value invalid.",
215 [ # # ][ # # ]: 0 : MsqError::INVALID_MESH );
216 : 0 : return;
217 : : }
218 [ - + ]: 510 : assert( fNewGrad.size() == (unsigned)num_vert );
219 : :
220 [ - + ]: 510 : if( conjGradDebug > 0 )
221 : : {
222 [ # # ][ # # ]: 0 : grad_norm = Linf( arrptr( fNewGrad ), num_vert );
223 : : MSQ_PRINT( 2 )
224 : : ( "\nCG's VALUE = %f, iter. = %i, grad_norm = %f, alp = %f", f, i, grad_norm, alp );
225 : : MSQ_PRINT( 2 )( "\n TIME %f", c_timer.since_birth() );
226 : : }
227 : 510 : double s11 = 0;
228 : 510 : double s12 = 0;
229 : 510 : double s22 = 0;
230 : : // free_iter.reset();
231 : : // while (free_iter.next()) {
232 : : // m=free_iter.value();
233 [ + + ]: 94393 : for( m = 0; m < num_vert; ++m )
234 : : {
235 [ + - ][ + - ]: 93883 : s11 += fGrad[m] % fGrad[m];
[ + - ]
236 [ + - ][ + - ]: 93883 : s12 += fGrad[m] % fNewGrad[m];
[ + - ]
237 [ + - ][ + - ]: 93883 : s22 += fNewGrad[m] % fNewGrad[m];
[ + - ]
238 : : }
239 : :
240 : : // Steepest Descent (takes 2-3 times as long as P-R)
241 : : // double bet=0;
242 : :
243 : : // Fletcher-Reeves (takes twice as long as P-R)
244 : : // double bet = s22/s11;
245 : :
246 : : // Polack-Ribiere
247 : : double bet;
248 [ + - ][ - + ]: 510 : if( !divide( s22 - s12, s11, bet ) ) return; // gradient is zero
249 : : // free_iter.reset();
250 : : // while (free_iter.next()) {
251 : : // m=free_iter.value();
252 [ + + ]: 94393 : for( m = 0; m < num_vert; ++m )
253 : : {
254 [ + - ][ + - ]: 93883 : pGrad[m] = ( -fNewGrad[m] + ( bet * pGrad[m] ) );
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
255 [ + - ][ + - ]: 93883 : fGrad[m] = fNewGrad[m];
[ + - ]
256 : : }
257 : 510 : if( conjGradDebug > 2 )
258 : : {
259 : : MSQ_PRINT( 2 )
260 : : ( " \nSEARCH DIRECTION INFINITY NORM = %e", Linf( arrptr( fNewGrad ), num_vert ) );
261 : : }
262 : :
263 : : } // end if on alp == 0
264 : :
265 [ + - ][ + - ]: 1330 : term_crit->accumulate_patch( pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
266 [ + - ][ + - ]: 1330 : term_crit->accumulate_inner( pd, f, arrptr( fGrad ), err );MSQ_ERRRTN( err );
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
267 : : } // end while
268 : 343 : if( conjGradDebug > 0 )
269 : : {
270 : : MSQ_PRINT( 2 )( "\nConjugate Gradient complete i=%i ", i );
271 : : MSQ_PRINT( 2 )( "\n- FINAL value = %f, alp=%4.2e grad_norm=%4.2e", f, alp, grad_norm );
272 : : MSQ_PRINT( 2 )( "\n FINAL TIME %f", c_timer.since_birth() );
273 : : }
274 : : }
275 : :
276 : 31 : void ConjugateGradient::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ )
277 : : {
278 : : // cout << "- Executing ConjugateGradient::iteration_complete()\n";
279 : 31 : }
280 : :
281 : 19 : void ConjugateGradient::cleanup()
282 : : {
283 : : // cout << "- Executing ConjugateGradient::iteration_end()\n";
284 : 19 : fGrad.clear();
285 : 19 : pGrad.clear();
286 : 19 : fNewGrad.clear();
287 : : // pMemento->~PatchDataVerticesMemento();
288 [ + - ]: 19 : delete pMemento;
289 : 19 : pMemento = NULL;
290 : 19 : }
291 : :
292 : : //! Computes a distance to move vertices given an initial position and search direction (stored in
293 : : //! data member pGrad).
294 : : /*!Returns alp, the double which scales the search direction vector
295 : : which when added to the old nodal positions yields the new nodal
296 : : positions.*/
297 : : /*!\todo Michael NOTE: ConjugateGradient::get_step's int &j is only
298 : : to remain consisitent with CUBIT for an initial test. It can be
299 : : removed.*/
300 : :
301 : 2151 : double ConjugateGradient::get_step( PatchData& pd, double f0, int& j, MsqError& err )
302 : : {
303 : : // get OF evaluator
304 [ + - ]: 2151 : OFEvaluator& objFunc = get_objective_function_evaluator();
305 : :
306 [ + - ]: 2151 : size_t num_vertices = pd.num_free_vertices();
307 : : // initial guess for alp
308 : 2151 : double alp = 1.0;
309 : 2151 : int jmax = 100;
310 : 2151 : double rho = 0.5;
311 : : // feasible=false implies the mesh is not in the feasible region
312 : 2151 : bool feasible = false;
313 : 2151 : int found = 0;
314 : : // f and fnew hold the objective function value
315 : 2151 : double f = 0;
316 : 2151 : double fnew = 0;
317 : : // Counter to avoid infinitly scaling alp
318 : 2151 : j = 0;
319 : : // save memento
320 [ + - ]: 2151 : pd.recreate_vertices_memento( pMemento, err );
321 : : // if we must check feasiblility
322 : : // while step takes mesh into infeasible region and ...
323 [ + - ][ + + ]: 4320 : while( j < jmax && !feasible && alp > MSQ_MIN )
[ + - ]
324 : : {
325 : 2169 : ++j;
326 [ + - ][ + - ]: 2169 : pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
327 [ + - ]: 2169 : feasible = objFunc.evaluate( pd, f, err );
328 [ + - ][ + + ]: 2169 : if( err.error_code() == err.BARRIER_VIOLATED )
329 [ + - ]: 14 : err.clear(); // barrier violation does not represent an actual error here
330 [ + - ][ - + ]: 2169 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
331 : : // if not feasible, try a smaller alp (take smaller step)
332 [ + + ]: 2169 : if( !feasible ) { alp *= rho; }
333 : : } // end while ...
334 : :
335 : : // if above while ended due to j>=jmax, no valid step was found.
336 [ - + ]: 2151 : if( j >= jmax )
337 : : {
338 : : MSQ_PRINT( 2 )( "\nFeasible Point Not Found" );
339 : 0 : return 0.0;
340 : : }
341 : : // Message::print_info("\nOriginal f %f, first new f = %f, alp = %f",f0,f,alp);
342 : : // if new f is larger than original, our step was too large
343 [ + + ]: 2151 : if( f >= f0 )
344 : : {
345 : 1728 : j = 0;
346 [ + + ][ + + ]: 166090 : while( j < jmax && found == 0 )
347 : : {
348 : 164362 : ++j;
349 : 164362 : alp *= rho;
350 [ + - ][ + - ]: 164362 : pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
351 : : // Get new obj value
352 : : // if patch is now invalid, then the feasible region is convex or
353 : : // we have an error. For now, we assume an error.
354 [ + - ][ - + ]: 164362 : if( !objFunc.evaluate( pd, f, err ) )
355 : : {
356 : : MSQ_SETERR( err )
357 [ # # ][ # # ]: 0 : ( "Non-convex feasiblility region found.", MsqError::INVALID_MESH );
358 : : }
359 [ + - ]: 164362 : pd.set_to_vertices_memento( pMemento, err );
360 [ + - ][ - + ]: 164362 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
361 : : // if our step has now improved the objective function value
362 [ + + ]: 164362 : if( f < f0 ) { found = 1; }
363 : : } // end while j less than jmax
364 : : // Message::print_info("\nj = %d found = %d f = %20.18f f0 = %20.18f\n",j,found,f,f0);
365 : : // if above ended because of j>=jmax, take no step
366 [ + + ]: 1728 : if( found == 0 )
367 : : {
368 : : // Message::print_info("alp = %10.8f, but returning zero\n",alp);
369 : 1641 : alp = 0.0;
370 : 1641 : return alp;
371 : : }
372 : :
373 : 87 : j = 0;
374 : : // while shrinking the step improves the objFunc value further,
375 : : // scale alp down. Return alp, when scaling once more would
376 : : // no longer improve the objFunc value.
377 [ + - ]: 132 : while( j < jmax )
378 : : {
379 : 132 : ++j;
380 : 132 : alp *= rho;
381 : : // step alp in search direction from original positions
382 [ + - ][ + - ]: 132 : pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
383 [ + - ][ - + ]: 132 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
384 : :
385 : : // get new objective function value
386 [ + - ][ - + ]: 132 : if( !objFunc.evaluate( pd, fnew, err ) ) MSQ_SETERR( err )
387 : : ( "Non-convex feasiblility region found while "
388 : : "computing new f.",
389 [ # # ][ # # ]: 0 : MsqError::INVALID_MESH );
390 [ + + ]: 132 : if( fnew < f ) { f = fnew; }
391 : : else
392 : : {
393 : : // Reset the vertices to original position
394 [ + - ]: 87 : pd.set_to_vertices_memento( pMemento, err );
395 [ + - ][ - + ]: 87 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
396 : 87 : alp /= rho;
397 : 87 : return alp;
398 : : }
399 : : }
400 : : // Reset the vertices to original position and return alp
401 [ # # ]: 0 : pd.set_to_vertices_memento( pMemento, err );
402 [ # # ][ # # ]: 0 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ # # ]
403 : 0 : return alp;
404 : : }
405 : : // else our new f was already smaller than our original
406 : : else
407 : : {
408 : 423 : j = 0;
409 : : // check to see how large of step we can take
410 [ + - ][ + + ]: 5845 : while( j < jmax && found == 0 )
411 : : {
412 : 5441 : ++j;
413 : : // scale alp up (rho must be less than 1)
414 : 5441 : alp /= rho;
415 : : // step alp in search direction from original positions
416 [ + - ][ + - ]: 5441 : pd.set_free_vertices_constrained( pMemento, arrptr( pGrad ), num_vertices, alp, err );
417 [ + - ][ - + ]: 5441 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
418 : :
419 [ + - ]: 5441 : feasible = objFunc.evaluate( pd, fnew, err );
420 [ + - ][ + + ]: 5441 : if( err.error_code() == err.BARRIER_VIOLATED )
421 [ + - ]: 18 : err.clear(); // evaluate() error does not represent an actual problem here
422 [ + - ][ - + ]: 5441 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
423 [ + + ]: 5441 : if( !feasible )
424 : : {
425 : 19 : alp *= rho;
426 : : // Reset the vertices to original position and return alp
427 [ + - ]: 19 : pd.set_to_vertices_memento( pMemento, err );
428 [ + - ][ - + ]: 19 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
429 : 19 : return alp;
430 : : }
431 [ + + ]: 5422 : if( fnew < f ) { f = fnew; }
432 : : else
433 : : {
434 : 404 : found = 1;
435 : 404 : alp *= rho;
436 : : }
437 : : }
438 : :
439 : : // Reset the vertices to original position and return alp
440 [ + - ]: 404 : pd.set_to_vertices_memento( pMemento, err );
441 [ + - ][ - + ]: 404 : MSQ_ERRZERO( err );
[ # # ][ # # ]
[ - + ]
442 : 2151 : return alp;
443 : : }
444 : : }
445 : :
446 : : /*!Quadratic one-dimensional line search.*/
447 : : /*
448 : : double ConjugateGradient::get_step(PatchData &pd,double f0,int &j,
449 : : MsqError &err){
450 : : const double CGOLD = 0.3819660;
451 : : const double ZEPS = 1.0e-10;
452 : : int n=pd.num_free_vertices();
453 : : MsqVertex* vertices=pd.get_vertex_array(err);
454 : : double a,b,d,etemp,fb,fu,fv,fw,fx,p,q,r,tol,tol1,tol2,u,v,w,x,xm;
455 : : double e=0.0;
456 : : d=0.0;
457 : : tol=.001;
458 : : int iter, maxiter;
459 : : maxiter=100;
460 : : a=0;
461 : : b=.125;
462 : : int m=0;
463 : : fb=f0-1.0;
464 : : iter=0;
465 : : //find b such that a b 'should' bracket the min
466 : : while (fb<=f0 && iter<maxiter){
467 : : ++iter;
468 : : b*=2.0;
469 : : for(m=0;m<n;++m){
470 : : mCoord[m]=mCoord[m] + (b*pGrad[m]);
471 : : vertices[m]=(mCoord[m]);
472 : : }
473 : : fb=objFunc->evaluate(pd,err);
474 : : }
475 : : iter=0;
476 : : x=w=v=(b/2.0);
477 : : for(m=0;m<n;++m){
478 : : mCoord[m]=mCoord[m] + (x*pGrad[m]);
479 : : vertices[m]=(mCoord[m]);
480 : : }
481 : : fw=fv=fx=objFunc->evaluate(pd,err);
482 : : for(iter=0;iter<maxiter;++iter){
483 : : //Message::print_info("a=%f,b=%f,x=%f,iter=%i\n",a,b,x,iter);
484 : : xm=(a+b)*.5;
485 : : tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
486 : : if(fabs(x-xm)<= (tol2-0.5*(b-a))){
487 : : return x;
488 : : }
489 : : if(fabs(e)>tol1){
490 : : r=(x-w)*(fx-fv);
491 : : q=(x-v)*(fx-fw);
492 : : p=(x-v)*q-(x-w)*r;
493 : : q=2.0*(q-r);
494 : : if(q>0.0)
495 : : p=-p;
496 : : q=fabs(q);
497 : : etemp=e;
498 : : e=d;
499 : : if(fabs(p)>=fabs(0.5*q*etemp)||(p<=q*(a-x))||(p>=q*(b-x))){
500 : : d=CGOLD*(e=(x>=xm?a-x:b-x));
501 : : }
502 : : else{
503 : : d=p/q;
504 : : u=x+d;
505 : : if(u-a<tol2||b-u<tol2)
506 : : {
507 : : if(tol1<0.0)
508 : : d=x-xm;
509 : : else
510 : : d=xm-x;
511 : : }
512 : : }
513 : : }
514 : :
515 : : else{
516 : : d=CGOLD*(e=(x>=xm?a-x:b-x));
517 : : }
518 : : if(tol<0.0)
519 : : u=(fabs(d)>=tol1?x+d:x-d);
520 : : else
521 : : u=(fabs(d)>=tol1?x+d:x+d);
522 : : for(m=0;m<n;++m){
523 : : mCoord[m]=mCoord[m] + (u*pGrad[m]);
524 : : vertices[m]=(mCoord[m]);
525 : : }
526 : : fu=objFunc->evaluate(pd,err);
527 : : if(fu<fx){
528 : : if(u>=x)
529 : : a=x;
530 : : else
531 : : b=x;
532 : : v=w;
533 : : w=x;
534 : : x=u;
535 : : fv=fw;
536 : : fw=fx;
537 : : fx=fu;
538 : : }
539 : : else{
540 : : if(u<x)
541 : : a=u;
542 : : else
543 : : b=u;
544 : : if(fu<=fw||w==x){
545 : : v=w;
546 : : w=u;
547 : : fv=fw;
548 : : fw=fu;
549 : : }
550 : : else if (fu<=fv||v==x||v==w){
551 : : v=u;
552 : : fv=fu;
553 : : }
554 : : }
555 : : }
556 : : for(m=0;m<n;++m){
557 : : vertices[m]=(mCoord[m]);
558 : : }
559 : : //PRINT_WARNING("TOO MANY ITERATIONS IN QUADRATIC LINE SEARCH");
560 : : return x;
561 : : }
562 : : */
563 : :
564 [ + - ][ + - ]: 48 : } // namespace MBMesquite
|
