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