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

Line File contents
1 #define EXTERN extern
2
3
4 #include "pcwin.h"
5 #include "pcmod.h"
6 #include "pot.h"
7 #include "energies.h"
8 #include "derivs.h"
9 #include "utility.h"
10 #include "gmmx.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 bounds(void);
33 void outminstat(int , double ,double );
34 void tncg(int,int,int *,double *,double *, double, double (*)(), double (*)());
35 double fgvalue(double *, double *);
36 void search(int,double *,double *,double *,double *,double,double *,int *,double (*)(),int *);
37 void solve(int,int,int,int, double *,double *,double *,double *,int *,
38 int *,int *,double *,int *,int *,int *, double (*)());
39 void hmatrix(int,double *,double *,int *, int *, int *, double *);
40 void inesc(char *);
41 void RefreshScreen(void);
42 void pcmfout(int);
43
44 EXTERN struct t_minvar{
45 double cappa, stpmin, stpmax, angmax;
46 int intmax;
47 } minvar;
48
49 double hesscut;
50
51 void tncg(int nvar,int imethod,int *iter,double *x,double *minimum, double grdmin,double (*fgvalue) (),double (*hmatrix) ())
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
63 maxerr = 3;
64 nerr = 0;
65 maxvar = 3*natom;
66 if (natom < 500)
67 maxhess = (3*natom*(3*natom-1))/2;
68 else if (natom < 800)
69 maxhess = (3*natom*(3*natom-1))/3;
70 else
71 maxhess = (3*natom*(3*natom-1))/20;
72
73
74 h_init = ivector(0,maxvar);
75 h_stop = ivector(0,maxvar);
76 h_index = ivector(0,maxhess);
77 x_old = dvector(0,maxvar);
78 g = dvector(0,maxvar);
79 p = dvector(0,maxvar);
80 h_diag = dvector(0,maxvar);
81 h = dvector(0,maxhess);
82
83 rms = sqrt((float) nvar)/ sqrt(3.0);
84
85 fctmin = -1000000.0;
86 nextiter = 1;
87 newhess = 1;
88 done = FALSE;
89 minvar.cappa = 0.1;
90 minvar.stpmin = 1.0e-20;
91 minvar.stpmax = 1.5;
92 minvar.angmax = 180.0;
93 minvar.intmax = 8;
94 negtest = TRUE;
95 maxiter = 1000;
96
97 iter_tn = nextiter - 1;
98 maxiter += iter_tn;
99
100
101 iter_cg = 0;
102 fg_call = 1;
103 f = fgvalue(x,g);
104 f_old = f;
105 g_norm = 0.0;
106 for (i=0; i < nvar; i++)
107 {
108 x_old[i] = x[i];
109 g_norm += g[i]*g[i];
110 }
111 g_norm = sqrt(g_norm);
112 f_move = 0.5 * minvar.stpmax * g_norm;
113 g_rms = g_norm / rms;
114
115 done = FALSE;
116
117 if (gmmx.run == TRUE && f > gmmx.eminim + 2000.)
118 {
119 done = TRUE;
120 *minimum = f;
121 }
122
123 if (g_rms <= grdmin)
124 {
125 done = TRUE;
126 *minimum = f;
127 } else if ( f <= fctmin)
128 {
129 done = TRUE;
130 *minimum = f;
131 } else if (iter_tn >= maxiter)
132 {
133 done = TRUE;
134 *minimum = f;
135 }
136 while ( ! done)
137 {
138 iter_tn++;
139
140
141 // mode = tcng method = diag
142 if ( g_rms > 3.0)
143 mode = TNCG;
144 else
145 mode = DTNCG;
146
147 hesscut = 0.0;
148 if (g_rms > 10.0)
149 {
150 method = DIAG;
151 hesscut = 1.0;
152 } else if (nvar < 10)
153 {
154 method = SSOR;
155 hesscut = 1.0;
156 }else if (g_rms < 1.0)
157 {
158 method = ICCG;
159 hesscut = .001*nvar;
160 if (hesscut > 0.1) hesscut = 0.1;
161 } else
162 {
163 method = ICCG;
164 hesscut = 0.001*nvar;
165 if (hesscut > 1.0) hesscut = 1.0;
166 }
167 hmode = FULL;
168
169 if ( (iter_tn%newhess) != 0) hmode = NONE;
170 if (mode == DTNCG && method == NONE) hmode = NONE;
171 if (mode == DTNCG && method == DIAG) hmode = DIAG;
172
173 hmatrix(hmode, x, h, h_init, h_stop, h_index, h_diag);
174
175 solve(mode,method,negtest,nvar, p,x,g,h,h_init,
176 h_stop,h_index,h_diag,&iter_tn,&iter_cg,&fg_call,fgvalue);
177 search(nvar,&f,g,x,p,f_move,&angle,&fg_call,fgvalue, &status);
178
179 if (status != Success)
180 {
181 nerr++;
182 if (nerr > maxerr)
183 {
184 // fprintf(pcmoutfile,"Error in tncg :: bad search\n");
185 done = TRUE;
186 }
187 }
188
189 f_move = f_old - f;
190 f_old = f;
191 x_move = 0.0;
192 g_norm = 0.0;
193 for(i=0; i < nvar; i++)
194 {
195 x_move += (x[i]-x_old[i])*(x[i]-x_old[i]);
196 x_old[i] = x[i];
197 g_norm += g[i]*g[i];
198 }
199
200 x_move = sqrt(x_move);
201 x_move = x_move / rms;
202 g_norm = sqrt(g_norm);
203 g_rms = g_norm / rms;
204
205 if (iter_tn >= maxiter)
206 done = TRUE;
207
208 if (f_move <= 0.0000001)
209 {
210 done = TRUE;
211 }
212
213 if (g_rms <= grdmin)
214 done = TRUE;
215 else if (f <= fctmin)
216 done = TRUE;
217
218
219 // print out some results
220 if ( (iter_tn%2) == 0)
221 {
222 if (gmmx.run && gmmx_data.ecut)
223 {
224 if (f_move <= .000001)
225 {
226 done = TRUE;
227 }
228 }
229
230 }
231
232 if (done)
233 {
234 *minimum = f;
235 *iter = iter_tn;
236 free_ivector(h_init ,0,maxvar);
237 free_ivector(h_stop ,0,maxvar);
238 free_ivector(h_index ,0,maxhess);
239 free_dvector(x_old ,0,maxvar);
240 free_dvector(g ,0,maxvar);
241 free_dvector(p ,0,maxvar);
242 free_dvector(h_diag ,0,maxvar);
243 free_dvector(h ,0,maxhess);
244
245 return;
246
247 }
248 }
249 free_ivector(h_init ,0,maxvar);
250 free_ivector(h_stop ,0,maxvar);
251 free_ivector(h_index ,0,maxhess);
252 free_dvector(x_old ,0,maxvar);
253 free_dvector(g ,0,maxvar);
254 free_dvector(p ,0,maxvar);
255 free_dvector(h_diag ,0,maxvar);
256 free_dvector(h ,0,maxhess);
257
258 return;
259 }