ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/search.c
Revision: 85
Committed: Mon Jan 12 15:05:44 2009 UTC (10 years, 9 months ago) by gilbertke
File size: 6547 byte(s)
Log Message:
more cleanup and code removal
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5
6 #include "energies.h"
7 #include "utility.h"
8
9 EXTERN struct t_minvar{
10 double cappa, stpmin, stpmax, angmax;
11 int intmax;
12 } minvar;
13
14 #define Success 0
15 #define ReSearch 1
16 #define WideAngle 2
17 #define BadIntpln 3
18 #define IntplnErr 4
19 #define blank 5
20 #define Failure 6
21
22
23 //#define min(a,b) (((a) < (b) ? (a) : (b))
24
25 double minimize1(double *, double *);
26 void search(int,double *,double *,double *,double *,double,double *,int *,double (*)(),int *);
27
28 void search(int nvar, double *f, double *g, double *x, double *p,
29 double f_move,double *angle,int *ncalls,double (*fgvalue)(),int *status)
30 {
31 int i, intpln;
32 int maxvar;
33 double *s; // s[maxvar];
34 double s_norm,g_norm,cosang;
35 double step,parab,cube,cubstp,sss,ttt;
36 double f_0,f_1,f_a,f_b,f_c,sg_0,sg_1,sg_a,sg_b,sg_c;
37
38 maxvar = 3*natom;
39 s = dvector(0,maxvar);
40
41 if (minvar.cappa == 0.0) minvar.cappa = 0.1;
42 if (minvar.stpmin == 0.0) minvar.stpmin = 1.0e-20;
43 if (minvar.stpmax == 0.0) minvar.stpmax = 2.0;
44 if (minvar.angmax == 0.0) minvar.angmax = 180.0;
45 if (minvar.intmax == 0) minvar.intmax = 5;
46
47 for (i=0; i < nvar; i++)
48 s[i] = p[i];
49
50 g_norm = 0.0;
51 s_norm = 0.0;
52 for (i=0; i < nvar; i++)
53 {
54 g_norm += g[i]*g[i];
55 s_norm += s[i]*s[i];
56 }
57 g_norm = sqrt(g_norm);
58 s_norm = sqrt(s_norm);
59 if (isnan(g_norm) != 0)
60 {
61 fprintf(pcmlogfile,"Found nan in search\n");
62 *status = Failure;
63 free_dvector(s,0,maxvar);
64 return;
65 }
66 f_0 = *f;
67 sg_0 = 0.0;
68 for (i=0; i < nvar; i++)
69 {
70 s[i] = s[i] / s_norm;
71 sg_0 += s[i]*g[i];
72 }
73 cosang = -sg_0 / g_norm;
74 if (cosang < -1.00)
75 cosang = -1.0;
76 if (cosang > 1.00)
77 cosang = 1.0;
78 *angle = radian * acos(cosang);
79 if (*angle > minvar.angmax)
80 {
81 *status = WideAngle;
82 free_dvector(s,0,maxvar);
83 return;
84 }
85 step = 2.0 * fabs(f_move/sg_0);
86 step = MinFun(step,s_norm);
87 if (step > minvar.stpmax) step = minvar.stpmax;
88 if (step < minvar.stpmin) step = minvar.stpmin;
89
90 L_10:
91 intpln = 0;
92 f_b = f_0;
93 sg_b = sg_0;
94
95 L_20:
96 f_a = f_b;
97 sg_a = sg_b;
98 for (i=0; i < nvar; i++)
99 x[i] = x[i] + step*s[i];
100
101 // get new function and gradient
102
103 ncalls++;
104 f_b = fgvalue(x,g);
105 sg_b = 0.0;
106
107 for (i=0; i < nvar; i++)
108 sg_b += s[i]*g[i];
109
110 if (fabs(sg_b/sg_0) <= minvar.cappa && f_b < f_a)
111 {
112 *f = f_b;
113 *status = Success;
114 free_dvector(s,0,maxvar);
115 return;
116 }
117
118 if (sg_a*sg_b < 0.0 || f_a < f_b) goto L_30;
119 if ( fabs(sg_a-sg_b) < 0.0000001)
120 {
121 *f = f_b;
122 *status = Success;
123 free_dvector(s,0,maxvar);
124 return;
125 }
126 step = 2.0 * step;
127 if (sg_b > sg_a)
128 {
129 parab = (f_a-f_b) / (sg_b-sg_a);
130 if (parab > 2.0*step) parab = 2.0 * step;
131 if (parab < 0.5*step) parab = 0.5 * step;
132 step = parab;
133 }
134 if (step > minvar.stpmax) step = minvar.stpmax;
135 goto L_20;
136
137 L_30:
138 intpln++;
139 sss = 3.0*(f_b-f_a)/step - sg_a - sg_b;
140 ttt = sss*sss - sg_a*sg_b;
141 if (ttt < 0.0)
142 {
143 *f = f_b;
144 *status = IntplnErr;
145 free_dvector(s,0,maxvar);
146 return;
147 }
148 ttt = sqrt(ttt);
149 cube = step * (sg_b+ttt+sss)/(sg_b-sg_a+2.0*ttt);
150 if (cube < 0.0 || cube > step)
151 {
152 *f = f_b;
153 *status = IntplnErr;
154 free_dvector(s,0,maxvar);
155 return;
156 }
157 for (i=0; i < nvar; i++)
158 x[i] -= cube*s[i];
159 // get new function and gradient and test
160 ncalls++;
161 f_c = fgvalue(x,g);
162 sg_c = 0.0;
163 for (i=0; i < nvar; i++)
164 sg_c += s[i]*g[i];
165
166 if (fabs(sg_c/sg_0) <= minvar.cappa)
167 {
168 *f = f_c;
169 *status = Success;
170 free_dvector(s,0,maxvar);
171 return;
172 }
173 if (f_c <= f_a || f_c <= f_b)
174 {
175 cubstp = MinFun(fabs(cube),fabs(step-cube));
176 if (cubstp >= minvar.stpmin && intpln < minvar.intmax)
177 {
178 if (sg_a*sg_b < 0.0)
179 {
180 if (sg_a*sg_c < 0.0)
181 {
182 f_b = f_c;
183 sg_b = sg_c;
184 step = step - cube;
185 }else
186 {
187 f_a = f_c;
188 sg_a = sg_c;
189 step = cube;
190 for(i=0; i < nvar; i++)
191 x[i] = x[i] + cube*s[i];
192
193 }
194 } else
195 {
196 if (sg_a*sg_c < 0.0 || f_a <= f_c)
197 {
198 f_b = f_c;
199 sg_b = sg_c;
200 step = step - cube;
201 } else
202 {
203 f_a = f_c;
204 sg_a = sg_c;
205 step = cube;
206 for(i=0; i < nvar; i++)
207 x[i] = x[i] + cube*s[i];
208
209 }
210 }
211 goto L_30;
212 }
213 }
214
215 // failed interpolation
216 if (f_a < f_b && f_a < f_c)
217 {
218 f_1 = f_a;
219 sg_1 = sg_a;
220 for (i=0; i < nvar; i++)
221 x[i] += (cube-step)*s[i];
222 } else if (f_b < f_a && f_b < f_c)
223 {
224 f_1 = f_b;
225 sg_1 = sg_b;
226 for (i=0; i < nvar; i++)
227 x[i] += cube*s[i];
228 } else
229 {
230 f_1 = f_c;
231 sg_1 = sg_c;
232 }
233 // try to restart
234 if (f_1 > f_0)
235 {
236 ncalls++;
237 *f = fgvalue(x,g);
238 *status = IntplnErr;
239 free_dvector(s,0,maxvar);
240 return;
241 }
242 f_0 = f_1;
243 sg_0 = sg_1;
244 if (sg_1 > 0.0)
245 {
246 for (i=0; i < nvar; i++)
247 s[i] = -s[i];
248 sg_0 = -sg_1;
249 }
250 if (cube > (step-cube) )
251 step = cube/10.0;
252 else
253 step = (step-cube)/10.0;
254
255 if (step < minvar.stpmin) step = minvar.stpmin;
256
257 if (*status == ReSearch)
258 {
259 ncalls++;
260 *f = fgvalue(x,g);
261 *status = BadIntpln;
262 free_dvector(s,0,maxvar);
263 return;
264 } else
265 {
266 *status = ReSearch;
267 goto L_10;
268 }
269 }