ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/solve.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (12 years, 4 months ago) by tjod
File size: 13368 byte(s)
Log Message:
test

Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     // mode
4     #define NONE 0
5     #define NEWTON 1
6     #define TNCG 2
7     #define DTNCG 3
8     #define AUTO 4
9     // method
10     #define DIAG 5
11     #define SSOR 6
12     #define ICCG 7
13     #define BLOCK 8
14    
15    
16     #include "pcwin.h"
17     #include "pcmod.h"
18    
19     #include "energies.h"
20     #include "derivs.h"
21     #include "utility.h"
22    
23     //#define min(a,b) (((a) < (b) ? (a) : (b))
24    
25     EXTERN struct t_minvar{
26     double cappa, stpmin, stpmax, angmax;
27     int intmax;
28     } minvar;
29    
30     void solve(int,int,int,int, double *,double *,double *,double *,int *,
31     int *,int *,double *,int *,int *,int *, double (*)());
32    
33    
34     void precond(int,int ,int ,double *,double *,double *,int *,int *,int *,double *);
35     void cholesky(int , double *, double *);
36     void column(int,int *,int *,int *, int *,int *,int *,int *);
37    
38     int stable;
39     double *f, *f_diag;
40     int *c_index, *c_value;
41    
42     void solve(int mode,int method,int negtest,int nvar,double *p,double *x,
43     double *g, double *h, int *h_init, int *h_stop, int *h_index,
44     double *h_diag,int *icycle,int *iter_cg,int *fg_call,double (*fgvalue) ())
45     {
46     int i,j,k, iter,maxiter,maxhess;
47     double *m, *d, *si, *r, *q; //
48     double alpha,beta,delta,sigma;
49     double f_sigma;
50     double *x_sigma;
51     double *g_sigma;
52     double g_norm,g_rms,eps,converge;
53     double hj,gg,dq,rr,dd,rs,rs_new,r_norm;
54     int maxvar = 3*natom;
55    
56     if (natom < 300)
57     maxhess = (3*natom*(3*natom-1))/2;
58     else if (natom < 800)
59     maxhess = (3*natom*(3*natom-1))/3;
60     else
61     maxhess = (3*natom*(3*natom-1))/20;
62    
63    
64     m = dvector(0,maxvar);
65     r = dvector(0,maxvar);
66     si = dvector(0,maxvar);
67     d = dvector(0,maxvar);
68     q = dvector(0,maxvar);
69     x_sigma = dvector(0,maxvar);
70     g_sigma = dvector(0,maxvar);
71     f_diag = dvector(0,maxvar);
72     f = dvector(0,maxhess);
73     c_index = ivector(0,maxhess);
74     c_value = ivector(0,maxhess);
75    
76     if (mode != DTNCG && method != NONE)
77     {
78     for (i=0; i < nvar; i++)
79     m[i] = 1.0 / sqrt(fabs(h_diag[i]));
80    
81     for (i=0; i < nvar; i++)
82     {
83     g[i] = g[i] * m[i];
84     h_diag[i] *= m[i] * m[i];
85     for (j = h_init[i]; j < h_stop[i]; j++)
86     {
87     k = h_index[j];
88     h[j] *= m[i] * m[k];
89     }
90     }
91     }
92    
93     iter = 0;
94     gg = 0.0;
95     for (i=0; i < nvar; i++)
96     {
97     p[i] = 0.0;
98     r[i] = -g[i];
99     gg += g[i]*g[i];
100     }
101     g_norm = sqrt(gg);
102     precond(method, iter,nvar, si, r, h, h_init, h_stop, h_index, h_diag);
103     rs = 0.0;
104     for (i=0; i < nvar; i++)
105     {
106     d[i] = si[i];
107     rs += r[i]*si[i];
108     }
109    
110     if (mode == NEWTON)
111     {
112     eps = 1.0e-10;
113     maxiter = nvar;
114     } else if (mode == TNCG || mode == DTNCG)
115     {
116     delta = 1.0;
117     eps = delta/(double) *icycle;
118     g_rms = g_norm/sqrt((double) nvar);
119     eps = MinFun(eps,g_rms);
120     converge = 1.0;
121     maxiter = (int) (10.0*sqrt((double)nvar));
122     }
123     iter = 1;
124    
125     // evaluate the matrix vector product
126     while (TRUE)
127     {
128     if (mode == TNCG || mode == NEWTON)
129     {
130     for (i=0; i < nvar; i++)
131     q[i] = 0.0;
132     for (i=0; i < nvar; i++)
133     {
134     q[i] += h_diag[i]*d[i];
135     for(j= h_init[i]; j < h_stop[i]; j++)
136     {
137     k = h_index[j];
138     hj = h[j];
139     q[i] += hj*d[k];
140     q[k] += hj*d[i];
141     }
142     }
143     }else if (mode == DTNCG)
144     {
145     dd = 0.0;
146     for (i=0; i < nvar; i++)
147     dd += d[i]*d[i];
148    
149     sigma = 1.0e-7 / sqrt(dd);
150    
151     for (i=0; i < nvar; i++)
152     x_sigma[i] = x[i] + sigma*d[i];
153    
154     fg_call = fg_call + 1;
155     f_sigma = fgvalue (x_sigma,g_sigma);
156    
157     for (i=0; i < nvar; i++)
158     q[i] = (g_sigma[i]-g[i]) / sigma;
159     }
160    
161     // test for direction of negative curvature
162     dq = 0.0;
163     for (i=0; i < nvar; i++)
164     dq += d[i]*q[i];
165    
166     if (negtest)
167     {
168     if (dq <= 0.0)
169     {
170     if (iter == 1)
171     {
172     for (i=0; i < nvar; i++)
173     p[i] = d[i];
174     }
175     goto L_10;
176     }
177     }
178    
179     // test truncated Newton termination criterion
180     alpha = rs / dq;
181     rr = 0.0;
182     for (i=0; i < nvar; i++)
183     {
184     p[i] += alpha*d[i];
185     r[i] -= alpha*q[i];
186     rr += r[i]*r[i];
187     }
188     r_norm = sqrt(rr);
189     if (r_norm/g_norm <= eps)
190     goto L_10;
191    
192     // precondition
193     precond(method,iter,nvar, si, r,h,h_init,h_stop,
194     h_index,h_diag);
195    
196     // update direction
197     rs_new = 0.0;
198     for (i=0; i < nvar; i++)
199     rs_new += r[i]*si[i];
200    
201     beta = rs_new/rs;
202     rs = rs_new;
203    
204     for (i=0; i < nvar; i++)
205     d[i] = si[i] + beta*d[i];
206    
207     // check for overlimit next iteration
208     if (iter >= maxiter)
209     goto L_10;
210    
211     iter++;
212     }
213    
214     // termination
215    
216     L_10:
217     if (mode != DTNCG && method != NONE)
218     {
219     for (i=0; i < nvar; i++)
220     {
221     p[i] *= m[i];
222     g[i] /= m[i];
223     }
224     }
225     iter_cg += iter;
226     free_dvector( m ,0,maxvar);
227     free_dvector( r ,0,maxvar);
228     free_dvector( si ,0,maxvar);
229     free_dvector( d ,0,maxvar);
230     free_dvector( q ,0,maxvar);
231     free_dvector( x_sigma ,0,maxvar);
232     free_dvector( g_sigma ,0,maxvar);
233     free_dvector( f_diag ,0,maxvar);
234     free_dvector( f ,0,maxhess);
235     free_ivector( c_index ,0,maxhess);
236     free_ivector( c_value ,0,maxhess);
237     }
238    
239     void precond(int method, int iter, int nvar, double *s, double *r,
240     double *h, int *h_init, int *h_stop,int *h_index,
241     double *h_diag)
242     {
243     int i,j,k,ii,kk,iii,kkk,nblock,ix,iy,iz,icount;
244     int *c_init, *c_stop; // c_init[maxvar],c_stop[maxvar];
245     double a[6],b[3];
246     double omega,factor,*diag; // diag[maxvar];
247     double maxalpha, alpha, f_i, f_k;
248     int maxvar;
249    
250     maxvar = 3*natom;
251     c_init = ivector(0,maxvar);
252     c_stop = ivector(0,maxvar);
253     diag = dvector(0,maxvar);
254    
255    
256     // no preconditioning
257     if (method == NONE)
258     {
259     for (i=0; i < nvar; i++)
260     s[i] = r[i];
261     }
262     // diagonal
263     if (method == DIAG)
264     {
265     for (i=0; i < nvar; i++)
266     s[i] = r[i]/ fabs(h_diag[i]);
267     }
268     // block diagonal
269     if (method == BLOCK)
270     {
271     nblock = 3;
272     for (i=1; i <= natom; i++)
273     {
274     iz = (i-1)*3+2;
275     iy = (i-1)*3+1;
276     ix = (i-1)*3;
277     a[0] = h_diag[ix];
278     if (h_index[h_init[ix]] == iy)
279     a[1] = h[h_init[ix]];
280     else
281     a[1] = 0.0;
282    
283     if (h_index[h_init[ix]+1] == iz)
284     a[2] = h[h_init[ix]+1];
285     else
286     a[2] = 0.0;
287    
288     a[3] = h_diag[iy];
289     if (h_index[h_init[iy]] == iz)
290     a[4] = h[h_init[iy]];
291     else
292     a[4] = 0.0;
293    
294     a[5] = h_diag[iz];
295     b[0] = r[ix];
296     b[1] = r[iy];
297     b[2] = r[iz];
298     cholesky (nblock,a,b);
299     s[ix] = b[0];
300     s[iy] = b[1];
301     s[iz] = b[2];
302     }
303     }
304    
305     // SSOR
306     if (method == SSOR)
307     {
308     omega = 1.0;
309     factor = 2.0 - omega;
310     for (i=0; i < nvar; i++)
311     {
312     s[i] = r[i] * factor;
313     diag[i] = h_diag[i] / omega;
314     }
315     for (i=0; i < nvar; i++)
316     {
317     s[i] = s[i] / diag[i];
318     for (j=h_init[i]; j< h_stop[i]; j++)
319     {
320     k = h_index[j];
321     s[k] = s[k] - h[j]*s[i];
322     }
323     }
324     for (i=nvar-1; i >= 0; i--)
325     {
326     s[i] = s[i] * diag[i];
327     for (j=h_init[i]; j< h_stop[i]; j++)
328     {
329     k = h_index[j];
330     s[i] = s[i] - h[j]*s[k];
331     }
332     s[i] = s[i] / diag[i];
333     }
334     }
335     // incomplete cholesky preconditioning
336     if (method == ICCG && iter == 0)
337     {
338     column (nvar,h_init,h_stop,h_index,
339     c_init,c_stop,c_index,c_value);
340     stable = TRUE;
341     icount = 0;
342     maxalpha = 2.1;
343     alpha = -0.001;
344     L_10:
345     if (alpha <= 0.0)
346     alpha = alpha + 0.001;
347     else
348     alpha = 2.0 * alpha;
349    
350     if (alpha > maxalpha)
351     {
352     stable = FALSE;
353     fprintf(pcmoutfile,"Precond - Incomplete Cholesky is Unstable\n");
354     } else
355     {
356     factor = 1.0 + alpha;
357     for (i=0; i < nvar; i++)
358     {
359     f_diag[i] = factor * h_diag[i];
360     for(j = c_init[i]; j<= c_stop[i]; j++)
361     {
362     k = c_index[j];
363     f_i = f[c_value[j]];
364     f_diag[i] = f_diag[i] - (f_i*f_i*f_diag[k]);
365     icount++;
366     }
367     if (f_diag[i] <= 0.0) goto L_10;
368     if (f_diag[i] < 1.0e-7) f_diag[i] = 1.0e-7;
369     f_diag[i] = 1.0 / f_diag[i];
370     for( j = h_init[i]; j < h_stop[i]; j++)
371     {
372     k = h_index[j];
373     f[j] = h[j];
374     ii = c_init[i];
375     kk = c_init[k];
376     while (ii <= c_stop[i] && kk <= c_stop[k])
377     {
378     iii = c_index[ii];
379     kkk = c_index[kk];
380     if (iii < kkk)
381     ii = ii + 1;
382     else if (kkk < iii)
383     kk = kk + 1;
384     else
385     {
386     f_i = f[c_value[ii]];
387     f_k = f[c_value[kk]];
388     f[j] = f[j] - (f_i*f_k*f_diag[iii]);
389     ii = ii + 1;
390     kk = kk + 1;
391     icount = icount + 1;
392     }
393     }
394     }
395     }
396     }
397     }
398     // incomplete cholesky preconditioning (solution phase)
399    
400     if (method == ICCG)
401     {
402     if (stable)
403     {
404     for (i=0; i < nvar; i++)
405     s[i] = r[i];
406    
407     for (i=0; i < nvar; i++)
408     {
409     s[i] = s[i] * f_diag[i];
410     for( j = h_init[i]; j < h_stop[i]; j++)
411     {
412     k = h_index[j];
413     s[k] = s[k] - f[j]*s[i];
414     }
415     }
416     for (i = nvar-1; i >= 0; i--)
417     {
418     s[i] = s[i] / f_diag[i];
419     for( j = h_init[i]; j < h_stop[i]; j++)
420     {
421     k = h_index[j];
422     s[i] = s[i] - f[j]*s[k];
423     }
424     s[i] = s[i] * f_diag[i];
425     }
426     } else
427     {
428     for (i=0; i < nvar; i++)
429     s[i] = r[i] / fabs(h_diag[i]);
430     }
431     }
432     free_ivector(c_init ,0,maxvar);
433     free_ivector(c_stop ,0,maxvar);
434     free_dvector(diag ,0,maxvar);
435     }
436    
437     void cholesky(int nvar, double *a, double *b)
438     {
439     int i,j,k,ii,ij,ik,im,jk,jm,ki,kk;
440     double r,sa,t;
441    
442     ii = 0;
443     for (i=0; i < nvar; i++)
444     {
445     im = i - 1;
446     if (i != 0)
447     {
448     ij = i;
449     for( j = 0; j < im; j++)
450     {
451     r = a[ij];
452     if (j != 0)
453     {
454     ik = i;
455     jk = j;
456     jm = j - 1;
457     for ( k = 0; k < jm; k++)
458     {
459     r = r - a[ik]*a[jk];
460     ik = nvar-1 - k + ik;
461     jk = nvar-1 - k + jk;
462     }
463     }
464     a[ij] = r;
465     ij = nvar-1 - j + ij;
466     }
467     }
468     r = a[ii];
469     if (i != 0)
470     {
471     kk = 0;
472     ik = i;
473     for (k = 0; k < im; k++)
474     {
475     sa = a[ik];
476     t = sa * a[kk];
477     a[ik] = t;
478     r = r - sa*t;
479     ik = nvar-1 - k + ik;
480     kk = nvar-1 - k +1 + kk;
481     }
482     }
483     a[ii] = 1.0/ r;
484     ii = nvar-1 - i + 1 + ii;
485     }
486    
487     // solve linear equations; first solve Ly = b for y
488    
489     for (i=0; i < nvar; i++)
490     {
491     if (i != 0)
492     {
493     ik = i;
494     im = i - 1;
495     r = b[i];
496     for (k = 0; k < im; k++)
497     {
498     r = r - b[k]*a[ik];
499     ik = nvar-1 - k + ik;
500     }
501     b[i] = r;
502     }
503     }
504    
505     // now, solve (D)(L transpose)(x) = y for x
506    
507     ii = (nvar*(nvar+1)/2)-1;
508     for ( j = 0; j < nvar; j++)
509     {
510     i = nvar-1 - j;
511     r = b[i] * a[ii];
512     if (j != 0)
513     {
514     im = i + 1;
515     ki = ii + 1;
516     for( k = im; k < nvar; k++)
517     {
518     r = r - a[ki]*b[k];
519     ki++;
520     }
521     }
522     b[i] = r;
523     ii = ii - j - 2;
524     }
525     }