mengine/src/tncg.c

#define EXTERN extern


#include "pcwin.h"
#include "pcmod.h"
#include "energies.h"
#include "utility.h"

// mode
#define NONE      0
#define NEWTON    1
#define TNCG      2
#define DTNCG     3
#define AUTO      4
// method
#define DIAG      5
#define SSOR      6
#define ICCG      7
#define BLOCK     8
#define FULL      9
// search results
#define  Success     0
#define  ReSearch    1
#define  WideAngle   2
#define  BadIntpln   3
#define  IntplnErr   4
#define  blank       5

void outminstat(int , double ,double );
void tncg(int,int,int *,double *,double *, double, double (*)(),  double (*)());
double fgvalue(double *, double *);
void search(int,double *,double *,double *,double *,double,double *,int *,double (*)(),int *);
void solve(int,int,int,int, double *,double *,double *,double *,int *,
              int *,int *,double *,int *,int *,int *, double (*)());
void hmatrix(int,double *,double *,int *, int *, int *, double *);
void pcmfout(int);

EXTERN struct t_minvar{
     double cappa, stpmin, stpmax, angmax;
     int   intmax;
     }  minvar;
     
double hesscut;

void tncg(int nvar,int imethod,int *iter,double *x,double *minimum, double grdmin,double (*fgvalue) (),double (*hmatrix) ())
{
    int i,  iter_tn, iter_cg, fg_call, newhess;
    int *h_init, *h_stop, *h_index;         // h_init[maxvar], h_stop[maxvar], h_index[MAXHESS];
    double *x_old, *g, *p, *h_diag;          // x_old[maxvar], g[maxvar], p[maxvar], h_diag[maxvar];
    double  *h;                   // h[MAXHESS];
    double f, angle, rms, x_move, f_move, f_old, g_norm, g_rms;
    int   done, negtest, nextiter, maxiter;
    double fctmin;
    int mode, hmode, status, maxvar, maxhess, method;
    int maxerr,nerr;

    maxerr = 3;
    nerr = 0;
    maxvar = 3*natom;
   if (natom < 500)
      maxhess = (3*natom*(3*natom-1))/2;
   else if (natom < 800)
      maxhess = (3*natom*(3*natom-1))/3;
   else
      maxhess = (3*natom*(3*natom-1))/20;


    h_init = ivector(0,maxvar);
    h_stop = ivector(0,maxvar);
    h_index = ivector(0,maxhess);
    x_old = dvector(0,maxvar);
    g = dvector(0,maxvar);
    p = dvector(0,maxvar);
    h_diag = dvector(0,maxvar);
    h = dvector(0,maxhess);
    
    rms = sqrt((float) nvar)/ sqrt(3.0);

    fctmin = -1000000.0;
    nextiter = 1;
    newhess = 1;
    done = FALSE;
    minvar.cappa = 0.1;
    minvar.stpmin = 1.0e-20;
    minvar.stpmax = 1.5;
    minvar.angmax = 180.0;
    minvar.intmax = 8;
    negtest = TRUE;
    maxiter = 1000;
    
    iter_tn = nextiter - 1;
    maxiter += iter_tn;


    iter_cg = 0;
    fg_call = 1;
    f = fgvalue(x,g);
    f_old = f;
    g_norm = 0.0;
    for (i=0; i < nvar; i++)
    {
         x_old[i] = x[i];
         g_norm += g[i]*g[i];
    }
    g_norm = sqrt(g_norm);
    f_move = 0.5 * minvar.stpmax * g_norm;
    g_rms = g_norm / rms;

    done = FALSE;


    if (g_rms <= grdmin)
    {
        done = TRUE;
        *minimum = f;
    } else if ( f <= fctmin)
    {
        done = TRUE;
        *minimum = f;
    } else if (iter_tn >= maxiter)
    {
        done = TRUE;
        *minimum = f;
    }
    while ( ! done)
    {
        iter_tn++;


//   mode = tcng   method = diag
       if ( g_rms > 3.0)
        mode = TNCG;
       else
        mode = DTNCG;
                   
        hesscut = 0.0;
        if (g_rms > 10.0)
        {
            method = DIAG;
            hesscut = 1.0;
        } else if (nvar < 10)
        {
            method = SSOR;
            hesscut = 1.0;
        }else if (g_rms < 1.0)
        {
            method = ICCG;
            hesscut = .001*nvar;
            if (hesscut > 0.1) hesscut = 0.1;
        } else
        {
            method = ICCG;
            hesscut = 0.001*nvar;
            if (hesscut > 1.0) hesscut = 1.0;
        }
        hmode = FULL;
        
        if ( (iter_tn%newhess) != 0) hmode = NONE;
        if (mode == DTNCG && method == NONE) hmode = NONE;
        if (mode == DTNCG && method == DIAG) hmode = DIAG;

        hmatrix(hmode, x, h, h_init, h_stop, h_index, h_diag);

        solve(mode,method,negtest,nvar, p,x,g,h,h_init,
              h_stop,h_index,h_diag,&iter_tn,&iter_cg,&fg_call,fgvalue);
        search(nvar,&f,g,x,p,f_move,&angle,&fg_call,fgvalue, &status);    
        
         if (status != Success)
         {
                nerr++;
                if (nerr > maxerr)
                {
                  //                    fprintf(pcmlogfile,"Error in tncg :: bad search\n");
                    done = TRUE;
                }
         }
         
         f_move = f_old - f;
         f_old = f;
         x_move = 0.0;
         g_norm = 0.0;
         for(i=0; i < nvar; i++)
         {
            x_move += (x[i]-x_old[i])*(x[i]-x_old[i]);
            x_old[i] = x[i];
            g_norm += g[i]*g[i];
         }
  
         x_move = sqrt(x_move);
         x_move = x_move / rms;
         g_norm = sqrt(g_norm);
         g_rms = g_norm / rms;

         if (iter_tn >= maxiter)
            done = TRUE;

         if (f_move <= 0.0000001)
         {
            done = TRUE;
         }

         if (g_rms <= grdmin)
            done = TRUE;
         else if (f <= fctmin)
            done = TRUE;
             
              
         if (done)
         {
             *minimum = f;
             *iter = iter_tn;
             free_ivector(h_init ,0,maxvar);
             free_ivector(h_stop ,0,maxvar);
             free_ivector(h_index ,0,maxhess);
             free_dvector(x_old ,0,maxvar);
             free_dvector(g ,0,maxvar);
             free_dvector(p ,0,maxvar);
             free_dvector(h_diag ,0,maxvar);
             free_dvector(h ,0,maxhess);

             return;

         }       
    }
    free_ivector(h_init ,0,maxvar);
    free_ivector(h_stop ,0,maxvar);
    free_ivector(h_index ,0,maxhess);
    free_dvector(x_old ,0,maxvar);
    free_dvector(g ,0,maxvar);
    free_dvector(p ,0,maxvar);
    free_dvector(h_diag ,0,maxvar);
    free_dvector(h ,0,maxhess);

    return;   
}
Revision:	93
Committed:	Sat Jan 17 23:24:55 2009 UTC (13 years, 4 months ago) by wdelano
File size:	5960 byte(s)
Log Message:	branch created
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 = (3natom(3*natom-1))/2;
62	else if (natom < 800)
63	maxhess = (3natom(3*natom-1))/3;
64	else
65	maxhess = (3natom(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	}