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