source: issm/trunk/src/c/parallel/cielo/GradJSearch.c@ 472

Last change on this file since 472 was 472, checked in by Eric.Larour, 16 years ago

New thermal solution. Still missing the thermal elements
File size: 7.7 KB
Line 
1/*
2 * GradJSearch.c:
3 */
4
5#include "../../../config.h"
6
7#if defined(_PARALLEL_) && defined(_HAVE_PETSC_)
8
9#include "../include/cielo.h"
10#include "../modules.h"
11#include "./parallel.h"
12
13#undef __FUNCT__
14#define __FUNCT__ "GradJSearch"
15#undef CLEANUP
16#define CLEANUP GradJSearchLocalCleanup();
17
18void GradJSearchLocalCleanup(void);
19
20
21int GradJSearch(double* search_vector,FemModel* femmodel,int step){
22
23 /*Error management: */
24 int noerr=1;
25 int i,n;
26 int dummy;
27 ParameterInputs* inputs=NULL;
28 int status;
29
30 //status=GoldenSearch(search_vector,femmodel->workspaceparams->J+step,-1,1,femmodel->workspaceparams->tolx,(int)femmodel->workspaceparams->maxiter[step],
31 // femmodel->workspaceparams->fit[step],femmodel->workspaceparams->optscal[step],&objectivefunctionC,femmodel); //do only one dimension search for now.
32
33 status=BrentSearch(search_vector,femmodel->workspaceparams->J+step,-1,1,femmodel->workspaceparams->tolx,(int)femmodel->workspaceparams->maxiter[step],
34 femmodel->workspaceparams->fit[step],femmodel->workspaceparams->optscal[step],&objectivefunctionC,femmodel); //do only one dimension search for now.
35
36 TESTEXIT(noerr);
37
38 return status;
39}
40
41void GradJSearchLocalCleanup(void){
42 return;
43}
44
45int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*,ParameterInputs*),FemModel* femmodel){
46
47 double xc, xd, fc, fd;
48 double oneminustau = 1 - (sqrt(5) - 1) / 2;
49 int iter = 0;
50 ParameterInputs* inputs=NULL;
51 int status;
52
53 inputs=NewParameterInputs();
54
55 xc = xa + oneminustau * (xb - xa);
56 fc = (*f)(&xc,fit,optscal,femmodel,inputs);
57 xd = xb - oneminustau * (xb - xa);
58 fd = (*f)(&xd,fit,optscal,femmodel,inputs);
59 do {
60 iter++;
61 if (fc < fd) {
62 xb = xd;
63 xd = xc;
64 xc = xa + oneminustau * (xb - xa);
65 fd = fc;
66 fc = (*f)(&xc,fit,optscal,femmodel,inputs);
67 }
68 else {
69 xa = xc;
70 xc = xd;
71 xd = xb - oneminustau * (xb - xa);
72 fc = fd;
73 fd = (*f)(&xd,fit,optscal,femmodel,inputs);
74 }
75 _printf_(" iter# %i f(x) %g x %g toler %g/%g iter %i/%i\n",iter,(fc+fd)/2,(xa+xb)/2,fabs(xb-xa),tolerance,iter,maxiter);
76 }
77 while (fabs(xb - xa) > tolerance && iter < maxiter);
78
79 if (fabs(xb-xa)<tolerance)status=0;
80 else status=1;
81
82 /*Assign output pointers: */
83 *psearch_scalar=(xa+xb)/2;
84 *pJ=(fc+fd)/2;
85
86 return status;
87}
88
89int BrentSearch(double* psearch_scalar,double* pJ,double a, double b, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*,ParameterInputs*),FemModel* femmodel){
90
91 /* This routine is optimizing a given function using Brent's method
92 * (Golden or parabolic procedure)*/
93
94 /*optimization variable: */
95 double si;
96 double gold;
97 double intervalgold;
98 double oldintervalgold;
99 double parab_num,parab_den;
100 double distance;
101
102 /*function values: */
103 double fxmax,fxmin,fxbest,fval;
104 double fx,fx1,fx2;
105
106 /*x : */
107 double xmax,xmin,xbest;
108 double x,x1,x2,xm,xval;
109
110 /*tolerances: */
111 double tol1,tol2,seps;
112
113 /*counters: */
114 int iter,goldenflag,loop;
115
116 /*inputs: */
117 int status;
118 ParameterInputs* inputs=NULL;
119
120 /*Recover inputs: */
121 inputs=NewParameterInputs();
122
123 //initialize counter and boundaries
124 iter=0;
125
126 //get the value of the function at the first boundary
127 fxmin = (*f)(&a,fit,optscal,femmodel,inputs);
128
129 //display result
130 _printf_("\n Iteration x f(x) Tolerance Procedure\n\n");
131 _printf_(" %s %12.6g %12.6g %s"," N/A",a,fxmin," N/A boundary\n");
132
133 //get the value of the function at the first boundary b and display result
134 fxmax = (*f)(&b,fit,optscal,femmodel,inputs);
135 _printf_(" %s %12.6g %12.6g %s"," N/A",b,fxmax," N/A boundary\n");
136
137 //initialize the other variables
138 seps=sqrt(DBL_EPSILON); //precision of a double
139 distance=0.0; //new_x=old_x + distance
140 gold=0.5*(3.0-sqrt(5.0)); //gold = 1 - golden ratio
141 intervalgold=0.0; //distance used by Golden procedure
142
143 //Compute initial point
144
145 //1: initialize the value of the 4 x needed (x1,x2,x,xbest)
146 x1=a+gold*(b-a);
147 x2=x1;
148 xbest=x1;
149 x=xbest;
150
151 //2: call the function to be evaluated
152 fxbest = (*f)(&x,fit,optscal,femmodel,inputs);
153 iter=iter+1;
154
155 //3: update the other variables
156 fx1=fxbest;
157 fx2=fxbest;
158 //xm is always in the middle of a and b
159 xm=0.5*(a+b);
160 //update tolerances
161 tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
162 tol2=2.0*tol1;
163
164 //4: print result
165 _printf_(" %5i %12.6g %12.6g %12.6g %s\n",iter,xbest,fxbest,pow(pow(xbest-xm,2),0.5)," initial");
166
167 //Main Loop
168 loop=1;
169 while(loop){
170
171 goldenflag=1;
172
173 // Is a parabolic fit possible ?
174 if (sqrt(pow(intervalgold,2))>tol1){
175
176 // Yes, so fit parabola
177 goldenflag=0;
178 parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
179 parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
180
181 //reverse p if necessary
182 if(parab_den>0.0){
183 parab_num=-parab_num;
184 }
185 parab_den=sqrt(pow(parab_den,2));
186 oldintervalgold=intervalgold;
187 intervalgold=distance;
188
189 // Is the parabola acceptable
190 if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
191 (parab_num>parab_den*(a-xbest)) &&
192 (parab_num<parab_den*(b-xbest))){
193
194 // Yes, parabolic interpolation step
195 distance=parab_num/parab_den;
196 x=xbest+distance;
197
198 // f must not be evaluated too close to min_x or max_x
199 if (((x-a)<tol2) || ((b-x)<tol2)){
200
201 if ((xm-xbest)<0.0){
202 si=-1;
203 }
204 else{
205 si=1;
206 }
207
208 //compute new distance
209 distance=tol1*si;
210 }
211 }
212 else{
213
214 // Not acceptable, must do a golden section step
215 goldenflag=1;
216 }
217 }
218
219 //Golden procedure
220 if(goldenflag){
221
222 // compute the new distance d
223 if(xbest>=xm){
224 intervalgold=a-xbest;
225 }
226 else{
227 intervalgold=b-xbest;
228 }
229 distance=gold*intervalgold;
230 }
231
232 // The function must not be evaluated too close to xbest
233 if(distance<0){
234 si=-1;
235 }
236 else{
237 si=1;
238 }
239 if (sqrt(pow(distance,2))>tol1){
240 x=xbest+si*sqrt(pow(distance,2));
241 }
242 else{
243 x=xbest+si*tol1;
244 }
245
246 //evaluate function on x
247 fx = (*f)(&x,fit,optscal,femmodel,inputs);
248 iter=iter+1;
249
250 // Update a, b, xm, x1, x2, tol1, tol2
251 if (fx<=fxbest){
252 if (x>=xbest){
253 a=xbest;
254 }
255 else{
256 b=xbest;
257 }
258 x1=x2; fx1=fx2;
259 x2=xbest; fx2=fxbest;
260 xbest=x; fxbest=fx;
261 }
262
263 else{ // fx > fxbest
264 if (x < xbest){
265 a=x;
266 }
267 else{
268 b=x;
269 }
270 if ((fx<=fx2) || (x2==xbest)){
271 x1=x2; fx1=fx2;
272 x2=x; fx2=fx;
273 }
274 else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
275 x1=x; fx1=fx;
276 }
277 }
278 xm = 0.5*(a+b);
279 tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
280 tol2=2.0*tol1;
281
282 //print result
283 if (goldenflag){
284 _printf_(" %5i %12.6g %12.6g %12.6g %s\n",iter,x,fx,pow(pow(xbest-xm,2),0.5)," golden");
285 }
286 else{
287 _printf_(" %5i %12.6g %12.6g %12.6g %s\n",iter,x,fx,pow(pow(xbest-xm,2),0.5)," parabolic");
288 }
289
290 //Stop the optimization?
291 if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(b-a))){
292 _printf_("\nOptimization terminated:\nthe current x satisfies the termination criteria using 'tolx' of %g \n", tolerance);
293 loop=0;
294 status=0;
295 }
296 else if (iter>=maxiter){
297 _printf_("\nExiting: Maximum number of iterations has been exceeded - increase 'maxiter'\n");
298 loop=0;
299 status=1;
300 }
301 else{
302 //continue
303 loop=1;
304 }
305 }//end while
306
307 //Now, check that the value on the boundaries are not better than current fxbest
308 if (fxbest>fxmin){
309 xval=xmin;
310 fval=fxmin;
311 }
312 else if (fxbest>fxmax){
313 xval=xmax;
314 fval=fxmax;
315 }
316 else{
317 xval=xbest;
318 fval=fxbest;
319 }
320
321 /*Assign output pointers: */
322 *psearch_scalar=xval;
323 *pJ=fval;
324
325 return status;
326}
327#endif //#if defined(_PARALLEL_) && defined(_HAVE_PETSC_)
328
Note: See TracBrowser for help on using the repository browser.