ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/tncg.c
Revision: 122
Committed: Tue Jun 16 22:56:52 2009 UTC (12 years, 7 months ago) by wdelano
File size: 6652 byte(s)
Log Message:
Line File contents
1 #define EXTERN extern
2
3
4 #include "pcwin.h"
5 #include "pcmod.h"
6 #include "utility.h"
7 #include "job_control.h"
8 #include <sys/timeb.h>
9 #include <time.h>
10
11 // mode
12 #define NONE 0
13 #define NEWTON 1
14 #define TNCG 2
15 #define DTNCG 3
16 #define AUTO 4
17 // method
18 #define DIAG 5
19 #define SSOR 6
20 #define ICCG 7
21 #define BLOCK 8
22 #define FULL 9
23 // search results
24 #define Success 0
25 #define ReSearch 1
26 #define WideAngle 2
27 #define BadIntpln 3
28 #define IntplnErr 4
29 #define blank 5
30
31 void outminstat(int , double ,double );
32 void tncg(double,int,int,int *,double *,double *, double, double (*)(), double (*)());
33 double fgvalue(double *, double *);
34 void search(int,double *,double *,double *,double *,double,double *,int *,double (*)(),int *);
35 void solve(int,int,int,int, double *,double *,double *,double *,int *,
36 int *,int *,double *,int *,int *,int *, double (*)());
37 void hmatrix(int,double *,double *,int *, int *, int *, double *);
38
39 EXTERN struct t_minvar{
40 double cappa, stpmin, stpmax, angmax;
41 int intmax;
42 } minvar;
43
44 double hesscut;
45
46 void tncg(double start_time,int nvar,int imethod,int *iter,double *x,double *minimum, double grdmin,double (*fgvalue) (),double (*hmatrix) ())
47 {
48 int i, iter_tn, iter_cg, fg_call, newhess;
49 int *h_init, *h_stop, *h_index; // h_init[maxvar], h_stop[maxvar], h_index[MAXHESS];
50 double *x_old, *g, *p, *h_diag; // x_old[maxvar], g[maxvar], p[maxvar], h_diag[maxvar];
51 double *h; // h[MAXHESS];
52 double f, angle, rms, x_move, f_move, f_old, g_norm, g_rms;
53 int done, negtest, nextiter, maxiter;
54 double fctmin;
55 int mode, hmode, status, maxvar, maxhess, method;
56 int maxerr,nerr;
57 double end_time,old_time,ttime;
58
59 old_time = start_time;
60 maxerr = 3;
61 nerr = 0;
62 maxvar = 3*natom;
63 if (natom < 500)
64 maxhess = (3*natom*(3*natom-1))/2;
65 else if (natom < 800)
66 maxhess = (3*natom*(3*natom-1))/3;
67 else
68 maxhess = (3*natom*(3*natom-1))/20;
69
70
71 h_init = ivector(0,maxvar);
72 h_stop = ivector(0,maxvar);
73 h_index = ivector(0,maxhess);
74 x_old = dvector(0,maxvar);
75 g = dvector(0,maxvar);
76 p = dvector(0,maxvar);
77 h_diag = dvector(0,maxvar);
78 h = dvector(0,maxhess);
79
80 rms = sqrt((float) nvar)/ sqrt(3.0);
81
82 fctmin = -1000000.0;
83 nextiter = 1;
84 newhess = 1;
85 done = FALSE;
86 minvar.cappa = 0.1;
87 minvar.stpmin = 1.0e-20;
88 minvar.stpmax = 1.5;
89 minvar.angmax = 180.0;
90 minvar.intmax = 8;
91 negtest = TRUE;
92 maxiter = 1000;
93
94 iter_tn = nextiter - 1;
95 maxiter += iter_tn;
96
97
98 iter_cg = 0;
99 fg_call = 1;
100 f = fgvalue(x,g);
101 f_old = f;
102 g_norm = 0.0;
103 for (i=0; i < nvar; i++)
104 {
105 x_old[i] = x[i];
106 g_norm += g[i]*g[i];
107 }
108 g_norm = sqrt(g_norm);
109 f_move = 0.5 * minvar.stpmax * g_norm;
110 g_rms = g_norm / rms;
111
112 done = FALSE;
113
114 if (g_rms <= grdmin)
115 {
116 done = TRUE;
117 *minimum = f;
118 } else if ( f <= fctmin)
119 {
120 done = TRUE;
121 *minimum = f;
122 } else if (iter_tn >= maxiter)
123 {
124 done = TRUE;
125 *minimum = f;
126 }
127 while ( ! done)
128 {
129 iter_tn++;
130
131
132 // mode = tcng method = diag
133 if ( g_rms > 3.0)
134 mode = TNCG;
135 else
136 mode = DTNCG;
137
138 hesscut = 0.0;
139 if (g_rms > 10.0)
140 {
141 method = DIAG;
142 hesscut = 1.0;
143 } else if (nvar < 10)
144 {
145 method = SSOR;
146 hesscut = 1.0;
147 }else if (g_rms < 1.0)
148 {
149 method = ICCG;
150 hesscut = .001*nvar;
151 if (hesscut > 0.1) hesscut = 0.1;
152 } else
153 {
154 method = ICCG;
155 hesscut = 0.001*nvar;
156 if (hesscut > 1.0) hesscut = 1.0;
157 }
158 hmode = FULL;
159
160 if ( (iter_tn%newhess) != 0) hmode = NONE;
161 if (mode == DTNCG && method == NONE) hmode = NONE;
162 if (mode == DTNCG && method == DIAG) hmode = DIAG;
163
164 hmatrix(hmode, x, h, h_init, h_stop, h_index, h_diag);
165
166 solve(mode,method,negtest,nvar, p,x,g,h,h_init,
167 h_stop,h_index,h_diag,&iter_tn,&iter_cg,&fg_call,fgvalue);
168 search(nvar,&f,g,x,p,f_move,&angle,&fg_call,fgvalue, &status);
169
170 if (status != Success)
171 {
172 nerr++;
173 if (nerr > maxerr)
174 {
175 // fprintf(pcmlogfile,"Error in tncg :: bad search\n");
176 done = TRUE;
177 }
178 }
179
180 f_move = f_old - f;
181 f_old = f;
182 x_move = 0.0;
183 g_norm = 0.0;
184 for(i=0; i < nvar; i++)
185 {
186 x_move += (x[i]-x_old[i])*(x[i]-x_old[i]);
187 x_old[i] = x[i];
188 g_norm += g[i]*g[i];
189 }
190
191 x_move = sqrt(x_move);
192 x_move = x_move / rms;
193 g_norm = sqrt(g_norm);
194 g_rms = g_norm / rms;
195
196 if (iter_tn >= maxiter)
197 done = TRUE;
198
199 if (f_move <= 0.0000001)
200 {
201 done = TRUE;
202 }
203
204 if (g_rms <= grdmin)
205 done = TRUE;
206 else if (f <= fctmin)
207 done = TRUE;
208
209 // check monitor
210 if (job_control.monitor)
211 {
212 #ifdef WIN32
213 struct _timeb timebuffer;
214 _ftime( &timebuffer );
215 end_time = (timebuffer.time+(timebuffer.millitm/((double)1000.0)));
216 #else
217 struct timeval tv;
218 gettimeofday(&tv,NULL);
219 end_time = tv.tv_sec + (tv.tv_usec/(double)1000000.0);
220 #endif
221 ttime = (end_time-old_time);
222 if ( ttime > 0 && (ttime > job_control.interval) )
223 {
224 // fprintf(pcmlogfile,"Times: %f %f %f\n",start_time,end_time,ttime);
225 old_time = end_time;
226 write_sdf(0);
227 }
228 }
229
230 if (done)
231 {
232 *minimum = f;
233 *iter = iter_tn;
234 free_ivector(h_init ,0,maxvar);
235 free_ivector(h_stop ,0,maxvar);
236 free_ivector(h_index ,0,maxhess);
237 free_dvector(x_old ,0,maxvar);
238 free_dvector(g ,0,maxvar);
239 free_dvector(p ,0,maxvar);
240 free_dvector(h_diag ,0,maxvar);
241 free_dvector(h ,0,maxhess);
242
243 return;
244
245 }
246 }
247 free_ivector(h_init ,0,maxvar);
248 free_ivector(h_stop ,0,maxvar);
249 free_ivector(h_index ,0,maxhess);
250 free_dvector(x_old ,0,maxvar);
251 free_dvector(g ,0,maxvar);
252 free_dvector(p ,0,maxvar);
253 free_dvector(h_diag ,0,maxvar);
254 free_dvector(h ,0,maxhess);
255
256 return;
257 }