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 (10 years, 9 months ago) by gilbertke
File size: 5960 byte(s)
Log Message:
more cleanup and code removal
Line File contents
1 #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 // fprintf(pcmlogfile,"Error in tncg :: bad search\n");
174 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
208 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 }