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