ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/search.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 6557 byte(s)
Log Message:
test

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