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, 10 months ago) by tjod
File size: 6469 byte(s)
Log Message:
test

Line User Rev File contents
1 tjod 3 #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     }