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