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:	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 = (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	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			}